11institutetext: Franz Achleitner 22institutetext: TU Wien, Institute of
Analysis and Scientific Computing, Wiedner Hauptstr. 8-10, A-1040 Wien,
Austria, 22email<EMAIL_ADDRESS>33institutetext: Anton Arnold
44institutetext: TU Wien, Institute of Analysis and Scientific Computing,
Wiedner Hauptstr. 8-10, A-1040 Wien, Austria, 44email:
<EMAIL_ADDRESS>55institutetext: Ansgar Jüngel 66institutetext: TU
Wien, Institute of Analysis and Scientific Computing, Wiedner Hauptstr. 8-10,
A-1040 Wien, Austria, 66email<EMAIL_ADDRESS>
# Necessary and sufficient conditions for strong stability of explicit
Runge–Kutta methods

Franz Achleitner  Anton Arnold  Ansgar Jüngel

###### Abstract

Strong stability is a property of time integration schemes for ODEs that
preserve temporal monotonicity of solutions in arbitrary (inner product)
norms. It is proved that explicit Runge–Kutta schemes of order $p\in
4{\mathbb{N}}$ with $s=p$ stages for linear autonomous ODE systems are not
strongly stable, closing an open stability question from [Z. Sun and C.-W.
Shu, SIAM J. Numer. Anal. 57 (2019), 1158–1182]. Furthermore, for explicit
Runge–Kutta methods of order $p\in{\mathbb{N}}$ and $s>p$ stages, we prove
several sufficient as well as necessary conditions for strong stability. These
conditions involve both the stability function and the hypocoercivity index of
the ODE system matrix. This index is a structural property combining the
Hermitian and skew-Hermitian part of the system matrix.

###### Keywords:

Strong stability, linear ordinary differential equations, hypocoercivity
(index), stability function; MSC 65L06, 65L20.

##  0.1 Introduction

Explicit Runge–Kutta methods are standard tools in the numerical solution of
ordinary differential equations and semi-discrete approximations of partial
differential equations. Of particular interest are strongly stable schemes for
which the norm of its numerical solution is nonincreasing in time,
guaranteeing that the numerical error in each time step is not amplified
during time integration. Surprisingly, the characterization of strongly stable
explicit schemes is not complete. For instance, while it is well-known
SunShu19  that explicit Runge–Kutta methods for linear systems of order $p$
are _not_ strongly stable if $p\in 4{\mathbb{N}}_{0}+1$ or $p\in
4{\mathbb{N}}_{0}+2$ and that they are strongly stable for $p\in
4{\mathbb{N}}_{0}+3$, the case $p\in 4{\mathbb{N}}$ is still open.

In this work, we show that the proof of strong stability can be reduced to
testing strong stability with respect to (w.r.t.) the set of asymptotically
stable, semi-dissipative matrices (see below for the definitions) and w.r.t.
the set of purely imaginary scalars. We deduce that explicit Runge–Kutta
methods for linear systems with order $p\in 4{\mathbb{N}}$ and $p=s$ stages
are _not_ strongly stable, thus filling the gap left open in (SunShu19, ,
Theorem 4.2). Furthermore, we prove that strong stability w.r.t.
asymptotically stable, semi-dissipative matrices holds if the so-called
hypocoercivity index of the system matrix is small enough.

We consider linear time-invariant systems of ordinary differential equations,

$\frac{\,\operatorname{d}\\!u}{\,\,\operatorname{d}\\!t}=\mathbf{L}u,\quad
t>0,\quad u(0)=u^{0}\in{\mathbb{C}}^{n},$ (1)

with Lyapunov stable matrices $\mathbf{L}\in{\mathbb{C}}^{n\times n}$; see
(Be18, , Definition 15.9.1), (HadChe08, , p.172). In this paper, we use
several notions of stability for matrices $\mathbf{L}\in{\mathbb{C}}^{n\times
n}$:

  * •

$\mathbf{L}$ is _Lyapunov stable_ if all eigenvalues have nonpositive real
part and those eigenvalues with vanishing real part are non-defective (i.e.,
the algebraic and geometric multiplicities coincide) (Be18, , Definition
15.9.1).

  * •

$\mathbf{L}$ is _asymptotically stable_ if all eigenvalues have negative real
part (Be18, , Definition 15.9.1).

  * •

$\mathbf{L}$ is _dissipative_ (resp.  _semi-dissipative_) if its Hermitian
part $(\mathbf{L}+\mathbf{L}^{*})/2$ is negative definite (resp. negative
semi-definite) (Be18, , Definition 4.1.1).

Dissipative matrices are asymptotically stable, and semi-dissipative matrices
are Lyapunov stable. Lyapunov stability characterizes all system matrices
$\mathbf{L}\in{\mathbb{C}}^{n\times n}$ such that $u\equiv 0$ is a Lyapunov
stable solution to (1).

We use the following notation: The conjugate transpose of a matrix
$\mathbf{L}\in{\mathbb{C}}^{n\times n}$ is denoted by $\mathbf{L}^{*}$.
Matrices $\mathbf{L}\in{\mathbb{C}}^{n\times n}$ have a unique decomposition
$\mathbf{L}=\mathbf{L}_{H}+\mathbf{L}_{S}$ into its Hermitian part
$\mathbf{L}_{H}:=(\mathbf{L}+\mathbf{L}^{*})/2$ and skew-Hermitian part
$\mathbf{L}_{S}:=(\mathbf{L}-\mathbf{L}^{*})/2$. We write $\mathbb{H}_{n}$ for
the set of all Hermitian matrices in ${\mathbb{C}}^{n\times n}$. Positive
definiteness (resp. semi-definiteness) of Hermitian matrices
$\mathbf{P}\in\mathbb{H}_{n}$ is denoted by $\mathbf{P}>0$ (resp.
$\mathbf{P}\geq 0$), and
$\mathbb{H}_{n}^{>}:=\\{\mathbf{H}\in\mathbb{H}_{n}:\mathbf{H}>0\\}$ is the
set of positive definite Hermitian matrices in ${\mathbb{C}}^{n\times n}$.
Moreover, $\langle\cdot,\cdot\rangle$ is the standard inner product in
${\mathbb{C}}^{n}$ with norm $\|\cdot\|=\sqrt{\langle\cdot,\cdot\rangle}$, and
$\|\cdot\|_{2}$ denotes the spectral norm in ${\mathbb{C}}^{n\times n}$.

It is known that a matrix $\mathbf{L}\in{\mathbb{C}}^{n\times n}$ is Lyapunov
stable if and only if there exists a positive definite Hermitian matrix
$\mathbf{P}\in\mathbb{H}_{n}^{>}$ such that

$\mathbf{L}^{*}\mathbf{P}+\mathbf{P}\mathbf{L}\leq 0;$ (2)

see, e.g., (Str75, , Corollary to Theorem 5), (HadChe08, , Theorem 3.18).
Under this condition, the solution $u(t)$ of (1) is nonincreasing in the norm
$\|\cdot\|_{\mathbf{P}}:=\sqrt{\langle\cdot,\mathbf{P}\cdot\rangle}$:

$\frac{\,\operatorname{d}\\!}{\,\,\operatorname{d}\\!t}\|u(t)\|_{\mathbf{P}}^{2}=\frac{\,\operatorname{d}\\!}{\,\,\operatorname{d}\\!t}\langle
u,\mathbf{P}u\rangle=\langle\mathbf{L}u,\mathbf{P}u\rangle+\langle
u,\mathbf{P}\mathbf{L}u\rangle=\langle
u,(\mathbf{L}^{*}\mathbf{P}+\mathbf{P}\mathbf{L})u\rangle\leq 0.$ (3)

Often, the aim of numerical schemes is to reproduce this property on the
discrete level, i.e. to show that a certain norm of an approximation $u^{k}$
of $u(k\tau)$, where $\tau>0$ is the uniform time step, is nonincreasing in
$k$. For this, we consider the explicit Runge–Kutta methods

$u^{k}=u^{k-1}+\tau\sum_{i=1}^{s}b_{i}K_{i}^{k},\quad
K_{i}^{k}=\mathbf{L}\bigg{(}u^{k-1}+\tau\sum_{j=1}^{i-1}a_{ij}K_{j}^{k}\bigg{)},\quad
i=1,\ldots,s,$ (4)

where $b_{i}\in{\mathbb{C}}$ are the weights, $a_{ij}\in{\mathbb{C}}$ are the
coefficients of the Runge–Kutta matrix, and $s\in{\mathbb{N}}$ is the number
of stages. Stability of explicit Runge–Kutta schemes can be expected only
under a restriction on the time step. In this article, we give necessary and
sufficient conditions such that explicit Runge–Kutta methods preserve the
monotonicity property (3) in the following sense.

######  Definition 1 (Strongly stable)

  1. (a)

The Runge–Kutta scheme (4) is strongly stable if for all matrix dimensions
$n\in{\mathbb{N}}$, for all Lyapunov stable matrices
$\mathbf{L}\in{\mathbb{C}}^{n\times n}$, and for all
$\mathbf{P}\in\mathbb{H}_{n}^{>}$ such that (2) holds, the numerical solution
to (1) satisfies $\|u^{1}\|_{\mathbf{P}}\leq\|u^{0}\|_{\mathbf{P}}$ for all
initial data $u^{0}\in{\mathbb{C}}^{n}$ and sufficiently small time steps.

  2. (b)

For practical reasons, we call the strongly stable Runge–Kutta scheme (4)
_strongly stable w.r.t. the set $\mathcal{L}$ of Lyapunov stable matrices_ of
any matrix dimension $n\in{\mathbb{N}}$,

$\mathcal{L}:=\bigcup_{n\in{\mathbb{N}}}\\{\mathbf{L}\in{\mathbb{C}}^{n\times
n}:\exists\mathbf{P}\in\mathbb{H}_{n}^{>}\text{ such
that~{}\eqref{ineq:Lyapunov} holds}\\}.$

  3. (c)

Similarly, the Runge–Kutta scheme (4) is called _strongly stable w.r.t. a
subset $\mathcal{L}_{0}$ of $\mathcal{L}$_ if the above property holds for all
$\mathbf{L}\in\mathcal{L}_{0}$.

###### Remark 1

  1. (a)

For explicit Runge–Kutta methods, strong stability w.r.t. $\mathcal{L}$ is
equivalent to strong stability w.r.t. _semi-dissipative_ matrices
$\mathbf{L}\in\mathcal{L}_{SD}$ with

$\mathcal{L}_{SD}:=\bigcup_{n\in{\mathbb{N}}}\\{\mathbf{L}\in{\mathbb{C}}^{n\times
n}:\mathbf{L}+\mathbf{L}^{*}\leq 0\\}\subset\mathcal{L};$

see the first step in the proof of Theorem 0.2.1 below.

  2. (b)

Strong stability of explicit Runge–Kutta schemes w.r.t. $\mathcal{L}$ (as
studied here and in (Tad02, , §4), SunShu19 ) is “more restrictive” than
strong stability w.r.t. (uniformly) coercive system matrices $\mathbf{L}$ as
studied in LeTa98 , (Tad02, , §3).

The goal of this paper is twofold: On the one hand, we give novel necessary
conditions for strong stability of explicit Runge–Kutta schemes. On the other
hand, we derive sufficient conditions for strong stability in a weaker sense:
In contrast to Definition 1, we then require the monotonicity property only in
some weighted norm $\|\cdot\|_{\mathbf{P}}$ with
$\mathbf{P}\in\mathbb{H}_{n}^{>}$.

######  Definition 2 (Strong stability in weak form)

  1. (a)

The Runge–Kutta scheme (4) is strongly stable in weak form if for all
$n\in{\mathbb{N}}$ and for all Lyapunov stable matrices
$\mathbf{L}\in{\mathbb{C}}^{n\times n}$, there exists a matrix
$\mathbf{P}\in\mathbb{H}_{n}^{>}$ satisfying (2), such that the numerical
solution (4) to (1) satisfies
$\|u^{1}\|_{\mathbf{P}}\leq\|u^{0}\|_{\mathbf{P}}$ for all initial data
$u^{0}\in{\mathbb{C}}^{n}$ and sufficiently small time steps.

  2. (b)

Similarly, the Runge–Kutta scheme (4) is called _strongly stable in weak form
w.r.t. a subset $\mathcal{L}_{0}$ of $\mathcal{L}$_ if the above property
holds for all $\mathbf{L}\in\mathcal{L}_{0}$.

This notion of _strong stability in weak sense_ is inspired by an example in
(LeTa98, , §3.5); see Example 1.

Most of the practical interest in strongly stable Runge–Kutta methods has been
for explicit schemes, e.g. in the integration of hyperbolic conservation laws
Shu88 ; ShuOsh88 . For implicit or implicit-explicit strongly stable
Runge–Kutta methods, we refer to, e.g., CGGS17 ; FeSp04 ; GST01 ; Hig04 ;
SunWeiWu22 . Furthermore, nonlinear problems are considered in, e.g., Ran21 ;
RaKe20 . Note that a Runge–Kutta method of linear order $p$ (i.e., the order
for linear systems) may possess a lower order when applied to nonlinear
systems.

Strong stability of the (implicit) midpoint rule can be seen as a corollary to
Remark 1(a) and (AAM23, , Lemma 47(ii)).

Let us summarize the main results of this paper:

  * •

An explicit Runge-Kutta scheme is strongly stable if and only if it is
strongly stable w.r.t. the set $\mathcal{L}_{AS}$ of asymptotically stable,
semi-dissipative matrices and w.r.t. the set $\mathcal{L}_{IM}$ of purely
imaginary scalars, see Theorem 0.2.1.

  * •

Explicit Runge–Kutta schemes of order $p\in 4{\mathbb{N}}$ and $s=p$ stages
are not strongly stable. While these Runge–Kutta schemes are strongly stable
w.r.t. $\mathcal{L}_{IM}$, they fail to be strongly stable w.r.t.
$\mathcal{L}_{AS}$: see Corollary 1. For this analysis, we derive novel
necessary conditions for strong stability of explicit Runge–Kutta schemes
using the precise short-time asymptotics for the spectral norm of the matrix
exponential as derived in (AAC22, , Theorem 2.7).

  * •

Concerning strong stability w.r.t. $\mathcal{L}_{AS}$, we show a remarkable
connection between the order $p\in{\mathbb{N}}$ of the Runge–Kutta scheme and
the hypocoercivity (HC) index of asymptotically stable, semi-dissipative
matrices; see Theorem 0.2.2. Roughly speaking, the HC-index of a matrix
$\mathbf{L}$ describes the structural complexity of the intertwining of the
Hermitian part $\mathbf{L}_{H}$ and skew-Hermitian part $\mathbf{L}_{S}$; see
Definition 3.

  * •

Finally, we show that each explicit Runge–Kutta method of order
$p\in{\mathbb{N}}$ is strongly stable in weak form if and only if it is
locally stable on the imaginary axis as defined in (KrSc92, , Definition 2.1);
see Theorem 0.2.4.

##  0.2 Main results

Before stating our main theorems, we recall that the Runge–Kutta scheme (4)
for ODE systems (1) can be written as $u^{k}=R(\tau\mathbf{L})u^{k-1}$, where
$R(z)$ is the stability function (HaWa96, , Definition 2.1). For linear time-
invariant systems, $s=p$ stages are sufficient to obtain a scheme of order $p$
(HaWa96, , § IV.2). Generally, if an explicit Runge–Kutta method is of order
$p$, its stability function can be written as

$R(z)=\sum_{j=0}^{p}\frac{z^{j}}{j!}+\sum_{j=p+1}^{s}c_{j}\frac{z^{j}}{j!},\quad
z\in{\mathbb{C}},\ c_{p+1}\neq 1,$ (5)

with some constants $c_{p+1},\ldots,c_{s}\in{\mathbb{R}}$ (HaWa96, , Theorem
IV.2.2) and $c_{j}=0$ for $j\geq s+1$. For an explicit Runge–Kutta scheme with
$s=p$, this expression reduces to the truncated exponential
$R(z)=\sum_{j=0}^{p}z^{j}/j!$; see, e.g., (HaWa96, , (2.12)).

Our first main result shows that, to establish strong stability for an
explicit Runge–Kutta scheme, it is (necessary and) sufficient to prove strong
stability w.r.t. two distinct subsets of Lyapunov stable matrices.

###### Theorem 0.2.1

Consider an explicit Runge–Kutta method with stability function (4). Then the
Runge–Kutta scheme is strongly stable if and only if it is strongly stable
w.r.t.

  1. (a)

the set $\mathcal{L}_{AS}$ of asymptotically stable, semi-dissipative
matrices, and

  2. (b)

the set $\mathcal{L}_{IM}$ of purely imaginary scalars $\mathbf{L}\in{\rm
i}{\mathbb{R}}\subset{\mathbb{C}}$.

This means that the Runge–Kutta scheme is strongly stable if and only if for
all linear systems (1) with $\mathbf{L}$ either in $\mathcal{L}_{AS}$ or
$\mathcal{L}_{IM}$, for sufficiently small time steps, the numerical solution
to (1) satisfies $\|u^{1}\|\leq\|u^{0}\|$ for all initial data
$u^{0}\in{\mathbb{C}}^{n}$.

The (strong) stability of Runge–Kutta schemes for scalar differential
equations is well studied; see, e.g., KrSc92 . In particular, strong stability
w.r.t. $\mathcal{L}_{IM}$ is equivalent to _local stability on the imaginary
axis_ (i.e., there exists $Z>0$ such that $|R(z)|\leq 1$ for all $z\in{\rm
i}{\mathbb{R}}$ with $|z|\leq Z$). For example, the Runge–Kutta method with
$p=s=4$ is locally stable on the imaginary axis.

Therefore, we shall focus here on strong stability w.r.t. asymptotically
stable, semi-dissipative matrices $\mathbf{L}\in\mathcal{L}_{AS}$. We shall
show that a Runge–Kutta method with $p=s=4$ fails to be strongly stable for
certain asymptotically stable, semi-dissipative matrices
$\mathbf{L}\in\mathcal{L}_{AS}$. We also discuss the (counter-)examples given
in (LeTa98, , §3.5) and (SunShu17, , Proposition 1.1).

Our second main result is concerned with strong stability w.r.t.
asymptotically stable, semi-dissipative matrices $\mathcal{L}_{AS}$. We show a
remarkable connection between the order $p\in{\mathbb{N}}$ of the Runge–Kutta
scheme and the hypocoercivity index (HC-index) of asymptotically stable, semi-
dissipative matrices $\mathbf{L}\in\mathcal{L}_{AS}$, which is defined as
follows; see (AAC22, , Definition 3.1), (AAM21, , Definition 3).

###### Definition 3

Let $\mathbf{L}\in{\mathbb{C}}^{n\times n}$ be semi-dissipative. The
_hypocoercivity index (HC-index) $m_{HC}=m_{HC}(\mathbf{L})$ of the (semi-
dissipative) matrix $\mathbf{L}$_ is defined as the smallest integer
$m\in{\mathbb{N}}_{0}$ (if it exists) such that

$T_{m}:=\sum_{j=0}^{m}\mathbf{L}_{S}^{j}\mathbf{L}_{H}(\mathbf{L}_{S}^{*})^{j}<0.$
(6)

###### Remark 2

Originally, the HC-index has been defined for accretive matrices
$\mathbf{C}\in{\mathbb{C}}^{n\times n}$ (i.e., matrices with positive semi-
definite Hermitian part $\mathbf{C}_{H}\geq 0$) in AAC18 . For practical
reasons, we present the hypocoercivity theory here for semi-dissipative
matrices (using the fact that for a semi-dissipative matrix
$\mathbf{L}\in{\mathbb{C}}^{n\times n}$, the matrix $\mathbf{C}:=-\mathbf{L}$
is accretive).

By definition, a semi-dissipative matrix $\mathbf{L}$ is dissipative if and
only if $m_{HC}(\mathbf{L})=0$. We recall from (AAC18, , Lemma 2.4) that every
matrix $\mathbf{L}\in\mathcal{L}_{AS}$ has a finite HC-index. For the next
theorem, we introduce the following subsets of $\mathcal{L}_{AS}$:

$\mathcal{L}_{AS}^{m}:=\\{\mathbf{L}\in\mathcal{L}_{AS}\>:\
{m_{HC}}(\mathbf{L})\leq m\\},\quad m\in{\mathbb{N}}_{0}.$

###### Theorem 0.2.2

Consider explicit Runge–Kutta schemes of order $p\in{\mathbb{N}}$ for (1) with
stability function (5). Then the following results for strong stability w.r.t.
asymptotically stable, semi-dissipative matrices
$\mathbf{L}\in\mathcal{L}_{AS}$ hold:

  1. (a)

All explicit Runge–Kutta schemes of order $p\in{\mathbb{N}}$ are strongly
stable w.r.t. $\mathcal{L}_{AS}^{m}$, if $m\in{\mathbb{N}}_{0}$ satisfies
$2m+1\leq p$.

  2. (b)

An explicit Runge–Kutta scheme of order $p\in{\mathbb{N}}$ is _not_ strongly
stable w.r.t. $\mathcal{L}_{AS}$ if the coefficient $c_{p+1}$ in (5) satisfies

$\displaystyle(-1)^{(p+1)/2}(1-c_{p+1})<0$ for $p$ odd; (7a)
$\displaystyle\text{or}\qquad 1+(-1)^{p/2}(c_{p+1}-1)\tbinom{p}{p/2}<0$ for
$p$ even. (7b)

  3. (c)

In case (7) holds and the matrix $\mathbf{L}\in\mathcal{L}_{AS}$ satisfies
$2{m_{HC}}(\mathbf{L})+1>p$, the considered Runge–Kutta scheme is _not_
strongly stable w.r.t. the singleton $\\{\mathbf{L}\\}$. This means that it is
strongly stable w.r.t. $\mathcal{L}_{AS}^{m}$ but _not_ strongly stable w.r.t.
$\mathcal{L}_{AS}\setminus\mathcal{L}_{AS}^{m}$ for any $m\in{\mathbb{N}}$
with $2m+1\leq p$.

The conditions (7) on the coefficients of the stability functions directly
yield the following result (which is also included in Table 1 and illustrated
in the examples in Section 0.3).

###### Corollary 1

Explicit Runge–Kutta schemes of order $p\in 4{\mathbb{N}}$ or $p\in
4{\mathbb{N}}_{0}+1$ with $s=p$ stages (and hence with $c_{p+1}=0$) are _not_
strongly stable.

Next, we shall combine two sufficient criteria for explicit Runge–Kutta
schemes to _fail_ strong stability: the criterion to fail strong stability
w.r.t. $\mathcal{L}_{AS}$ (given in Theorem 0.2.2(b)), and w.r.t.
$\mathcal{L}_{IM}$ (given in Theorem 0.3.1 — and taken from the literature:
(KrSc92, , Theorem 3.1)).

###### Theorem 0.2.3

Explicit Runge–Kutta schemes of order $p\in{\mathbb{N}}$ for (1) are _not_
strongly stable if

$\displaystyle(-1)^{(p+1)/2}(1-c_{p+1})<0,$ for $p$ odd; (8a)
$\displaystyle\begin{rcases}(-1)^{p/2}\big{(}c_{p+2}-(p+2)c_{p+1}+(p+1)\big{)}<0\\\
\mbox{or }\quad 1+(-1)^{p/2}(c_{p+1}-1)\tbinom{p}{p/2}<0\end{rcases}$ for $p$
even. (8b)

###### Remark 3

The negation of Theorem 0.2.3 represents necessary conditions for strong
stability. As such it complements the sufficient condition for strongly
stable, explicit Runge–Kutta schemes of even order $p\in 2{\mathbb{N}}$ given
in (SunShu19, , Theorem 4.5). More precisely, the negation of the first
condition of (8b) coincides with the first condition in (SunShu19, , Theorem
4.5), except of the (non-)strict relation sign. The negation of the second
condition from (8b) and its analog from (SunShu19, , Theorem 4.5) both give a
bound on $c_{p+1}$, but they leave a gap between the sufficient and necessary
conditions. We notice that condition (8a) is sufficient and necessary in
SunShu19 . In Table 1, we summarize the results for explicit Runge–Kutta
schemes of order $p\in{\mathbb{N}}$ and $s=p$ stages.

$p$ | $\mathcal{L}_{AS}$ | $\mathcal{L}_{IM}$ | $\mathcal{L}$
---|---|---|---
$4{\mathbb{N}}_{0}+1$ | No | No | No
$4{\mathbb{N}}_{0}+2$ | ? | No | No
$4{\mathbb{N}}_{0}+3$ | Yes | Yes | Yes
$4{\mathbb{N}}$ | No | Yes | No
Table 1: shows whether or not explicit Runge–Kutta methods of order
$p\in{\mathbb{N}}$ with $s=p$ stages are strongly stable w.r.t. asymptotically
stable, semi-dissipative matrices $\mathcal{L}_{AS}$, purely imaginary scalars
$\mathcal{L}_{IM}$, and Lyapunov stable matrices $\mathcal{L}$, respectively;
new results are in bold face. Theorem 0.2.3 shows that explicit Runge–Kutta
schemes for linear systems of order $p\in 4{\mathbb{N}}$ and $s=p$ stages are
not strongly stable. This closes the open case in (SunShu19, , Section 4.2).

Our third main result is concerned with sufficient conditions for strong
stability in weak form.

###### Theorem 0.2.4

An explicit Runge–Kutta scheme of order $p\in{\mathbb{N}}$ is strongly stable
in weak form if and only if it is strongly stable w.r.t. $\mathcal{L}_{IM}$
(i.e., if and only if it is locally stable on the imaginary axis).

Due to Theorem 0.2.1, an explicit Runge–Kutta method is strongly stable if and
only if it is strongly stable w.r.t. to the sets $\mathcal{L}_{AS}$ and
$\mathcal{L}_{IM}$. For asymptotically stable, semi-dissipative matrices
$\mathbf{L}\in\mathcal{L}_{AS}$, using a suitable, modified inner product norm
ensures monotonicity/strong stability.

This paper is organized as follows. We prove Theorem 0.2.1 in Section 0.3,
while Section 0.4 is concerned with the proof of Theorem 0.2.2. Finally, we
present the proof of Theorem 0.2.4 in Section 0.5.

##  0.3 Proof of Theorem 0.2.1

The proof is based on a reduction of the matrix $\mathbf{L}$ to block diagonal
form. First, we recall some preparatory results and then detail the reduction
strategy.

Preparations. It is well-known that a similarity transformation of a Lyapunov
stable matrix yields a semi-dissipative transformed matrix.

###### Lemma 1

Let $\mathbf{L}\in{\mathbb{C}}^{n\times n}$ be Lyapunov stable. Then there
exists $\mathbf{P}\in\mathbb{H}_{n}^{>}$ such that
$\widehat{\mathbf{L}}:=\mathbf{P}^{1/2}\mathbf{L}\mathbf{P}^{-1/2}$ is semi-
dissipative and $\widehat{u}:=\mathbf{P}^{1/2}u$ transforms (1) into the semi-
dissipative ODE system
$\mathrm{d}\widehat{u}/\mathrm{d}t=\widehat{\mathbf{L}}\widehat{u}$. Here,
$\mathbf{P}^{-1/2}:=(\mathbf{P}^{-1})^{1/2}$.

###### Proof

Since $\mathbf{L}$ is Lyapunov stable, there exists
$\mathbf{P}\in\mathbb{H}_{n}^{>}$ such that
$\mathbf{L}^{*}\mathbf{P}+\mathbf{P}\mathbf{L}\leq 0$. A congruence
transformation with the Hermitian matrix $\mathbf{P}^{-1/2}$ yields

$0\geq\mathbf{P}^{-1/2}(\mathbf{L}^{*}\mathbf{P}+\mathbf{P}\mathbf{L})\mathbf{P}^{-1/2}=\mathbf{P}^{-1/2}\mathbf{L}^{*}\mathbf{P}^{1/2}+\mathbf{P}^{1/2}\mathbf{L}\mathbf{P}^{-1/2}=2(\mathbf{P}^{1/2}\mathbf{L}\mathbf{P}^{-1/2})_{H},$

proving that $\widehat{\mathbf{L}}$ is semi-dissipative. Finally, we multiply
(1) from the left by $\mathbf{P}^{1/2}$ to find that
$\mathrm{d}\widehat{u}/\mathrm{d}t=\mathbf{P}^{1/2}\mathbf{L}u=(\mathbf{P}^{1/2}\mathbf{L}\mathbf{P}^{-1/2})\widehat{u}=\widehat{\mathbf{L}}\widehat{u}$.
$\Box$

We can characterize semi-dissipative matrices that are not asymptotically
stable; see (MehMS16, , Lemma 3.1), (AAC18, , Lemma 2.4 with Prop. 1 (B2),
(B4)).

###### Proposition 1

Let $\mathbf{L}\in{\mathbb{C}}^{n\times n}$ be semi-dissipative. Then
$\mathbf{L}$ has an eigenvalue on the imaginary axis (and is hence _not_
asymptotically stable) if and only if $\mathbf{L}_{H}v=0$ for some eigenvector
$v$ of $\mathbf{L}_{S}$.

Note that, due to the assumption on $\mathbf{L}$, purely imaginary eigenvalues
of semi-dissipative matrices are necessarily semi-simple, i.e., their
algebraic and geometric multiplicities coincide; see MehMW18 ; MehMW20 .

_Unitary congruence transformations_
$\mathbf{L}\mapsto\mathbf{U}\mathbf{L}\mathbf{U}^{*}$ with unitary matrices
$\mathbf{U}$, preserve (asymptotic) stability, semi-dissipativity, and the HC-
index (cf. (6)). For semi-dissipative matrices
$\mathbf{L}\in{\mathbb{C}}^{n\times n}$, the pair
$(\mathbf{J},\mathbf{R}):=(\mathbf{L}_{S},\mathbf{L}_{H})$ can be transformed
into a (so-called) _staircase form_ by using unitary congruence
transformations (AAM23, , Lemma 57).

######  Lemma 2 (Staircase form for $(\mathbf{J},\mathbf{R})$)

Let $\mathbf{J}\in{\mathbb{C}}^{n\times n}$ be a skew-Hermitian matrix and
$\mathbf{R}\in{\mathbb{C}}^{n\times n}$ be a nonzero Hermitian matrix. Then
there exists a unitary matrix $\mathbf{V}\in{\mathbb{C}}^{n\times n}$ such
that $\mathbf{V}\mathbf{J}\mathbf{V}^{*}$ and
$\mathbf{V}\mathbf{R}\mathbf{V}^{*}$ are block tridiagonal matrices of the
form

$\displaystyle\mathbf{V}\mathbf{J}\mathbf{V}^{*}\\!\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\\!\\!\begin{array}[]{l}\left[\begin{array}[]{ccccccc|c}\mathbf{J}_{1,1}&-\mathbf{J}_{2,1}^{*}&&&\cdots&&0&0\\\
\mathbf{J}_{2,1}&\mathbf{J}_{2,2}&-\mathbf{J}_{3,2}^{*}&&&&&\\\
&\ddots&\ddots&\ddots&&&\vdots&\\\
&&\mathbf{J}_{k,k-1}&\mathbf{J}_{k,k}&-\mathbf{J}_{k+1,k}^{*}&&&\vdots\\\
\vdots&&&\ddots&\ddots&\ddots&&\\\
&&&&\mathbf{J}_{r-2,r-3}&\mathbf{J}_{r-2,r-2}&-\mathbf{J}_{r-1,r-2}^{*}&\\\
0&&\cdots&&&\mathbf{J}_{r-1,r-2}&\mathbf{J}_{r-1,r-1}&0\\\ \hline\cr
0&&&\cdots&&&0&\mathbf{J}_{r,r}\end{array}\right]\begin{array}[]{c}n_{1}\\\
n_{2}\\\\[2.0pt] \vdots\\\\[4.0pt] n_{k}\\\\[2.0pt] \vdots\\\\[2.0pt]
n_{r-2}\\\ n_{r-1}\\\\[1.0pt] n_{r}\end{array}\\\ \quad\hskip 5.0ptn_{1}\hskip
130.0ptn_{r-2}\hskip 22.0ptn_{r-1}\hskip 18.0ptn_{r},\end{array}$ (9s)
$\displaystyle\mathbf{V}\mathbf{R}\mathbf{V}^{*}\\!\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\\!\\!\begin{array}[]{l}\left[\begin{array}[]{cccc|c}\mathbf{R}_{1}&0&\cdots&0&0\\\
0&0&&\vdots&\vdots\\\ \vdots&&\ddots&\vdots&\vdots\\\ 0&\cdots&\cdots&0&0\\\
\hline\cr
0&\cdots&\cdots&0&0\end{array}\right]\begin{array}[]{c}n_{1}\\\\[4.0pt]
n_{2}\\\ \vdots\\\ n_{r-1}\\\ n_{r}\end{array}\\\ \quad\hskip 1.0ptn_{1}\hskip
8.0ptn_{2}\hskip 4.0pt\cdots\hskip 4.0ptn_{r-1}\hskip 6.0ptn_{r},\end{array}$
(9af)

where $n_{1}\geq n_{2}\geq\cdots\geq n_{r-1}>0$, $n_{r}\geq 0$, and
$\mathbf{R}_{1}\in{\mathbb{C}}^{n_{1},n_{1}}$ is nonsingular.

If $\mathbf{R}$ is nonsingular, then $r=2$ and $n_{2}=0$. For example,
$\mathbf{V}=\mathbf{I}$, $\mathbf{J}_{1,1}=\mathbf{J}$, and
$\mathbf{R}_{1}=\mathbf{R}$ is an admissible choice.

If $\mathbf{R}$ is singular, then $r\geq 3$ and the matrices
$\mathbf{J}_{i,i-1}$, $i=2,\ldots,r-1$, in the subdiagonal have full row rank
and are of the form

$\mathbf{J}_{i,i-1}=\begin{bmatrix}\Sigma_{i,i-1}&0\end{bmatrix},\quad
i=2,\ldots,r-1,$

with nonsingular matrices $\Sigma_{i,i-1}\in{\mathbb{C}}^{n_{i}\times n_{i}}$,
and moreover $\Sigma_{r-1,r-2}$ is a real-valued diagonal matrix.

For semi-dissipative matrices $\mathbf{L}\in{\mathbb{C}}^{n\times n}$, the
staircase form of $(\mathbf{J},\mathbf{R})=(\mathbf{L}_{S},\mathbf{L}_{H})$
can be used to derive a unitary congruence transformation to a block diagonal
form.

###### Lemma 3

Let $\mathbf{L}\in{\mathbb{C}}^{n\times n}$ be semi-dissipative. Then there
exists a unitary matrix $\mathbf{V}\in{\mathbb{C}}^{n\times n}$ such that
$\mathbf{V}\mathbf{L}\mathbf{V}^{*}$ is a block diagonal matrix of the form

$\mathbf{V}\ \mathbf{L}\
\mathbf{V}^{*}=\left[\begin{array}[]{c|c}\mathbf{L}_{1}&0\\\ \hline\cr
0&\mathbf{L}_{2}\end{array}\right],$ (10)

where $\mathbf{L}_{1}\in{\mathbb{C}}^{\tilde{n}_{1}\times\tilde{n}_{1}}$,
$\mathbf{L}_{2}\in{\mathbb{C}}^{\tilde{n}_{2}\times\tilde{n}_{2}}$ with
$\tilde{n}_{1},\tilde{n}_{2}\in{\mathbb{N}}_{0}$ such that
$\tilde{n}_{2}=n-\tilde{n}_{1}$. Moreover, the semi-dissipative matrix
$\mathbf{L}_{1}\in{\mathbb{C}}^{\tilde{n}_{1}\times\tilde{n}_{1}}$ is
asymptotically stable, and the semi-dissipative matrix
$\mathbf{L}_{2}\in{\mathbb{C}}^{\tilde{n}_{2}\times\tilde{n}_{2}}$ is skew-
Hermitian.

In the block-diagonal form (10), not both matrices
$\mathbf{L}_{1},\mathbf{L}_{2}$ have to be present. Indeed, if $\mathbf{L}$ is
skew-Hermitian then $\tilde{n}_{1}=0$, $\tilde{n}_{2}=n$ and
$\mathbf{L}_{2}=\mathbf{L}$, and if $\mathbf{L}$ is asymptotically stable then
$\tilde{n}_{1}=n$, $\tilde{n}_{2}=0$ and $\mathbf{L}_{1}=\mathbf{L}$.

######  Proof (of Lemma 3)

If the Hermitian part $\mathbf{L}_{H}$ of $\mathbf{L}$ is the null matrix then
we set $\mathbf{V}=\mathbf{I}$, $\mathbf{L}_{2}=\mathbf{L}=\mathbf{L}_{S}$,
$\tilde{n}_{1}=0$, and $\tilde{n}_{2}=n$.

If the Hermitian part $\mathbf{L}_{H}$ is not the null matrix, we deduce from
Lemma 2 for the pair
$(\mathbf{J},\mathbf{R}):=(\mathbf{L}_{S},\mathbf{L}_{H})$ the existence of a
unitary matrix $\mathbf{V}\in{\mathbb{C}}^{n\times n}$ such that
$\mathbf{V}\mathbf{L}\mathbf{V}^{*}=\mathbf{V}\mathbf{J}\mathbf{V}^{*}+\mathbf{V}\mathbf{R}\mathbf{V}^{*}$
given in (9) holds with $\mathbf{R}_{1}<0$ (although the matrix
$\mathbf{J}_{r,r}$ may not be present). If the matrix $\mathbf{J}_{r,r}$ is
present, then $\mathbf{L}_{2}=\mathbf{J}_{r,r}$, $\tilde{n}_{2}=n_{r}$,
$\tilde{n}_{1}=n-n_{r}>0$, otherwise $\tilde{n}_{2}=0$ and $\tilde{n}_{1}=n$.
If the matrix $\mathbf{L}_{2}$ is present, it is skew-Hermitian, since it is a
(non-leading) principal minor of the skew-Hermitian matrix
$\mathbf{V}\mathbf{L}_{S}\mathbf{V}^{*}$.

Because of $\mathbf{L}_{H}\neq 0$, the matrix $\mathbf{L}_{1}\neq 0$ is the
leading principal minor of $\mathbf{V}\mathbf{L}\mathbf{V}^{*}$ given as

$\displaystyle\mathbf{L}_{1}:=\begin{array}[]{l}\left[\begin{array}[]{ccccccc}\mathbf{J}_{1,1}+\mathbf{R}_{1}&-\mathbf{J}_{2,1}^{*}&&&\cdots&&0\\\
\mathbf{J}_{2,1}&\mathbf{J}_{2,2}&-\mathbf{J}_{3,2}^{*}&&&&\\\
&\ddots&\ddots&\ddots&&&\vdots\\\
&&\mathbf{J}_{k,k-1}&\mathbf{J}_{k,k}&-\mathbf{J}_{k+1,k}^{*}&&\\\
\vdots&&&\ddots&\ddots&\ddots&\\\
&&&&\mathbf{J}_{r-2,r-3}&\mathbf{J}_{r-2,r-2}&-\mathbf{J}_{r-1,r-2}^{*}\\\
0&&\cdots&&&\mathbf{J}_{r-1,r-2}&\mathbf{J}_{r-1,r-1}\end{array}\right]\begin{array}[]{c}n_{1}\\\
n_{2}\\\\[2.0pt] \vdots\\\\[2.0pt] n_{k}\\\\[1.0pt] \vdots\\\ n_{r-2}\\\
n_{r-1}\end{array},\\\ \quad\hskip 5.0ptn_{1}\hskip 140.0ptn_{r-2}\hskip
25.0ptn_{r-1}\end{array}$

and $\mathbf{L}_{1}\in{\mathbb{C}}^{\tilde{n}_{1}\times\tilde{n}_{1}}$ is
semi-dissipative.

Finally, we apply Proposition 1 in ${\mathbb{C}}^{\tilde{n}_{1}}$ to deduce
that $\mathbf{L}_{1}$ is asymptotically stable. If this were not true, there
would exist a vector $\tilde{v}_{1}=(0,v_{2},...,v_{r-1})^{T}$ in the kernel
of $(\mathbf{L}_{1})_{H}=\operatorname{diag}(\mathbf{R}_{1},0,...,0)$ (see
(9af)) that is an eigenvector of $(\mathbf{L}_{1})_{S}$. The latter has the
block tridiagonal shape in (9s). Since $\mathbf{J}^{*}_{i,i-1}$,
$i=2,...,r-1$, has full column rank, we obtain iteratively $\tilde{v}_{1}=0$,
which cannot be an eigenvector of $(\mathbf{L}_{1})_{S}$. $\Box$

###### Remark 4

A semi-dissipative matrix $\mathbf{L}=\mathbf{L}_{S}+\mathbf{L}_{H}$ is
asymptotically stable if and only if $n_{r}=0$ holds in the staircase form
(9af) of $(\mathbf{J},\mathbf{R})=(\mathbf{L}_{S},\mathbf{L}_{H})$; see the
final step in the proof of Lemma 3. If this is the case, the HC-index
satisfies ${m_{HC}}(\mathbf{L})=r-2$; see (AAM21, , Lemma 4).

Reduction steps. We prove the structural result in Theorem 0.2.1 by performing
several reduction steps.

_Step 1. Reduction to semi-dissipative matrices._ First, we show that strong
stability of explicit Runge–Kutta schemes w.r.t. Lyapunov stable matrices
$\mathbf{L}\in\mathcal{L}$ is equivalent to strong stability w.r.t. semi-
dissipative matrices $\widehat{\mathbf{L}}\in\mathcal{L}_{SD}$ (as already
stated in Remark 1). To this end, we consider a Lyapunov stable matrix
$\mathbf{L}\in{\mathbb{C}}^{n\times n}$ and $\mathbf{P}\in\mathbb{H}_{n}^{>}$
satisfying (2). By Lemma 1,
$\widehat{\mathbf{L}}:=\mathbf{P}^{1/2}\mathbf{L}\mathbf{P}^{-1/2}$ is semi-
dissipative. We define $\widehat{u}^{0}:=\mathbf{P}^{1/2}u^{0}$,
$u^{1}:=R(\tau\mathbf{L})u^{0}$, and
$\widehat{u}^{1}:=R(\tau\widehat{\mathbf{L}})\widehat{u}^{0}$. Then
$\|u^{0}\|_{\mathbf{P}}=\|\mathbf{P}^{1/2}u^{0}\|=\|\widehat{u}^{0}\|$ and

$\displaystyle\|u^{1}\|_{\mathbf{P}}$
$\displaystyle=\|R(\tau\mathbf{L})u^{0}\|_{\mathbf{P}}=\|\mathbf{P}^{1/2}R(\tau\mathbf{L})u^{0}\|$
$\displaystyle=\|R(\tau\mathbf{P}^{1/2}\mathbf{L}\mathbf{P}^{-1/2})\mathbf{P}^{1/2}u^{0}\|=\|R(\tau\widehat{\mathbf{L}})\widehat{u}^{0}\|=\|\widehat{u}^{1}\|.$

Hence, $\|u^{1}\|_{\mathbf{P}}\leq\|u^{0}\|_{\mathbf{P}}$ if and only if
$\|\widehat{u}^{1}\|\leq\|\widehat{u}^{0}\|$.

_Step 2. Reduction to semi-dissipative matrices in block diagonal form._
Considering the unitary congruence transformation from Lemma 3, we note that
$\|\mathbf{V}u^{j}\|=\|u^{j}\|$, with the unitary matrix $\mathbf{V}$ from
(10). Therefore, an explicit Runge–Kutta scheme is strongly stable if and only
if it is strongly stable w.r.t. semi-dissipative matrices $\mathbf{L}$ of the
form (10).

_Step 3. Reduction to two distinct subsets of semi-dissipative matrices._ Due
to the block diagonal structure of $\mathbf{L}$ in (10), an explicit
Runge–Kutta scheme is strongly stable if and only if it is strongly stable
w.r.t. two distinct subsets of semi-dissipative matrices, namely

  1. (a)

asymptotically stable, semi-dissipative matrices
$\mathbf{L}\in{\mathbb{C}}^{n\times n}$, and

  2. (b)

skew-Hermitian matrices $\mathbf{L}\in{\mathbb{C}}^{n\times n}$ (hence
satisfying $\mathbf{L}+\mathbf{L}^{*}=0$).

_Step 4. Reduction of skew-Hermitian matrices to (purely imaginary) scalar
problems._ Skew-Hermitian matrices $\mathbf{L}\in{\mathbb{C}}^{n\times n}$ are
unitarily congruent to a diagonal matrix with purely imaginary eigenvalues,
i.e., there exists a unitary matrix $\mathbf{U}\in{\mathbb{C}}^{n\times n}$
such that $\mathbf{U}\mathbf{L}\mathbf{U}^{*}=\Lambda$ with
$\Lambda=\operatorname{diag}(\lambda_{1}^{\mathbf{L}},\ldots,\lambda_{n}^{\mathbf{L}})$,
where $\lambda_{j}^{\mathbf{L}}\in{\rm i}{\mathbb{R}}$. Therefore, in case
(b), the analysis can be reduced to scalar ODEs of the form

$\frac{\mathrm{d}u}{\mathrm{d}t}=\lambda u,\quad t\geq 0;\quad\lambda\in{\rm
i}{\mathbb{R}}.$

This finishes the proof of Theorem 0.2.1. $\Box$

We already mentioned that the (strong) stability of Runge–Kutta schemes for
scalar ODEs is well studied; see, e.g., KrSc92 . By definition, strong
stability w.r.t. $\mathcal{L}_{IM}$ is equivalent to _local stability on the
imaginary axis_ as defined in (KrSc92, , Definition 2.1) (and recalled in
§0.2). Thus, we deduce from (KrSc92, , Theorem 3.1) the following statement.

###### Theorem 0.3.1

The explicit Runge–Kutta method (4) of order $p\in{\mathbb{N}}$ is strongly
stable w.r.t. $\mathcal{L}_{IM}$ if

$\displaystyle\gamma_{p+1}:=(-1)^{(p+1)/2}(1-c_{p+1})>0$ for $p$ odd,
$\displaystyle\delta_{p+1}:=(-1)^{p/2}\big{(}c_{p+2}-(p+2)c_{p+1}+(p+1)\big{)}>0$
for $p$ even.

If $\gamma_{p+1}<0$ or $\delta_{p+1}<0$, the Runge–Kutta method is not
strongly stable.

Notice that the case $\delta_{p+1}=0$ is left open in (KrSc92, , Theorem 3.1).
Complementary, $\gamma_{p+1}=0$ implies $c_{p+1}=1$, this would contradict
that the Runge–Kutta method has order $p$, see (5).

Examples to illustrate Corollary 1 and Theorem 0.2.2(b). The following two
examples were used in the literature to investigate the strong stability of
the explicit Runge–Kutta method with $s=p=4$. Due to the new result in
Corollary 1, it is clear that this scheme is _not_ strongly stable. However,
in view of Theorem 0.3.1, since $c_{p+1}=c_{p+2}=0$ and hence
$\delta_{p+1}>0$, the method is strongly stable w.r.t. $\mathcal{L}_{IM}$. We
deduce from Theorem 0.2.2(a) that the method is strongly stable w.r.t.
$\mathcal{L}_{AS}^{m}$ with $m=1$. Since (7b) holds, the considered
Runge–Kutta scheme is _not_ strongly stable w.r.t. $\\{\mathbf{L}\\}$ if the
semi-dissipative matrix $\mathbf{L}\in\mathcal{L}_{AS}$ satisfies
$2{m_{HC}}(\mathbf{L})+1>p$, i.e., if ${m_{HC}}(\mathbf{L})>3/2$; see Theorem
0.2.2(b).

###### Example 1

Levy and Tadmor (LeTa98, , §3.5) use the example

$\mathbf{L}=-5\begin{pmatrix}1&2&2&2&2\\\ 0&1&2&2&2\\\ 0&0&1&2&2\\\
0&0&0&1&2\\\ 0&0&0&0&1\end{pmatrix}\quad\text{with
}\mathbf{L}_{H}=-5\begin{pmatrix}1&1&1&1&1\\\ 1&1&1&1&1\\\ 1&1&1&1&1\\\
1&1&1&1&1\\\ 1&1&1&1&1\end{pmatrix}$

to test numerically the strong stability of the explicit Runge–Kutta method
with $s=p=4$. The matrix $\mathbf{L}$ is semi-dissipative, but
$\|R(\tau\mathbf{L})\|>1$ for some $\tau>0$. The staircase form of
$(\mathbf{J},\mathbf{R})=(\mathbf{L}_{S},\mathbf{L}_{H})$ is

$\mathbf{V}\mathbf{L}\mathbf{V}^{*}=-5\left[\begin{array}[]{ccccc}5&-2\sqrt{2}&0&0&0\\\
2\sqrt{2}&0&\frac{\sqrt{35}}{5}&0&0\\\
0&-\frac{\sqrt{35}}{5}&0&-\frac{4\sqrt{35}}{35}&0\\\
0&0&\frac{4\sqrt{35}}{35}&0&\frac{\sqrt{7}}{7}\\\
0&0&0&-\frac{\sqrt{7}}{7}&0\end{array}\right],\quad\mathbf{V}=\left[\begin{array}[]{ccccc}\frac{\sqrt{5}}{5}&\frac{\sqrt{5}}{5}&\frac{\sqrt{5}}{5}&\frac{\sqrt{5}}{5}&\frac{\sqrt{5}}{5}\\\
\frac{\sqrt{10}}{5}&\frac{\sqrt{10}}{10}&0&-\frac{\sqrt{10}}{10}&-\frac{\sqrt{10}}{5}\\\
-\frac{\sqrt{14}}{7}&\frac{\sqrt{14}}{14}&\frac{\sqrt{14}}{7}&\frac{\sqrt{14}}{14}&-\frac{\sqrt{14}}{7}\\\
-\frac{\sqrt{10}}{10}&\frac{\sqrt{10}}{5}&0&-\frac{\sqrt{10}}{5}&\frac{\sqrt{10}}{10}\\\
\frac{\sqrt{70}}{70}&-\frac{2\sqrt{70}}{35}&\frac{3\sqrt{70}}{35}&-\frac{2\sqrt{70}}{35}&\frac{\sqrt{70}}{70}\end{array}\right].$

In particular, $r=6$, $n_{1}=\cdots=n_{5}=1$, $n_{6}=0$. Hence, the matrix
$\mathbf{L}$ has the HC-index ${m_{HC}}=r-2=4$ (see Remark 4).

###### Example 2

To show that the explicit Runge–Kutta method with $s=p=4$ is not strongly
stable, Sun and Shu (SunShu17, , Prop. 1.1) use the counter-example

$\mathbf{L}=-\begin{pmatrix}1&2&2\\\ 0&1&2\\\
0&0&1\end{pmatrix}\quad\text{with }\mathbf{L}_{H}=-\begin{pmatrix}1&1&1\\\
1&1&1\\\ 1&1&1\end{pmatrix}.$

The matrix $\mathbf{L}$ is semi-dissipative, but $\|R(\tau\mathbf{L})\|>1$ for
$\tau>0$ sufficiently small. The staircase form of
$(\mathbf{J},\mathbf{R})=(\mathbf{L}_{S},\mathbf{L}_{H})$ is

$\mathbf{V}\mathbf{L}\mathbf{V}^{*}=-\left[\begin{array}[]{ccc}3&\frac{2\sqrt{3}\,\sqrt{2}}{3}&0\\\
-\frac{2\sqrt{3}\,\sqrt{2}}{3}&0&-\frac{\sqrt{2}\,\sqrt{6}}{6}\\\
0&\frac{\sqrt{2}\,\sqrt{6}}{6}&0\end{array}\right],\quad\mathbf{V}=\left[\begin{array}[]{ccc}\frac{\sqrt{3}}{3}&\frac{\sqrt{3}}{3}&\frac{\sqrt{3}}{3}\\\
-\frac{\sqrt{2}}{2}&0&\frac{\sqrt{2}}{2}\\\
-\frac{\sqrt{6}}{6}&\frac{\sqrt{6}}{3}&-\frac{\sqrt{6}}{6}\end{array}\right].$

Here, we have $r=4$, $n_{1}=n_{2}=n_{3}=1$, $n_{4}=0$, and $\mathbf{L}$ has
the HC-index ${m_{HC}}=r-2=2$. Due to Theorem 0.2.2(b), matrix $\mathbf{L}$
has the minimal HC-index ${m_{HC}}(\mathbf{L})$ to detect the failure of
strong stability, such that ${m_{HC}}(\mathbf{L})>3/2$.

##  0.4 Proof of Theorem 0.2.2

The proof is based on the short-time decay of the spectral norm of the matrix
exponential $t\mapsto e^{t\mathbf{L}}\in{\mathbb{C}}^{n\times n}$,
characterized by the HC-index of $\mathbf{L}$. We recall from (AAC22, ,
Theorem 2.7(a)):

###### Proposition 2

Let the matrix $\mathbf{L}\in{\mathbb{C}}^{n\times n}$ be semi-dissipative.
Then $\mathbf{L}$ is asymptotically stable (with HC-index
${m_{HC}}\in{\mathbb{N}}_{0}$) if and only if

$\|e^{t\mathbf{L}}\|_{2}=1-ct^{a}+\mathcal{O}(t^{a+1})\quad\text{ for
}t\in[0,\epsilon),$ (11)

for some $a,c,\epsilon>0$. In this case, necessarily $a=2m_{HC}+1$.

The sharp multiplicative factor $c$ in (11) has been determined explicitly in
(AAC22, , Theorem 2.7(b)). Definition 3 immediately yields the following
observation.

###### Remark 5

Let $\mathbf{L}$ be semi-dissipative with HC-index
${m_{HC}}\in{\mathbb{N}}_{0}$. In view of (6), there exists a normalized
vector $u_{0}$ such that

$u_{0}\in\ker\big{(}T_{{m_{HC}}-1}\big{)}=\ker\bigg{(}\sum_{j=0}^{{m_{HC}}-1}\mathbf{L}_{S}^{j}\mathbf{L}_{H}(\mathbf{L}_{S}^{*})^{j}\bigg{)},$
(12)

but none that would satisfy instead also
$(-\mathbf{L}_{H})^{1/2}(\mathbf{L}_{S})^{{m_{HC}}}u_{0}=0$. Hence, there
exists a normalized vector $u_{0}$ such that

$\|(-\mathbf{L}_{H})^{1/2}(\mathbf{L}_{S})^{j}u_{0}\|_{2}=0\ \text{ for }\
0\leq
j\leq{m_{HC}}-1,\quad(-\mathbf{L}_{H})^{1/2}(\mathbf{L}_{S})^{{m_{HC}}}u_{0}\neq
0.$ (13)

Preparation. We claim that $t\mapsto\|e^{t\mathbf{L}}\|_{2}$ and
$\tau\mapsto\|R(\tau\mathbf{L})\|_{2}$ are real analytic on a sufficiently
small time interval. In particular, the Taylor expansions of order $p+1$ of
$\|e^{t\mathbf{L}}\|_{2}$ and $\|R(\tau\mathbf{L})\|_{2}$ about zero exist.
Indeed, according to (Ko01, , Lemma 1), for sufficiently small time $t_{0}>0$,
there exists a real analytic function $\Phi:[0,t_{0}]\to{\mathbb{R}}$ such
that $\|e^{t\mathbf{L}}\|_{2}=\Phi(t)$ for all $t\in[0,t_{0}]$. This statement
can also be derived from (Ka95, , Chap. 2, §6) or (HiPr05, , Theorem 4.3.17).

Similarly, we prove that $\|R(\tau\mathbf{L})\|_{2}$ is (at least) locally
real analytic for sufficiently small $\tau\geq 0$. The stability function
$R(\tau\mathbf{L})$ is an analytic matrix function w.r.t. $\tau$, hence
$G(\tau):=R(\tau\mathbf{L}^{*})R(\tau\mathbf{L})$ is again an analytic matrix
function. Due to (HiPr05, , Theorem 4.3.17), there exists a neighborhood of
$\tau_{0}=0$ and locally analytic functions
$\lambda_{j}(\tau)=\sigma_{j}(G(\tau))$, where $\sigma_{j}$ denotes the
singular values of (the self-adjoint) matrix $G(\tau)$, and $\lambda_{j}(0)=1$
for $j=1,...,n$. Hence,

$\|R(\tau\mathbf{L})\|_{2}=\max_{j=1,\ldots,n}\Big{\\{}\sqrt{\lambda_{j}(R(\tau\mathbf{L}))}\Big{\\}}.$

If two real analytic functions are the same on some converging sequence then
they are identical. Since $\lambda_{j}(0)=1$, we use null-sequences to deduce
that any two real analytic functions $\lambda_{k}$, $\lambda_{\ell}$ on
$[0,\infty)$ are either identical or there exists $\tau_{*}^{k\ell}>0$ such
that they do not intersect on a finite interval $[0,\tau_{*}^{k\ell})$. Since
there are only finitely many $\lambda_{j}$, there exists
$\tau_{*}:=\min\\{\tau_{*}^{k\ell}\\}>0$ such that $\|R(\tau\mathbf{L})\|_{2}$
is real analytic on $[0,\tau_{*})$, see also the discussion in (Ka95, , Chap.
2, §6.4).

Proof of Theorem 0.2.2, statement (a). Following the proof of (AAC22, ,
Theorem 2.7), we consider the spectral norm
$\|e^{t\mathbf{L}}\|_{2}^{2}=\|Q(t)\|_{2}=\lambda_{\rm max}(Q(t))$ for small
$t>0$, where

$Q(t):=e^{t\mathbf{L}^{*}}e^{t\mathbf{L}}=\sum_{j=0}^{\infty}\frac{t^{j}}{j!}\
U_{j},\quad
U_{j}=(\mathbf{L}+\mathbf{L}^{*})^{j}=\sum_{k=0}^{j}\binom{j}{k}(\mathbf{L}^{*})^{k}\mathbf{L}^{j-k},$

satisfying $\|U_{j}\|_{2}\leq(2\|\mathbf{L}\|_{2})^{j}$ for
$j\in{\mathbb{N}}_{0}$.

First, we compute the stability function $R(\tau\mathbf{L})$ in (5), setting
$c_{k}=0$ for $k>s$ and $c_{k}=1$ for $k\leq p$,

$\displaystyle G(\tau)$
$\displaystyle=R(\tau\mathbf{L}^{*})R(\tau\mathbf{L})=\sum_{j=0}^{2s}\frac{1}{j!}\sum_{k=0}^{j}\binom{j}{k}c_{k}c_{j-k}(\tau\mathbf{L}^{*})^{k}(\tau\mathbf{L})^{j-k}$
$\displaystyle=\sum_{j=0}^{2s}\frac{\tau^{j}}{j!}\sum_{k=0}^{j}\binom{j}{k}c_{k}c_{j-k}(\mathbf{L}^{*})^{k}\mathbf{L}^{j-k}$
$\displaystyle=\sum_{j=0}^{p}\frac{\tau^{j}}{j!}U_{j}+\sum_{j=p+1}^{2s}\frac{\tau^{j}}{j!}\sum_{k=0}^{j}\binom{j}{k}c_{k}c_{j-k}(\mathbf{L}^{*})^{k}\mathbf{L}^{j-k}$
$\displaystyle=\sum_{j=0}^{p}\frac{\tau^{j}}{j!}U_{j}+\frac{\tau^{p+1}}{(p+1)!}\bigg{(}\sum_{k=1}^{p}\binom{p+1}{k}(\mathbf{L}^{*})^{k}\mathbf{L}^{p+1-k}$
$\displaystyle\phantom{xx}+c_{p+1}\mathbf{L}^{p+1}+c_{p+1}(\mathbf{L}^{*})^{p+1}\bigg{)}+\mathcal{O}(\tau^{p+2}).$

We subtract

$\displaystyle Q(\tau)$
$\displaystyle=\sum_{j=0}^{p}\frac{\tau^{j}}{j!}U_{j}+\frac{\tau^{p+1}}{(p+1)!}\bigg{(}\sum_{k=1}^{p}\binom{p+1}{k}(\mathbf{L}^{*})^{k}\mathbf{L}^{p+1-k}$
$\displaystyle\phantom{xX}+\mathbf{L}^{p+1}+(\mathbf{L}^{*})^{p+1}\bigg{)}+\mathcal{O}(\tau^{p+2})$

from $G(\tau)$ and set $\widetilde{c}_{k}=c_{k}-1$ to find that

$G(\tau)-Q(\tau)=\frac{\tau^{p+1}}{(p+1)!}\widetilde{c}_{p+1}\big{(}\mathbf{L}^{p+1}+(\mathbf{L}^{*})^{p+1}\big{)}+\mathcal{O}(\tau^{p+2}).$
(14)

Consequently, the Taylor expansions for $G(\tau)$ and $Q(\tau)$ differ
starting from the $\tau^{p+1}$-term. Now, consider an explicit Runge–Kutta
scheme of order $p\in{\mathbb{N}}$ with stability function $R(z)$, and
matrices $\mathbf{L}\in\mathcal{L}_{AS}$ with HC-index ${m_{HC}}$ satisfying
$2{m_{HC}}+1\leq p$. Then, taking into account Proposition 2, there exist $c$,
$\epsilon>0$ such that

$\displaystyle\|R(\tau\mathbf{L})\|_{2}$
$\displaystyle=\|G(\tau)\|_{2}^{1/2}=\|Q(\tau)\|_{2}^{1/2}+\mathcal{O}(\tau^{p+1})$
$\displaystyle=\|e^{\tau\mathbf{L}}\|_{2}+\mathcal{O}(\tau^{p+1})=1-c\tau^{2{m_{HC}}+1}+\mathcal{O}(\tau^{2{m_{HC}}+2})$

for all $\tau\in[0,\epsilon)$. We conclude that $\|R(\tau\mathbf{L})\|_{2}\leq
1$ for sufficiently small $\tau$, hence proving that all explicit Runge–Kutta
schemes of order $p\in{\mathbb{N}}$ are strongly stable w.r.t. matrices
$\mathbf{L}\in\mathcal{L}_{AS}$ with HC-index ${m_{HC}}$ satisfying
$2{m_{HC}}+1\leq p$.

Proof of Theorem 0.2.2, statement (b) for even order $p$. Let $p\in
2{\mathbb{N}}$ and let $\mathbf{L}\in\mathcal{L}_{AS}$ with HC-index
${m_{HC}}=p/2$. The proof uses the following result, which follows from
(AAC22, , Lemma A.4), when using the simple observation from (AAC22, , Lemma
2.4) that
$\ker(T_{m})=\ker[\sum_{j=0}^{m}(\mathbf{L}^{*})^{j}\mathbf{L}_{H}\mathbf{L}^{j}]$.

###### Lemma 4

Let $\mathbf{L}\in{\mathbb{C}}^{n\times n}$ be semi-dissipative with HC-index
${m_{HC}}\in{\mathbb{N}}$. Then, for each $u_{0}\in\ker(T_{{m_{HC}}-1})$,
there exists a polynomial vector function $u_{\tau}\in{\mathbb{C}}^{n}$,
$\tau\in[0,1]$, of the form

$u_{\tau}=u_{0}+\sum_{\ell=1}^{{m_{HC}}}b_{\ell}\tau^{\ell}(-\mathbf{L})^{\ell}u_{0}$
(15)

for suitable numbers $b_{\ell}\in{\mathbb{R}}$, such that for $\tau\in[0,1]$,

$g(u_{\tau};\tau):=u_{\tau}^{*}\bigg{(}\sum_{j=1}^{\infty}\frac{\tau^{j}}{j!}U_{j}\bigg{)}u_{\tau}=-2c_{1}(u_{0})\tau^{a}+\mathcal{O}(\tau^{a+1}),$
(16)

where $a=2{m_{HC}}+1$ and
$c_{1}(u_{0}):=\|(-\mathbf{L}_{H})^{1/2}\mathbf{L}^{m_{HC}}u_{0}\|_{2}^{2}/((2{m_{HC}}+1)!\binom{2{m_{HC}}}{{m_{HC}}})$.

The vector function $u_{\tau}$ of this lemma is used as an approximation of an
eigenfunction of $G(\tau)=R(\tau\mathbf{L})^{*}R(\tau\mathbf{L})$, whose
associated eigenvalue is larger than one for sufficiently small $\tau>0$.
Since $g(u_{\tau};\tau)=u_{\tau}^{*}Q(\tau)u_{\tau}-\|u_{\tau}\|_{2}^{2}$, we
can write

$u_{\tau}^{*}G(\tau)u_{\tau}=\|u_{\tau}\|_{2}^{2}+g(u_{\tau};\tau)+u_{\tau}^{*}\big{(}G(\tau)-Q(\tau)\big{)}u_{\tau},$
(17)

where here and in the following, the function $g$ is considered on
${\mathbb{C}}^{n}\times[0,\infty)$. As $\mathbf{L}$ is semi-dissipative,
Remark 5 yields the existence of some
$u_{0}\in\ker(T_{{m_{HC}}-1})\setminus\\{0\\}$ satisfying (13). This shows
that

$\displaystyle
u_{0}^{*}\big{(}\mathbf{L}^{p+1}+(\mathbf{L}^{*})^{p+1}\big{)}u_{0}$
$\displaystyle=u_{0}^{*}\big{(}(\mathbf{L}_{H}+\mathbf{L}_{S})^{p+1}+(\mathbf{L}_{H}-\mathbf{L}_{S})^{p+1}\big{)}u_{0}$
$\displaystyle=2(-1)^{{m_{HC}}}\
u_{0}^{*}(\mathbf{L}_{S}^{*})^{{m_{HC}}}\mathbf{L}_{H}(\mathbf{L}_{S})^{{m_{HC}}}u_{0}$
$\displaystyle=-2(-1)^{{m_{HC}}}\|(-\mathbf{L}_{H})^{1/2}\mathbf{L}_{S}^{{m_{HC}}}u_{0}\|_{2}^{2}\neq
0,$ (18)

where we used from (13) that $\mathbf{L}_{H}(\mathbf{L}_{S})^{j}u_{0}=0$ and
$u_{0}^{*}(\mathbf{L}_{S}^{*})^{j}\mathbf{L}_{H}=0$ for $0\leq
j\leq{m_{HC}}-1$. Therefore, we infer from (14) and (15) that

$\displaystyle u_{\tau}^{*}\big{(}G(\tau)-Q(\tau)\big{)}u_{\tau}$
$\displaystyle=-2(-1)^{m_{HC}}\frac{\tau^{2m_{HC}+1}}{(2m_{HC}+1)!}\widetilde{c}_{p+1}\|(-\mathbf{L}_{H})^{1/2}\mathbf{L}_{S}^{{m_{HC}}}u_{0}\|_{2}^{2}$
$\displaystyle\phantom{xX}+\mathcal{O}(\tau^{2{m_{HC}}+2}).$

Inserting this expression and (16) into (17) yields

$\displaystyle u_{\tau}^{*}$ $\displaystyle
G(\tau)u_{\tau}=\|u_{\tau}\|_{2}^{2}+g(u_{\tau};\tau)+u_{\tau}^{*}\big{(}G(\tau)-Q(\tau)\big{)}u_{\tau}$
$\displaystyle=\|u_{\tau}\|_{2}^{2}-\frac{2\tau^{2{m_{HC}}+1}}{(2{m_{HC}}+1)!}\Big{(}\tbinom{2{m_{HC}}}{{m_{HC}}}^{-1}+(-1)^{{m_{HC}}}\widetilde{c}_{p+1}\Big{)}\|(-\mathbf{L}_{H})^{1/2}\mathbf{L}_{S}^{{m_{HC}}}u_{0}\|_{2}^{2}$
$\displaystyle\phantom{xX}+\mathcal{O}(\tau^{2{m_{HC}}+2}).$ (19)

By condition (7b), the expression in the round brackets is negative,

$\tbinom{2{m_{HC}}}{{m_{HC}}}^{-1}+(-1)^{{m_{HC}}}\widetilde{c}_{p+1}=\tbinom{p}{p/2}^{-1}+(-1)^{p/2}(c_{p+1}-1)<0,$

implying that $u_{\tau}^{*}G(\tau)u_{\tau}>\|u_{\tau}\|_{2}^{2}$ for all
$\tau\in(0,\tau_{*})$ for some $\tau_{*}>0$. We conclude that
$\|R(\tau)\|_{2}>1$ on that interval. Thus, the Runge–Kutta scheme with even
order $p$ is not strongly stable under condition (7b).

Proof of Theorem 0.2.2, statement (b) for odd order $p$. Let
$p\in~{}2{\mathbb{N}}_{0}+1$ and let $\mathbf{L}\in\mathcal{L}_{AS}$ with HC-
index ${m_{HC}}=(p+1)/2$. As before, Remark 5 yields the existence of some
$u_{0}\in\ker(T_{{m_{HC}}-1})\setminus\\{0\\}$ satisfying (13) and hence,
proceeding similarly as in (18),

$\displaystyle
u_{0}^{*}\big{(}\mathbf{L}^{p+1}+(\mathbf{L}^{*})^{p+1}\big{)}u_{0}$
$\displaystyle=u_{0}^{*}\big{(}(\mathbf{L}_{H}+\mathbf{L}_{S})^{2{m_{HC}}}+(\mathbf{L}_{H}-\mathbf{L}_{S})^{2{m_{HC}}}\big{)}u_{0}$
$\displaystyle=2u_{0}^{*}\mathbf{L}_{S}^{2{m_{HC}}}u_{0}=2(-1)^{{m_{HC}}}\|(\mathbf{L}_{S})^{{m_{HC}}}u_{0}\|_{2}^{2}\neq
0.$

Then, with the vector function $u_{\tau}$ as in Lemma 4 and taking into
account (14) and $2m_{HC}=p+1$,

$\displaystyle u_{\tau}^{*}G(\tau)u_{\tau}$
$\displaystyle=\|u_{\tau}\|_{2}^{2}+g(u_{\tau};\tau)+u_{\tau}^{*}\big{(}G(\tau)-Q(\tau)\big{)}u_{\tau}$
$\displaystyle=\|u_{\tau}\|_{2}^{2}-2c_{1}(u_{0})\tau^{2{m_{HC}}+1}+\mathcal{O}(\tau^{2{m_{HC}}+2})$
$\displaystyle\phantom{xX}+\frac{\tau^{p+1}}{(p+1)!}\widetilde{c}_{p+1}u_{0}^{*}\big{(}\mathbf{L}^{p+1}+(\mathbf{L}^{*})^{p+1}\big{)}u_{0}+\mathcal{O}(\tau^{p+2})$
$\displaystyle=\|u_{\tau}\|_{2}^{2}+\frac{\tau^{p+1}}{(p+1)!}\widetilde{c}_{p+1}u_{0}^{*}\big{(}\mathbf{L}^{p+1}+(\mathbf{L}^{*})^{p+1}\big{)}u_{0}+\mathcal{O}(\tau^{p+2})$
$\displaystyle=\|u_{\tau}\|_{2}^{2}+\frac{2\tau^{p+1}}{(p+1)!}(-1)^{(p+1)/2}\widetilde{c}_{p+1}\|\mathbf{L}_{S}^{(p+1)/2}u_{0}\|_{2}^{2}+\mathcal{O}(\tau^{p+2}).$
(20)

Since $\|\mathbf{L}_{S}^{(p+1)/2}u_{0}\|_{2}^{2}\neq 0$, condition (7a)
implies that

$(-1)^{(p+1)/2}\widetilde{c}_{p+1}=(-1)^{(p+1)/2}(c_{p+1}-1)>0$

and consequently, $u_{\tau}^{*}G(\tau)u_{\tau}>\|u_{\tau}\|_{2}^{2}$ for all
$\tau\in(0,\tau_{*})$ for some $\tau_{*}>0$. Thus, the Runge–Kutta scheme with
odd order $p$ is not strongly stable under condition (7a).

Proof of Theorem 0.2.2, statement (c). Consider an explicit Runge–Kutta scheme
of order $p\in{\mathbb{N}}$ which is not strongly stable and let (7) hold. If
a matrix $\mathbf{L}\in\mathcal{L}_{AS}$ with HC-index
${m_{HC}}={m_{HC}}(\mathbf{L})$ satisfies $2{m_{HC}}+1>p$, then we can repeat
the construction in the proof of statement (b): Consider $m\in{\mathbb{N}}$
such that $2m+1>p$ but $2m-1\leq p$. Due to ${m_{HC}}\geq m$, there exists a
vector $u_{0}\in{\mathbb{C}}^{n}$ such that

$\|(-\mathbf{L}_{H})^{1/2}\mathbf{L}_{S}^{j}u_{0}\|_{2}=0\text{ for }0\leq
j\leq m-1,\quad\|(-\mathbf{L}_{H})^{1/2}\mathbf{L}_{S}^{m}u_{0}\|_{2}\neq 0,$

see (13), but with $m$ replacing ${m_{HC}}$. Using a vector function
$u_{\tau}$ for $\tau\in[0,1]$ as constructed in (15), we can use the above
equations (19), (20), again with $m$ replacing ${m_{HC}}$, to prove statement
(c) in Theorem 0.2.2. $\Box$

##  0.5 Proof of Theorem 0.2.4

As pointed out in (LeTa98, , §3.5), an explicit Runge–Kutta method of order
$p=4$ with $s=4$ stages can be made strongly stable w.r.t. $\mathcal{L}_{AS}$
when restricting to some adapted $\mathbf{P}$-norm. Motivated by this
observation, we prove Theorem 0.2.4.

First, we observe that for strong stability in weak form, the analogous
statement as in Theorem 0.2.1 holds, i.e., an explicit Runge–Kutta method is
strongly stable in weak form if and only if it is strongly stable in weak form
w.r.t. $\mathcal{L}_{AS}$ and $\mathcal{L}_{IM}$.

Second, let an asymptotically stable, semi-dissipative matrix
$\mathbf{L}\in\mathcal{L}_{AS}$ be given. A matrix
$\mathbf{L}\in{\mathbb{C}}^{n\times n}$ is asymptotically stable if and only
if there exists a matrix $\mathbf{P}\in\mathbb{H}_{n}^{>}$ such that

$\mathbf{L}^{*}\mathbf{P}+\mathbf{P}\mathbf{L}<0,$ (21)

see e.g. (Be18, , Corollary 15.10.1), (HadChe08, , Theorem 3.17). Following
the proof of Lemma 1, we find that
$\widehat{\mathbf{L}}:=\mathbf{P}^{1/2}\mathbf{L}\mathbf{P}^{-1/2}$ is
dissipative. Hence, $\widehat{\mathbf{L}}$ has the HC-index
${m_{HC}}(\widehat{\mathbf{L}})=0$. Consider the solution
$(u^{k})_{k\in{\mathbb{N}}_{0}}$ in (4) in the weighted norm
$\|\cdot\|_{\mathbf{P}}$ with $\mathbf{P}\in\mathbb{H}_{n}^{>}$ satisfying
(21). By Step 1 of the proof of Theorem 0.2.1, strong stability in weak form
for the Runge–Kutta scheme, i.e.
$\|u^{1}\|_{\mathbf{P}}\leq\|u^{0}\|_{\mathbf{P}}$, is equivalent to

$\|\widehat{u}^{1}\|\leq\|\widehat{u}^{0}\|,\quad\mbox{where
}\widehat{u}^{1}=R(\tau\widehat{\mathbf{L}})\widehat{u}^{0}.$ (22)

The inequality in (22) indeed holds due to Theorem 0.2.2(a), since
$\widehat{\mathbf{L}}$ has the HC-index ${m_{HC}}=0$. Thus,
$\|u^{1}\|_{\mathbf{P}}\leq\|u^{0}\|_{\mathbf{P}}$ holds without (further)
constraints on $\mathbf{L}$.

Finally, consider strong stability w.r.t. $\mathcal{L}_{IM}$, the set of
purely imaginary scalars $\lambda\in{\rm i}{\mathbb{R}}$. Introducing a
“weighted norm” $\|\cdot\|_{\mathbf{P}}$ for some positive constant
$\mathbf{P}\in{\mathbb{R}}^{+}$ does not change the stability criterion
$\|u^{1}\|_{\mathbf{P}}\leq\|u^{0}\|_{\mathbf{P}}$ on the considered explicit
Runge–Kutta scheme. This finishes the proof. $\Box$

######  Example 3 (Continuation of Example 2)

While the explicit Runge-Kutta scheme with $s=p=4$ steps is _not_ strongly
stable (due to Corollary 1), it is strongly stable in weak form (due to
Theorem 0.2.4, since it is strongly stable w.r.t. $\mathcal{L}_{IM}$; see
Table 1). Notice that the authors of SunShu17  study this problem in the
Euclidean norm, whereas we use an adapted $\mathbf{P}$-norm to guarantee the
decay estimate $\|u^{1}\|_{\mathbf{P}}\leq\|u^{0}\|_{\mathbf{P}}$.

In fact, the Jordan normal form
$-\mathbf{L}=(\mathbf{W}^{*})^{-1}\mathbf{J}\mathbf{W}^{*}$ of $-\mathbf{L}$
holds with

$\mathbf{J}=\begin{pmatrix}1&1&0\\\ 0&1&1\\\
0&0&1\end{pmatrix}\quad\text{and}\quad\mathbf{W}=\begin{pmatrix}1&0&0\\\
0&2&0\\\ 0&2&4\end{pmatrix}.$

Moreover, a suitable matrix $\mathbf{P}=\mathbf{P}(\mathbf{L})$ is given as
$\mathbf{P}=\mathbf{W}\mathbf{W}^{*}>0$. Then, the explicit Runge-Kutta scheme
with $s=p=4$ stages is strongly stable w.r.t. $\\{\mathbf{L}\\}$ and this
specific $\mathbf{P}$-norm.

Acknowledgement: The authors were supported by the Austrian Science Fund (FWF)
via the FWF-funded SFB # F65. The third author acknowledges support from the
FWF, # P33010. This work received funding from the European Research Council
(ERC) under the European Union’s Horizon 2020 research and innovation
programme, ERC Advanced Grant NEUROMORPH, no. 101018153.

## References

  * (1) F. Achleitner, A. Arnold, and E. Carlen.  On multi-dimensional hypocoercive BGK models. Kinet. Relat. Models 11 (2018), 953–1009.
  * (2) F. Achleitner, A. Arnold, and E. A. Carlen.  The hypocoercivity index for the short time behavior of linear time-invariant ODE systems. J. Differential Equations 371 (2023), 83–115.
  * (3) F. Achleitner, A. Arnold, and V. Mehrmann.  Hypocoercivity and controllability in linear semi-dissipative Hamiltonian ordinary differential equations and differential-algebraic equations.  Z. Angew. Math. Mech. 103 (2023), e202100171, 30 pages.
  * (4) F. Achleitner, A. Arnold, and V. Mehrmann.  Hypocoercivity and hypocontractivity concepts for linear dynamical systems.  Electronic J. Linear Algebra 39 (2023), 33–61. DOI: https://doi.org/10.13001/ela.2023.7531
  * (5) D. S. Bernstein,  Scalar, vector, and matrix mathematics. Theory, facts, and formulas. Revised and expanded edition. Princeton University Press, Princeton, NJ, 2018.
  * (6) S. Conda, S. Gottlieb, Z. J. Grant, and J. N. Shadid.  Implicit and implicit-explicit strong stability preserving Runge–Kutta methods with high linear order. J. Sci. Comput. 73 (2017), 667–690.
  * (7) L. Ferracina and M. N. Spijker.  Stepsize restrictions for the total-variation-diminishing property in general Runge–Kutta methods. SIAM J. Numer. Anal. 42 (2004), 1073–1093.
  * (8) S. Gottlieb, C.-W. Shu, and E. Tadmor.  Strong-stability-preserving high-order time discretization methods. SIAM Rev. 43 (2001), 89–112.
  * (9) W. M. Haddad and V. Chellaboina.  Nonlinear dynamical systems and control. A Lyapunov-based approach. Princeton University Press, Princeton, NJ, 2008.
  * (10) E. Hairer and G. Wanner.  Solving Ordinary Differential Equations. II. Stiff and Differential-algebraic Problems.  Second edition. Springer, Berlin, 1996.
  * (11) D. Hinrichsen and A. J. Pritchard.  Mathematical systems theory. I, volume 48 of Texts in Applied Mathematics.  Springer, Berlin, 2005.
  * (12) I. Higueras.  On strong stability preserving time discretization methods. J. Sci. Comput. 21 (2004), 193–223.
  * (13) T. Kato.  Perturbation theory for linear operators.  Classics in Mathematics. Springer-Verlag, Berlin, 1995.  Reprint of the 1980 edition.
  * (14) L. Kohaupt.  Differential calculus for some $p$-norms of the fundamental matrix with applications.  J. Comput. Appl. Math. 135 (2001), 1–21.
  * (15) H.-O. Kreiss and G. Scherer.  Method of lines for hyperbolic differential equations. SIAM J. Numer. Anal. 29 (1992), 640–646.
  * (16) D. Levy and E. Tadmor.  From semidiscrete to fully discrete: Stability of Runge–Kutta schemes by the energy method. SIAM Rev. 40 (1998), 40–73.
  * (17) C. Mehl, V. Mehrmann, and P. Sharma, Stability radii for linear Hamiltonian systems with dissipation under structure-preserving perturbations. SIAM J. Matrix Anal. Appl. 37 (2016), 1625–1654.
  * (18) C. Mehl, V. Mehrmann, and M. Wojtylak. Linear algebra properties of dissipative Hamiltonian descriptor systems, SIAM J. Matrix Anal. Appl. 39 (2018), 1489–1519.
  * (19) C. Mehl, V. Mehrmann, and M. Wojtylak, Distance problems for dissipative hamiltonian systems and related matrix polynomials, Linear Algebra Appl., 623 (2021), 335–366.
  * (20) H. Ranocha.  On strong stability of explicit Runge–Kutta methods for nonlinear semibounded operators. IMA J. Numer. Anal. 41 (2021), 654–682.
  * (21) H. Ranocha and D. I. Ketcheson.  Energy stability of explicit Runge–Kutta methods for nonautonomous or nonlinear problems. SIAM J. Numer. Anal. 58 (2020), 3382–3405.
  * (22) C.-W. Shu.  Total-variation-diminishing time discretizations. SIAM J. Sci. Statist. Comput. 9 (1988), 1073–1084.
  * (23) C.-W. Shu and S. Osher.  Efficient implementation of essentially nonoscillatory shock-capturing schemes. J. Comput. Phys. 77 (1988), 439–471.
  * (24) T. Ström.  On logarithmic norms. SIAM J. Numer. Anal. 12 (1975), 741–753.
  * (25) Z. Sun and C.-W. Shu.  Stability of the fourth order Runge–Kutta method for time-dependent partial differential equations. Ann. Math. Sci. Appl. 2 (2017), 255–284.
  * (26) Z. Sun and C.-W. Shu.  Strong stability of explicit Runge–Kutta time discretizations. SIAM J. Numer. Anal. 57 (2019), 1158–1182.
  * (27) Z. Sun, Y. Wei and K. Wu.  On energy laws and stability of Runge-Kutta methods for linear seminegative problems. SIAM J. Numer. Anal. 60 (2022), 2448–2481.
  * (28) E. Tadmor.  From semidiscrete to fully discrete: Stability of Runge–Kutta schemes by the energy method. II. In: D. Estep and S. Tavener (eds.), Collected Lectures on the Preservation of Stability under Discretization. Proc. Appl. Math. 109, pp. 25–49. SIAM, Philadelphia, 2002.


# Shocks and solitary waves in series connected discrete Josephson
transmission lines

Eugene Kogan<EMAIL_ADDRESS>Jack and Pearl Resnick Institute,
Department of Physics, Bar-Ilan University, Ramat-Gan 52900, Israel

###### Abstract

We analytically study shocks and solitary waves in the discrete Josephson
transmission line (JTL), constructed from Josephson junctions (JJs) and ground
capacitors. Our approach is based on the quasi-continuum approximation, which
we discuss in details. We find both the profiles of the waves, and their
propagation velocities. We also study the effects of the shunting capacitors
and of the resistors on the waves profiles. We find that the resistors
shunting the JJs and/or in series with the ground capacitors lead to the
existence of shock waves of more general type than those existing in the
lossless JTL. We suggest the generalized quasi-continuum approximation.

##  I Introduction

The concept that in a nonlinear wave propagation system the various parts of
the wave travel with different velocities, and that wave fronts (or tails) can
sharpen into shock waves, is deeply imbedded in the classical theory of fluid
dynamics whitham . The methods developed in that field can be profitably used
to study signal propagation in nonlinear transmission lines french ; nouri ;
neto ; nikoo ; silva ; wang ; rangel ; kyuregyan ; akem ; fairbanks . In the
early studies of shock waves in transmission lines, the origin of the
nonlinearity was due to nonlinear capacitance in the circuit landauer ; peng ;
rabinovich .

Interesting and potentially important examples of nonlinear transmission lines
are circuits containing Josephson junctions (JJs) josephson  \- Josephson
transmission lines (JTLs) barone ; pedersen ; tinkham ; kadin . The unique
nonlinear properties of JTLs allow to construct soliton propagators, microwave
oscillators, mixers, detectors, parametric amplifiers, and analog amplifiers
pedersen ; kadin ; tinkham .

Transmission lines formed by JJs connected in series were studied beginning
from 1990s, though much less than transmission lines formed by JJs connected
in parallel solitons . However, the former began to attract quite a lot of
attention recently yaakobi ; brien ; macklin ; kochetov ; zorin ; basko ;
dixon ; goldstein , especially in connection with possible JTL traveling wave
parametric amplification white ; miano ; pekker .

In our previous publication kogan  we studied the propagation of shock waves
in the discrete JTLs with resistive elements. The emphasis of that paper was
analysis of the influence of the resistive elements on the width (and the
shape) of the shocks. Now we want to analyse the structure of the shocks and
solitary waves in the dissipationless JTL. We improve the continuum
approximation used previously, thus taking into account the intrinsic
dispersion of the discrete JTL.

The rest of the paper is constructed as follows. In Section II we formulate
quasi-continuum approximation for the equations describing discrete lossless
JTL. These equations are solved in Section III for the case of shock and
solitary waves. In Section IV we rigorously justify the quasi-continuum
approximation for the shock and solitary waves in certain limiting cases. In
Section V we discuss the effect of the shunting capacitors and resistive
elements on the shock and solitary waves propagation in the discrete JTL.
Application of the results obtained in the paper and opportunities for their
generalization are very briefly discussed in Section VI. We conclude in
Section VII. In the Appendix we illustrate the advantage of the quasi-
continuum approximation over the continuum one. The Appendix B formulates the
generalized quasi-continuum approximation, by introducing the appropriate
integral operator.

##  II Quasi-continuum approximation for the Discrete Josephson transmission
line

Consider the dissipationless model of JTL constructed from identical JJs and
capacitors, which is shown on Fig. 1.

Figure 1: Discrete JTL.

We take as dynamical variables the phase differences $\varphi_{n}$ across the
JJs and the charges $q_{n}$ which have passed through the JJs. The circuit
equations are

$\displaystyle\frac{\hbar}{2e}\frac{d\varphi_{n}}{dt}$
$\displaystyle=\frac{1}{C}\left(q_{n+1}-2q_{n}+q_{n-1}\right)\;,$ (1a)
$\displaystyle\frac{dq_{n}}{dt}$ $\displaystyle=I_{c}\sin\varphi_{n}\;,$ (1b)

where $C$ is the ground capacitor, and $I_{c}$ is the critical current of the
JJ. Differentiating Eq. (1a) with respect to $t$ and substituting $dq_{n}/dt$
from Eq. (1b), we obtain closed equation for $\varphi_{n}$ mohebbi ; chen ;
katayama ; kogan

$\displaystyle
L_{J}C\frac{d^{2}\varphi_{n}}{dt^{2}}=\sin\varphi_{n+1}-2\sin\varphi_{n}+\sin\varphi_{n-1}\,,$
(2)

where $L_{J}=\hbar/(2eI_{c})$.

Treating $\varphi$ as a function of the continuous variable $z$, we can write
down the r.h.s. of Eq. (2) as

$\displaystyle\text{r.h.s.}=\Lambda^{2}\frac{\partial^{2}\sin\varphi}{\partial
z^{2}}+\frac{\Lambda^{4}}{12}\frac{\partial^{4}\sin\varphi}{\partial
z^{4}}+\frac{\Lambda^{6}}{360}\frac{\partial^{6}\sin\varphi}{\partial
z^{6}}+\dots\,,$ (3)

where $\Lambda$ is the period of the line. In this paper we will truncate the
expansion, taking into account only the first two terms of the expansion, thus
obtaining Eq. (3) in the form

$\displaystyle\frac{L_{J}C}{\Lambda^{2}}\frac{\partial^{2}\varphi}{\partial
t^{2}}=\frac{\partial^{2}\sin\varphi}{\partial
z^{2}}+\frac{\Lambda^{2}}{12}\frac{\partial^{4}\sin\varphi}{\partial
z^{4}}\,.$ (4)

We will call such truncation the quasi-continuum approximation. To illustrate
the approximation, we apply it in the Appendix to the signalling problem for
the transmission line with linear inductances instead of the JJs, and analyse
how well it reproduces the exact analytical solution, which can be found in
that case.

Looking for the running wave solution of Eq. (4) in the form
$\varphi(z,t)=\varphi(x)$, where $x=Vt-z$, we obtain an ordinary differential
equation

$\displaystyle\frac{\Lambda^{2}}{12}\frac{d^{4}\sin\varphi}{dx^{4}}+\frac{d^{2}\sin\varphi}{dx^{2}}-\overline{V}^{2}\frac{d^{2}\varphi}{dx^{2}}=0\,,$
(5)

where $\overline{V}^{2}\equiv V^{2}L_{J}C/\Lambda^{2}$.

##  III Shocks and solitary waves

###  III.1 Newtonian equation

We are interested in the propagation of waves characterised by the boundary
conditions

$\displaystyle\lim_{x\to-\infty}\varphi=\varphi_{2}\;,\hskip
28.45274pt\lim_{x\to+\infty}\varphi=\varphi_{1}$ (6)

Integrating Eq. (5) with respect to $x$ twice, we obtain

$\displaystyle\frac{\Lambda^{2}}{12}\frac{d^{2}\sin\varphi}{dx^{2}}=-\sin\varphi+\overline{V}^{2}\varphi+F\,,$
(7)

where $F$ is the constant of integration (the other constant of integration is
equal to zero because of the boundary conditions). We can think about $x$ as
time and about $\sin\varphi$ as coordinate of the fictitious particle, thus
considering (7) as the Newtonian equation. The problem of finding the profile
of the wave is reduced to studying the the motion of the particle which is at
rest both initially at $x\to-\infty$ and finally at $x\to+\infty$.

Going in Eq. (5) to the limits $x\to+\infty$ and $x\to-\infty$, we obtain

$\displaystyle\overline{V}^{2}\varphi_{1}$
$\displaystyle=\sin\varphi_{1}-F\;,$ (8a)
$\displaystyle\overline{V}^{2}\varphi_{2}$
$\displaystyle=\sin\varphi_{2}-F\;.$ (8b)

Solving (8) relative to $\overline{V}^{2}$ and $F$ we obtain kogan

$\displaystyle\overline{V}^{2}$
$\displaystyle=\overline{V}^{2}_{sh}(\varphi_{1},\varphi_{2})\equiv\frac{\sin\varphi_{2}-\sin\varphi_{1}}{\varphi_{2}-\varphi_{1}}\;,$
(9a) $\displaystyle F$
$\displaystyle=\frac{\varphi_{2}\sin\varphi_{1}-\varphi_{1}\sin\varphi_{2}}{\varphi_{2}-\varphi_{1}}\;.$
(9b)

Multiplying Eq. (7) by $d\sin\varphi/dx$ and integrating once again we obtain

$\displaystyle\frac{\Lambda^{2}}{24}\left(\frac{d\sin\varphi}{dx}\right)^{2}+\Pi(\sin\varphi)=0\,,$
(10)

where

$\displaystyle\Pi(\sin\varphi)=\frac{1}{2}\sin^{2}\varphi-\overline{V}^{2}(\varphi\sin\varphi+\cos\varphi)-F\sin\varphi-E\,,$

where $E$ is another constant of integration.

The motion in question either connects two maxima of equal height of the
”potential energy” $\Pi(\sin\varphi)$ at $\sin\varphi_{1}$ and at
$\sin\varphi_{2}$, or starts from a maximum at $\sin\varphi_{1}$ and returns
there, after reflection from a potential wall. In both cases
$E=\Pi(\sin\varphi_{1})$.

In the first case the solution describes the shock wave, in the second - the
solitary wave. Introducing the term ”shock wave” we actually jumped ahead and
took into account that the length scale at which $\varphi$ changes
substantially for the solution in question will be shown to be of order of
$\Lambda$, which is a ”microscopical” length in our case. (The same will be
true for the solitary wave.)

###  III.2 Shock wave

In this case from equal height of the maxima follows

$\displaystyle\varphi_{2}$ $\displaystyle=-\varphi_{1}\;,$ (12a)
$\displaystyle\overline{V}^{2}$
$\displaystyle=\overline{V}^{2}_{sh}(\varphi_{1},-\varphi_{1})=\frac{\sin\varphi_{1}}{\varphi_{1}}\;,$
(12b) $\displaystyle F$ $\displaystyle=0\;,$ (12c)

and for the ”potential energy” we obtain

$\displaystyle\Pi(\sin\varphi)$ $\displaystyle=$
$\displaystyle\frac{1}{2}(\sin\varphi-\sin\varphi_{1})^{2}$ (13)
$\displaystyle-$
$\displaystyle\frac{\sin\varphi_{1}}{\varphi_{1}}[(\varphi-\varphi_{1})\sin\varphi+\cos\varphi-\cos\varphi_{1}]\,.$

Equation (13) and the results of integration of Eq. (10) for this ”potential
energy” are graphically presented on Figs. 2, and 3.

Figure 2: The ”potential energy” (13). We have chosen $\varphi_{1}=.5$. Figure
3: The shock wave profile calculated according to Eq. (10). We have chosen
$\varphi_{1}=.5$.

Consider specifically the limiting case $|\varphi_{1}|\ll 1$. Expanding the
”potential energy” with respect to $\varphi$ and $\varphi_{1}$ and keeping
only the lowest order terms we obtain the approximation to Eq. (10) in the
form

$\displaystyle\Lambda^{2}\left(\frac{d\varphi}{dx}\right)^{2}=\left(\varphi_{1}^{2}-\varphi^{2}\right)^{2}\,.$
(14)

The solution of Eq. (14) is

$\displaystyle\varphi(x)=\pm\varphi_{1}\tanh\frac{\varphi_{1}x}{\Lambda}\,,$
(15)

where plus/minus signs correspond to two solutions: plus - with negative
$\varphi_{1}$, and minus - with positive $\varphi_{1}$. Equations (15) and Eq.
(12b), expanded in series with respect to $\varphi_{1}$ and truncated after
the first two terms
$\overline{V}^{2}_{sh}(\varphi_{1}-\varphi_{1})=1-\frac{\varphi_{1}^{2}}{6}$,
coincide with those obtained by Katayama et al. katayama .

###  III.3 Solitary wave

For a solitary wave two equations of (8) become one equation. As an additional
parameter we take the maximally different from $\varphi_{1}$ value of
$\varphi$, which we will designate as $\varphi_{0}$. This brings another
equation, and $\overline{V}^{2}$ and $F$ are found from the system

$\displaystyle\overline{V}^{2}\varphi_{1}$
$\displaystyle=\sin\varphi_{1}-F\;,$ (16a) $\displaystyle\Pi(\sin\varphi_{0})$
$\displaystyle=0\;.$ (16b)

Solving (16) we obtain the velocity

$\displaystyle\overline{V}^{2}_{sol}(\varphi_{0},\varphi_{1})=\frac{\left(\sin\varphi_{0}-\sin\varphi_{1}\right)^{2}}{2[(\varphi_{0}-\varphi_{1})\sin\varphi_{0}+\cos\varphi_{0}-\cos\varphi_{1}]}\,,$

and the ”potential energy”

$\displaystyle\Pi(\sin\varphi)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\sin\varphi-\sin\varphi_{1}\right)^{2}-\overline{V}^{2}_{sol}(\varphi_{0},\varphi_{1})$
(18) $\displaystyle\cdot$
$\displaystyle\left[(\varphi-\varphi_{1})\sin\varphi+\cos\varphi-\cos\varphi_{1}\right]\,.$

Equation (18) is graphically presented on Fig. 4,

Figure 4: The ”potential energy” (18). We have chosen $\varphi_{1}=-1.0$, and
$\varphi_{0}=.5$.

Considering the limiting case $|\varphi_{1}|,|\varphi_{0}|\ll 1$, expanding
Eq. (18) with respect to all the phases and keeping only the lowest order
terms we obtain Eq. (10) in the form

$\displaystyle\Lambda^{2}\left(\frac{d\varphi}{dx}\right)^{2}=\left(\varphi-\varphi_{1}\right)^{2}\left(\varphi-\varphi_{0}\right)\left(\varphi+2\varphi_{1}+\varphi_{0}\right)\,.$
(19)

Equation (19) can be integrated in elementary functions and we obtain

$\displaystyle\sqrt{\frac{(3\varphi_{1}+\varphi_{0})(\varphi-\varphi_{0})}{(\varphi+2\varphi_{1}+\varphi_{0})(\varphi_{1}-\varphi_{0})}}$
$\displaystyle=\tanh\frac{\sqrt{(3\varphi_{1}+\varphi_{0})(\varphi_{1}-\varphi_{0})}|x|}{2\Lambda}\,.$
(20)

Equation (III.3) is graphically presented on Fig. 5.

Figure 5: Solitary wave profile according to Eq. (III.3). We have chosen
$\varphi_{1}$ to be negative and $\varphi_{0}=-.5\varphi_{1}$.

In an another limiting case $|\varphi_{1}-\varphi_{0}|\ll|\varphi_{1}|$, Eq.
(10) takes the form

$\displaystyle\Lambda^{2}\left(\frac{d\varphi}{dx}\right)^{2}=4\tan\varphi_{1}\cdot\left(\varphi-\varphi_{1}\right)^{2}\left(\varphi-\varphi_{0}\right)\,.$
(21)

The solution of Eq. (21) is

$\displaystyle\sqrt{\frac{\varphi-\varphi_{0}}{\varphi_{1}-\varphi_{0}}}=\tanh\frac{\sqrt{\tan\varphi_{1}\cdot(\varphi_{1}-\varphi_{0})}|x|}{\Lambda}\,.$
(22)

Note that Eq. (22), for $|\varphi_{1}|\ll 1$, coincides with Eq. (III.3), for
$\varphi\approx\varphi_{0}\approx\varphi_{1}$.

##  IV Controlled quasi-continuum approximation

Let us return to Eq. (3). Looking at Eqs. (15) and (III.3) we realize with
hindsight, that in the description of shock and solitary waves with
$|\varphi_{1}|\ll 1$, the expansion parameter in the r.h.s. of Eq. (3) is
$\varphi_{1}^{2}$; thus the quasi-continuum approximation (4) can be
rigorously justified. However, strictly speaking, the truncation of the
expansion should be performed in accordance with the truncation of the series
expansion of the sine function, and Eq. (4), in the consistent approximation
should be written as

$\displaystyle\frac{L_{J}C}{\Lambda^{2}}\frac{\partial^{2}\varphi}{\partial
t^{2}}=\frac{\partial^{2}\varphi}{\partial
z^{2}}-\frac{1}{6}\frac{\partial^{2}\varphi^{3}}{\partial
z^{2}}+\frac{\Lambda^{2}}{12}\frac{\partial^{4}\varphi}{\partial z^{4}}\,.$
(23)

Note that shock wave solution of Eq. (23) for the running wave with the
boundary conditions (6) gives (15).

Looking at Eq. (22) we realize that alternatively, the quasi-continuum
approximation can be rigorously justified when it is applied to the
description of the solitary wave with $|\varphi_{1}-\varphi_{0}|\ll 1$. The
expansion parameter in the r.h.s. of Eq. (3) in this case is
$\sin\varphi_{1}\cdot(\varphi_{1}-\varphi_{0})$. Again, strictly speaking,
truncation of the expansion should be performed in accordance with the
truncation of the series expansion of the sine function, and Eq. (4) in this
case in the consistent approximation should be written as

$\displaystyle\frac{L_{J}C}{\Lambda^{2}}\frac{\partial^{2}\psi}{\partial
t^{2}}=\cos\varphi_{1}\frac{\partial^{2}\psi}{\partial
z^{2}}-\frac{1}{2}\sin\phi_{1}\frac{\partial^{2}\psi^{2}}{\partial z^{2}}$
$\displaystyle+\cos\varphi_{1}\frac{\Lambda^{2}}{12}\frac{\partial^{4}\psi}{\partial
z^{4}}\,,$ (24)

where $\psi=\varphi-\varphi_{1}$. Note that solitary solution of Eq. (IV) for
the running wave with the boundary conditions (6) and the amplitude
$\varphi_{0}$ gives (22).

Here we would like to attract the attention of the reader to the following
issue. Common wisdom says that the continuum approximation and the small
amplitude approximation are independent - there could be a wave with small
amplitude, which allows to expand the sine function, but which varies fast in
space (wavelength comparable to lattice spacing), so the continuum limit is
not justified. And there could be the opposite situation (large amplitude,
long wavelength), in which the sine needs to be retained but the continuum
limit is allowed.

However, for the shocks and solitary waves these approximations are not
independent. Parametrically, the length scale of the waves is of the order of
the lattice spacing $\Lambda$, so, naively, the continuum (or even the quasi-
continuum) limit is not justified. What we have shown above, is that for the
waves with small amplitude $|\varphi_{1}|$
($\tan\varphi_{1}(\varphi_{1}-\varphi_{0})$), the length scale is
$\Lambda/|\varphi_{1}|$
($\Lambda/(\tan\varphi_{1}(\varphi_{1}-\varphi_{0}))$), thus justifying the
quasi-continuum approximation.

We were able to rigorously justify the results of Section III in the two
limiting cases. The answer to the question, whether Eq. (4), taken at its face
value, gives qualitatively correct description of the shock and solitary wave
with non-small amplitude, we postpone until later time.

##  V The effect of shunting capacitors and resistive elements

Consider JTL with the capacitor and resistance shunting the JJ and another
resistance in series with the ground capacitor, shown on Fig. 6.

Figure 6: JTL with the capacity and resistance shunting the JJ and another
resistance in series with the ground capacitor

In this case we chose as the dynamical variables the charges $q_{n}$ which
have passed through all the three branches. As the result, Eq. (1) changes to

$\displaystyle\frac{\hbar}{2e}\frac{d\varphi_{n}}{dt}$
$\displaystyle=\left(\frac{1}{C}+R\frac{\partial}{\partial
t}\right)\left(q_{n+1}-2q_{n}+q_{n-1}\right)\;,$ (25a)
$\displaystyle\frac{dq_{n}}{dt}$
$\displaystyle=I_{c}\sin\varphi_{n}+\frac{\hbar}{2eR_{J}}\frac{d\varphi_{n}}{dt}+C_{J}\frac{\hbar}{2e}\frac{d^{2}\varphi_{n}}{dt^{2}}\;,$
(25b)

where $R$ is the ohmic resistor in series with the ground capacitor, and
$C_{J}$ and $R_{J}$ are the capacitor and the ohmic resistor shunting the JJ.
Equation (4) changes to

$\displaystyle\frac{L_{J}C}{\Lambda^{2}}\frac{\partial^{2}\varphi}{\partial
t^{2}}=\left(1+RC\frac{\partial}{\partial
t}\right)\left(\frac{\partial^{2}}{\partial
z^{2}}+\frac{\Lambda^{2}}{12}\frac{\partial^{4}}{\partial z^{4}}\right)$
$\displaystyle\left(\sin\varphi+\frac{L_{J}}{R_{J}}\frac{\partial\varphi}{\partial
t}+L_{J}C_{J}\frac{\partial^{2}\varphi}{\partial t^{2}}\right)\,.$ (26)

Opening the parenthesis in Eq. (V) and considering the running wave solutions
we obtain instead of Eq. (7) a ”Newtonian” equation

$\displaystyle\left(\frac{C_{J}}{C}+\frac{R}{R_{J}}\right)\overline{V}^{2}\Lambda^{2}\frac{d^{2}\varphi}{dx^{2}}+\frac{\Lambda^{2}}{12}\frac{d^{2}\sin\varphi}{dx^{2}}$
$\displaystyle+V\left(RC\cos\varphi+\frac{L_{J}}{R_{J}}\right)\frac{d\varphi}{dx}=-\sin\varphi+\overline{V}^{2}\varphi+F$
(27)

(we ignored the terms with the derivatives higher than of the forth order).

Note that Eq. (V) can be considered as the balance between the dispersion
effects, described by the l.h.s. of the equation, and the nonlinear effects
described by the (nonlinear terms of the) r.h.s. of the equation. Note also
that when we go to the continuum limit ($\Lambda\to 0$), $C,L_{J}$ and $R_{J}$
scale as $\Lambda$, and $C_{J}$ and $R$ scale as inverse $\Lambda$, so the
first and the third term in the l.h.s. of Eq. (V) don’t disappear, contrary to
what happens with the second terms. In this limit, equation (V) again becomes
Newtonian, with $\varphi$ as the coordinate kogan .

Now let us compare (V) with (7). We realize that Eq. (8) and following from
the former Eq. (9) is still valid for (V). We realize that the shunting
capacitors (at least qualitatively) just renormalize the ”mass” of the
particle. On the other hand, the resistive elements by introducing the
effective ”friction force” break the (analog of the) ”energy” conservation law
and make possible motion with $|\varphi_{2}|\neq|\varphi_{1}|$, thus allowing
shocks of more general type than in the lossless JTL.

The r.h.s. of Eq. (V) can be written as $-d\widetilde{\Pi}(\varphi)/d\varphi$,
where

$\displaystyle\widetilde{\Pi}(\varphi)=-\cos\varphi-\frac{\sin\varphi_{2}-\sin\varphi_{1}}{\varphi_{2}-\varphi_{1}}\frac{\varphi^{2}}{2}$
$\displaystyle+\frac{\varphi_{1}\sin\varphi_{2}-\varphi_{2}\sin\varphi_{1}}{\varphi_{2}-\varphi_{1}}\varphi\,.$
(28)

Figure 7: The ”potential energy” (V). We have chosen one of
$\varphi_{1},\varphi_{2}$ to be equal to $.5$ and the other - equal to $1.0$

The ”potential energy” (V), which is graphically presented on Fig. 7, is
symmetric with respect to interchange of $\varphi_{1}$ and $\varphi_{2}$.
However, we are tempted to state that the motion with the initial and final
conditions (6) is possible only if $\varphi_{2}$ is a local maximum of the
”potential energy” $\widetilde{\Pi}(\varphi)$, and $\varphi_{1}$ is a local
minimum. Thus on Fig. 7, $\varphi_{1}=.5$ and $\varphi_{2}=1$. That
maximum/minimum distinction means that

$\displaystyle\left.\frac{d^{2}\widetilde{\Pi}(\varphi)}{d\varphi^{2}}\right|_{\varphi=\varphi_{2}}<0,\hskip
28.45274pt\left.\frac{d^{2}\widetilde{\Pi}(\varphi)}{d\varphi^{2}}\right|_{\varphi=\varphi_{1}}>0\,.$
(29)

This is certainly true both in the continuum limit and for the considered in
Section IV limiting cases, when the second derivative term was linearised.

Nor do the criteria (29) contradict to the approach of Section III based on
the properties of $\Pi(\varphi)$. In fact, the two ”potential energies” are
connected by the relation

$\displaystyle\frac{d\widetilde{\Pi}(\varphi)}{d\varphi}=\frac{d\Pi(\sin\varphi)}{d\sin\varphi}\,,$
(30)

hence their critical points coincide. Also

$\displaystyle\frac{d^{2}\widetilde{\Pi}(\varphi)}{d\varphi^{2}}=\cos\varphi\frac{d^{2}\Pi(\varphi)}{d(\sin\varphi)^{2}}\,,$
(31)

hence the second derivatives of $\Pi(\sin\varphi)$ and
$\widetilde{\Pi}(\varphi)$ have the same sign.

From Eq. (V) follows

$\displaystyle\frac{d\widetilde{\Pi}(\varphi)}{d\varphi}=\frac{(\varphi-\varphi_{1})(\varphi-\varphi_{2})}{\varphi_{2}-\varphi_{1}}\left[\overline{V}^{2}_{sh}(\varphi,\varphi_{2})-\overline{V}^{2}_{sh}(\varphi_{1},\varphi)\right]\,.$

Taking into account that

$\displaystyle
V^{2}_{sh}(\varphi_{1},\varphi_{1})=\frac{\Lambda^{2}}{L_{J}C}\cos\varphi_{1}=u^{2}(\varphi_{1})\,,$
(33)

where $u(\varphi_{1})$ is the velocity of propagation along the lossless JTL
of small amplitude smooth disturbances of $\varphi$ on a homogeneous
background $\varphi_{1}$ kogan , we can present the inequalities (29) as

$\displaystyle
u^{2}(\varphi_{1})>V^{2}_{sh}(\varphi_{1},\varphi_{2})>u^{2}(\varphi_{2})\,,$
(34)

thus establishing the connection with the well known fact: the shock velocity
is lower than the sound velocity in the region behind the shock, but higher
than the sound velocity in the region before the shock whitham .

##  VI Discussion

We hope that the results obtained in the paper are applicable to kinetic
inductance based traveling wave parametric amplifiers based on a coplanar
waveguide architecture. Onset of shock-waves in such amplifiers is an
undesirable phenomenon. Therefore, shock waves in various JTL should be
further studied, which was one of motivations of the present work.

Recently, quantum mechanical description of JTL in general and parametric
amplification in such lines in particular started to be developed, based on
quantisation techniques in terms of discrete mode operators reep , continuous
mode operators fasolo , a Hamiltonian approach in the Heisenberg and
interaction pictures greco , the quantum Langevin method yuan , or on
partitions a quantum device into compact lumped or quasi-distributed cells
minev . It would be interesting to understand in what way the results of the
present paper are changed by quantum mechanics. Particularly interesting looks
studying of quantum ripples over a semi-classical shock glazman  and fate of
quantum shock waves at late times glazman2 . Closely connected problem of
classical and quantum dispersion-free coherent propagation in waveguides and
optical fibers was studied recently in Ref. costas .

##  VII Conclusions

We have formulated the quasi-continuum approximation for the discrete JTL. The
approximation becomes controllable for the case of Josephson phase difference
across the JJs being small and for the case of small amplitude disturbances of
the phase on a homogeneous background. The quasi-continuum approximation is
applied to study the shocks and solitary waves which can propagate along the
line. The width of such waves (in the lossless JTL) turns out to be of the
order of the line period. We have found that the resistors shunting the JJs
and/or in series with the ground capacitors lead to the existence of the shock
waves of more general type than those existing in the lossless JTL. We have
suggested the generalization of the quasi-continuum approximation.

###### Acknowledgements.

The main idea of the present work was born in the discussions with M.
Goldstein. T. H. A. van der Reep read the manuscript and suggested numerous
textual corrections which were accepted with gratitude by us. We also thank B.
Ya. Shapiro for his commentaries.

##  Appendix A Quasi-continuum approximation for a linear transmission line

We consider here a linear discrete transmission line, that is the transmission
line presented on Fig. 1) modified by substituting linear inductance $L$ for
the JJ. Thus the obtained transmission line is described by

$\displaystyle
4\frac{d^{2}q_{n}}{d\tau^{2}}=\left(q_{n+1}-2q_{n}+q_{n-1}\right)\;,$ (35)

where we introduced dimensionless time $\tau$ by the equation
$t^{2}=4LC\tau^{2}$.

We start from reproducing the solution of the signalling problem for a semi-
infinite line, described by Eq. (35), presented in our previous publication
kogan . The problem is characterised by the boundary conditions
$q_{0}(\tau)=1$, $\lim_{n\to\infty}q_{n}(t)=0$, and the initial conditions
$q_{n}(0)=\dot{q}_{n}(0)=0$ for $n\geq 1$. To solve (35) we used the Laplace
transform kogan . For a given time-dependent function $f(\tau)$, we define the
Laplace transform

$Q_{n}(s)=\int_{0}^{\infty}{\rm d}\tau\,e^{-s\tau}q_{n}(\tau)\,.$ (36)

Laplace transforming (35) we obtain the difference equation

$Q_{n+1}-2(1+2s^{2})Q_{n}+Q_{n-1}=0\ ,$ (37)

with the boundary conditions $Q_{0}(s)=1/s$, $\lim_{n\to\infty}Q_{n}(s)=0$.
Solving (37):

$Q_{n}(s)=\frac{1}{s}\left[1-2s\left(\sqrt{s^{2}+1}-s\right)\right]^{n}\;,$
(38)

Taking into account the inverse Laplace transform correspondence tables abram
, we obtain

$\displaystyle\frac{dq_{n}}{d\tau}=\sum_{k=1}^{n}C^{n}_{k}(-1)^{k}2^{k}k\frac{d^{k}}{d\tau^{k}}\left(\frac{J_{k}(\tau)}{\tau}\right)\,.$
(39)

where $J_{k}$ is the Bessel function.

The solution for $dq_{n}/d\tau$ can be alternatively presented by computing
the inverse Laplace transform of (38) using the Bromwich integral

$\displaystyle\frac{dq_{n}}{d\tau}=\frac{1}{2\pi
i}\int_{a-i\infty}^{a+i\infty}ds$
$\displaystyle\exp\left\\{s\tau+n\ln\left[1-2s\left(\sqrt{s^{2}+1}-s\right)\right]\right\\}\;.$
(40)

Expanding the logarithm in (A) with respect to $s$ and keeping only the lowest
order terms we get

$\displaystyle\frac{dq_{n}}{d\tau}=\frac{1}{2\pi
i}\int_{a-i\infty}^{a+i\infty}ds\exp\left[s(\tau-2n)+\frac{n}{3}s^{3}\right]\;.$
(41)

The contour of integration in (41) can be deformed so that it will start at
the point at infinity with argument $-\pi/3$ and will end at the point at
infinity with argument $\pi/3$. Hence the integral determines the Airy
function abram

$\displaystyle\frac{dq_{n}}{d\tau}=n^{-1/3}\text{Ai}\left[n^{-1/3}(2n-\tau)\right]\;.$
(42)

Equation (42) describes the shock front at $n\sim\tau/2$, exponential decrease
of the signal with increasing $n$ for $n>\tau/2$, and oscillations and power
law decrease of the signal with decreasing $n$ for $n<\tau/2$.

Now let us solve the signalling problem for Eq. (35) approximately. In the
continuum approximation we obtain the equation

$\displaystyle\frac{\partial^{2}q}{\partial
z^{2}}=\frac{4}{\Lambda^{2}}\frac{\partial^{2}q}{\partial\tau^{2}}\,.$ (43)

The solution of Eq. (43) with the boundary conditions $q(0,\tau)=1$,
$\lim_{z\to\infty}q_{(}z,t)=0$ and the initial conditions
$q(z,0)=\dot{q}(z,0)=0$, is obvious:

$\displaystyle\frac{\partial
q(z,\tau)}{\partial\tau}=\delta(\tau-2z/\Lambda)\,.$ (44)

It is also obvious, that the solution (44) completely misses the structure of
the exact solution.

In the quasi-continuum approximation we modify Eq. (43) to

$\displaystyle\frac{\partial^{2}q}{\partial
z^{2}}+\frac{\Lambda^{2}}{12}\frac{\partial^{4}q}{\partial
z^{4}}=\frac{4}{\Lambda^{2}}\frac{\partial^{2}q}{\partial\tau^{2}}\,,$ (45)

keeping the same boundary condition initial conditions as above. To solve the
signalling problem we will again use the Laplace transform, obtaining the
equation for $Q(z,s)$

$\frac{d^{2}{\cal
Q}}{dz^{2}}+\frac{\Lambda^{2}}{12}\frac{d^{4}Q}{dz^{4}}=\frac{4s^{2}Q}{\Lambda^{2}}\,,$
(46)

with the boundary conditions $Q(0,s)=1/s$, $\lim_{z\to\infty}Q(z,s)=0$. Taking
into account the boundary conditions we obtain the solution of Eq. (46) as

$Q(z,s)=\frac{1}{s}e^{kz}\,.$ (47)

where $k$ is the negative real root of the polynomial equation

$\displaystyle
k^{2}+\frac{\Lambda^{2}}{12}k^{4}=\frac{4s^{2}}{\Lambda^{2}}\,.$ (48)

Computing the Bromwich integral we get

$\displaystyle\frac{\partial q(z,\tau)}{\partial\tau}=\frac{1}{2\pi
i}\int_{a-i\infty}^{a+i\infty}dse^{s\tau+k(s)z}\;.$ (49)

Solving Eq. (48) for $s/\Lambda\ll 1$

$\displaystyle k=-2\frac{s}{\Lambda}+\frac{s^{3}}{3\Lambda}\,,$ (50)

we obtain

$\displaystyle\frac{\partial q(z,\tau)}{\partial\tau}=\frac{1}{2\pi
i}\int_{a-i\infty}^{a+i\infty}ds\exp\left[s\left(\tau-\frac{2z}{\Lambda}\right)+\frac{z}{\Lambda}\frac{s^{3}}{3}\right]\;,$

which coincides with Eq. (41).

##  Appendix B The integral quasi-continuum approximation

In an attempt to generalize the approach of Section II, let us, after
upgrading $\varphi_{n}$ to the function of the continuous variable $z$,
approximate the r.h.s. of Eq. (2) as

$\displaystyle\text{r.h.s.}=\Lambda^{2}\int_{-\infty}^{+\infty}dz^{\prime}f(z-z^{\prime})\frac{d^{2}\sin\varphi(z^{\prime})}{dz^{\prime
2}}\,,$ (52)

where $f(z)$ is some given positive even function with zero moment equal to
one

$\displaystyle\int_{-\infty}^{+\infty}dzf(z)=1\,.$ (53)

The approximation being done, Eq. (2) is reduced to (compare with (4))

$\displaystyle\frac{L_{J}C}{\Lambda^{2}}\frac{\partial^{2}\varphi(z,t)}{\partial
t^{2}}=\Lambda^{2}\int_{-\infty}^{+\infty}dz^{\prime}f(z-z^{\prime})\frac{d^{2}\sin\varphi(z^{\prime},t)}{dz^{\prime
2}}\,.$

Looking for the solution of the equation in the form
$\varphi(z,t)=\varphi(x)$, where $x=Vt-z$, we obtain the integro-differential
equation (the counterpart of (5)) for the function $\varphi(x)$

$\displaystyle\overline{V}^{2}\frac{d^{2}\varphi(x)}{dx^{2}}=\int_{-\infty}^{+\infty}dx^{\prime}f(x-x^{\prime})\frac{d^{2}\sin\varphi(x^{\prime})}{dx^{\prime
2}}\,.$ (55)

Imposing the boundary conditions (6) and integrating Eq. (55) with respect to
$x$ twice we obtain the counterpart of (7)

$\displaystyle\overline{V}^{2}\varphi(x)=\int_{-\infty}^{+\infty}dx^{\prime}f(x-x^{\prime})\sin\varphi(x^{\prime})-F\,.$
(56)

Going in Eq. (56) to the limits $x\to+\infty$ and $x\to-\infty$, we recover
Eq. (8). Thus for the shock wave with $\varphi_{2}=-\varphi_{1}$ we recover
Eq. (12). The equation for the shock wave velocity turns out to be more
general than the quasi-continuum approximation we used in Section II to derive
it. Actually, this fact could have been anticipated, because Eq. (12b) doesn’t
contain the ”microscopic” parameter $\Lambda$. Finally, Eq. (56) takes the
form of Hammerstein equation of the second kind tricomi

$\displaystyle\frac{\sin\varphi_{1}}{\varphi_{1}}\varphi(x)=\int_{-\infty}^{+\infty}dx^{\prime}f(x-x^{\prime})\sin\varphi(x^{\prime})\,.$
(57)

This equation is the counterpart of Eq. (7).

Detailed study of Eq. (57) will be presented elsewhere, but a few features of
the solution can be seen immediately. (We consider now Eq. (57) per se,
forgetting the properties of $\varphi(x)$ which were used to derive the
equation.)

First, for any $x$,

$\displaystyle-\varphi_{1}\leq\varphi(x)\leq\varphi_{1}$ (58)

(everywhere in our paper $\varphi\in(-\pi/2,\pi/2)$; here for the sake of
definiteness we consider $\varphi_{1}$ to be positive). In fact, let
$\sin\varphi(x)$ reaches maximum value at some point $x_{0}$. Then

$\displaystyle\int_{-\infty}^{+\infty}dx^{\prime}f(x-x^{\prime})\sin\varphi(x^{\prime})<\sin\varphi(x_{0})\,.$
(59)

Let also $\varphi(x_{0})>\varphi_{1}$. Then

$\displaystyle\frac{\sin\varphi_{1}}{\varphi_{1}}\sin\varphi(x_{0})>\sin\varphi(x_{0})\,,$
(60)

because $\sin\varphi/\varphi$ decreases when $\varphi$ increases for positive
$\varphi$. So we came to a contradiction. Similar for the minimum value of
$\sin\varphi$.

Second, if $\varphi(x)$ goes to some limits when $x\to+\infty$ and
$x\to-\infty$, each of these limits is either $\varphi_{1}$ or $-\varphi_{1}$
(this obviously follows from Eq. (53)). Because of (58), one of the limits is
$-\varphi_{1}$, and the other $\varphi_{1}$.

Third, for $\varphi_{1}\sim 1$, the width of the region where the solution
effectively changes from one asymptotic value to the other is of the order of
$\Lambda$ (there are no other scales in the problem). This justifies the term
shock wave in the context.

It would probably be quite difficult to find the solution (15) just by
inspection of Eq. (57), but with the wisdom after the event we can check up
that (15) really solves (57) for $\varphi_{1}\ll 1$, if we impose on $f$, in
addition to (53), the condition

$\displaystyle\int_{-\infty}^{+\infty}dxx^{2}f(x)=\frac{\Lambda^{2}}{6}\,.$
(61)

In fact, we can verify the equality

$\displaystyle\varphi_{1}\left(1-\frac{\varphi_{1}^{2}}{6}\right)\tanh\frac{\varphi_{1}x}{\Lambda}=\int_{-\infty}^{+\infty}dx^{\prime}\Big{[}\delta(x-x^{\prime})$
$\displaystyle+\frac{\Lambda^{2}}{12}\delta^{\prime\prime}(x-x^{\prime})\Big{]}\left(\varphi_{1}\tanh\frac{\varphi_{1}x}{\Lambda}-\frac{\varphi_{1}^{3}}{6}\tanh^{3}\frac{\varphi_{1}x}{\Lambda}\right)$

to the orders $\varphi_{1}$ and $\varphi_{1}^{3}$. (While approximating Eq.
(57), we took into account (53) and (LABEL:moments) defining $f$, and we
expanded the sine functions in both sides of the equation up to the third
order terms.)

For $\varphi_{1}\ll 1$, the integral quasi-continuum approximation (52) (and
following from it Eq. (57)) look like an awkward way to arrive to the result,
which we obtained much easier previously. It is possible, however, that in the
general case, where the approximation (4) is uncontrollable, the integral
quasi-continuum approximation and following from it equations give better
results than (4) and (10). The question, how to chose $f(z)$ in this case,
should be additionally addressed.

## References

  * (1) G. B. Whitham, Linear and Nonlinear Waves, John Wiley & Sons Inc., New York (1999).
  * (2) D. M. French and B. W. Hoff, IEEE Trans. Plasma Sci. 42, 3387 (2014).
  * (3) B. Nouri, M. S. Nakhla, and R. Achar, IEEE Trans. Microw. Theory Techn. 65, 673 (2017).
  * (4) L. P. S. Neto, J. O. Rossi, J. J. Barroso, and E. Schamiloglu, IEEE Trans. Plasma Sci. 46, 3648 (2018).
  * (5) M. S. Nikoo, S. M.-A. Hashemi, and F. Farzaneh, IEEE Trans. Microw. Theory Techn. 66, 3234 (2018); 66, 4757 (2018).
  * (6) L. C. Silva, J. O. Rossi, E. G. L. Rangel, L. R. Raimundi, and E. Schamiloglu, Int. J. Adv. Eng. Res. Sci. 5, 121 (2018).
  * (7) Y. Wang, L.-J. Lang, C. H. Lee, B. Zhang, and Y. D. Chong, Nat. Comm. 10, 1102 (2019).
  * (8) E. G. L. Range, J. O. Rossi, J. J. Barroso, F. S. Yamasaki, and E. Schamiloglu, IEEE Trans. Plasma Sci. 47, 1000 (2019).
  * (9) A. S. Kyuregyan, Semiconductors 53, 511 (2019).
  * (10) N. A. Akem, A. M. Dikande, and B. Z. Essimbi, Social Netw. Appl. Sci. 2, 21 (2020).
  * (11) A. J. Fairbanks, A. M. Darr, A. L. Garner, IEEE Access 8, 148606 (2020).
  * (12) R. Landauer, IBM J. Res. Develop. 4, 391 (1960).
  * (13) S. T. Peng and R. Landauer, IBM J. Res. Develop. 17(1973).
  * (14) M. I. Rabinovich and D. I. Trubetskov, Oscillations and Waves, Kluwer Academic Publishers, Dordrecht / Boston / London (1989).
  * (15) B. D. Josephson, Phys. Rev. Lett. 1, 251 (1962).
  * (16) A. Barone and G. Paterno, Physics and Applications of the Josephson Effect, John Wiley & Sons, Inc, New York (1982).
  * (17) N. F. Pedersen, Solitons in Josephson Transmission lines, in Solitons, North-Holland Physics Publishing, Amsterda (1986).
  * (18) C. Giovanella and M. Tinkham, Macroscopic Quantum Phenomena and Coherence in Superconducting Networks, World Scientific, Frascati (1995).
  * (19) A. M. Kadin, Introduction to Superconducting Circuits, Wiley and Sons, New York (1999).
  * (20) M. Remoissenet, Waves Called Solitons: Concepts and Experiments, Springer-Verlag Berlin Heidelberg GmbH (1996).
  * (21) O. Yaakobi, L. Friedland, C. Macklin, and I. Siddiqi, Phys. Rev. B 87, 144301 (2013).
  * (22) K. O’Brien, C. Macklin, I. Siddiqi, and X. Zhang, Phys. Rev. Lett. 113, 157001 (2014).
  * (23) C. Macklin, K. O’Brien, D. Hover, M. E. Schwartz, V. Bolkhovsky, X. Zhang, W. D. Oliver, and I. Siddiqi, Science 350, 307 (2015).
  * (24) B. A. Kochetov, and A. Fedorov, Phys. Rev. B. 92, 224304 (2015).
  * (25) A. B. Zorin, Phys. Rev. Applied 6, 034006 (2016); Phys. Rev. Applied 12, 044051 (2019).
  * (26) D. M. Basko, F. Pfeiffer, P. Adamus, M. Holzmann, and F. W. J. Hekking, Phys. Rev. B 101, 024518 (2020).
  * (27) T. Dixon, J. W. Dunstan, G. B. Long, J. M. Williams, Ph. J. Meeson, C. D. Shelly, Phys. Rev. Applied 14, 034058 (2020)
  * (28) A. Burshtein, R. Kuzmin, V. E. Manucharyan, and M. Goldstein, Phys. Rev. Lett. 126, 137701 (2021).
  * (29) T. C. White et al., Appl. Phys. Lett. 106, 242601 (2015).
  * (30) A. Miano and O. A. Mukhanov, IEEE Trans. Appl. Supercond. 29, 1501706 (2019).
  * (31) Ch. Liu, Tzu-Chiao Chien, M. Hatridge, D. Pekker, Phys. Rev. A 101, 042323 (2020).
  * (32) E. Kogan, Journal of Applied Physics 130, 013903 (2021).
  * (33) H. R. Mohebbi and A. H. Majedi, IEEE Trans. Appl. Supercond. 19, 891 (2009); IEEE Transactions on Microwave Theory and Techniques 57, 1865 (2009).
  * (34) G. J. Chen and M. R. Beasley, IEEE Trans. Appl. Supercond. 1, 140 (1991).
  * (35) H. Katayama, N. Hatakenaka, and T. Fujii, Phys. Rev. D 102, 086018 (2020).
  * (36) T. H. A. van der Reep, Phys. Rev. A 99, 063838 (2019).
  * (37) L Fasolo, A Greco, E Enrico, in Advances in Condensed-Matter and Materials Physics: Rudimentary Research to Topical Technology, (ed. J. Thirumalan and S. I. Pokutny), Sceence (2019).
  * (38) A. Greco, L. Fasolo, A. Meda, L. Callegaro, and E. Enrico, Phys. Rev. B 104, 184517 (2021).
  * (39) Y. Yuan, M. Haider, J. A. Russer, P. Russer and C. Jirauschek, 2020 XXXIIIrd General Assembly and Scientific Symposium of the International Union of Radio Science, Rome, Italy (2020).
  * (40) Z. K. Minev, Th. G. McConkey, M. Takita, A. D. Corcoles, and J. M. Gambetta, arXiv:2103.10344 (2021).
  * (41) E. Bettelheim and L. I. Glazman, Phys. Rev. Lett. 109, 260602 (2012).
  * (42) Th. Veness and L. I. Glazman, Phys. Rev. B 100, 235125 (2019).
  * (43) A. Mandilara, C. Valagiannopoulos, and V. M. Akulin, Phys. Rev. A99, 023849 (2019).
  * (44) M. Abramowitz, I. A. Stegun eds., Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, (National Bureau of Standards, Washington, 1964).
  * (45) F. G. Tricomi, Integral equations, (Mineola, NY: DoverPublications, 1982).