\section{Introduction}
Light-by-light scattering is a purely quantum process, hence, in the standard model,
only possible in the context of quantum electrodynamics (QED).
QED is the most strictly tested dynamics in its perturbative regime.
The scattering amplitude in the lowest order is described by the box diagram
where a virtual $e^+e^-$ pair loops between four external photons.
So far the diagram is indirectly confirmed with off-shell external photons via
$\gamma A\rightarrow \gamma A$~\cite{DelbrueckExp}
known as Delbr\"uck scattering~\cite{Delbrueck} and
also proposed to be further tested via $A A \rightarrow A A\gamma \gamma$
in the higher energy domain at the Large Hadron Collider~\cite{LHCLbyL}.
In these scatterings, however, photons are emitted from nuclei $A$ and
have finite off-shell masses. With the exact on-shell condition, namely,
real-photon - real-photon scattering has not been demonstrated hitherto
despite of its explicit predictions~\cite{Halpern:1934,Toll:1952,R5,R6,R7}.
Moreover, its dependence on the photon polarization states
has not been tested at all. By means of off-shell incident photons from nuclei,
it is difficult to test the polarization property.
Photon-photon interactions also contain rich information on
the two-photon coupling to standard model / non-standard model
resonance states depending on the center of mass energy.
In the MeV range, no direct search for resonance states
coupling to two photons has been performed.
Therefore, it is indispensable to first verify the purely QED-based
scattering amplitude
with specified polarization states in a pristine experimental condition,
because a significant deviation from the prediction indicates
the existence of physics beyond the standard model.

\begin{figure}
\includegraphics[width=0.7\linewidth]{Fig1.eps}
\caption{
Helicity configuration in $\gamma\gamma$ scattering via inverse Compton
scatterings between cicularly polarized laser pulses and
unpolarized electrons.
R and L denote right- and left-handed helicity
states of photons, respectively.
}
\label{Fig0}
\end{figure}

In lower energies, searches for photon-photon
interactions have been performed at keV~\cite{R8}, eV~\cite{R9,R10} and
sub-eV energy scales attempting to search for dark matter
~\cite{R11,R12,R13,R14,R15,R16,R17,R18,R19} and also
to put the upper limit on the QED photon-photon interaction~\cite{R20}.
There are several proposals~\cite{Lundstrom,Heinzl:2006xc,R21,R22}
to probe photon-photon interactions
in the context of the low energy limit of the QED interactions
~\cite{Euler:1935zz,Heisenberg:1935qt} based on only optical photons.
These proposals commonly suffer from its
extremely weak interaction. Therefore, currently state-of-the-art high-power
lasers are supposed to be indispensable. However, if such lasers are available,
we may also use them to generate $\gamma$-rays. This idea leads us to try to
consider the direct measurement of the scattering with the maximized
cross-section at the center-of-mass system energy of 1-2 MeV
rather than based on 1 eV photons.

In this paper we discuss the feasibility to measure
the helicity dependent $\gamma\gamma$-scattering
based on the helicity configuration as illustrated in Fig.\ref{Fig0}
where circularly polarized laser pulses are reflected upon unpolarized
electron bunches and $\gamma$-rays in the same helicity state collide head-on.
We then consider to measure the elastic scattering
without specifying the helicity states of the final state $\gamma$-rays.
In contrast to the helicity dependence emphasizing on the forward scattering
amplitude~\cite{Dinu:2013gaa,Dinu:2014tsa}, we rather aim
at larger angle scattering events in order to verify
the purely perturbative prediction. In the following paragraphs, we discuss
the QED-based helicity dependent $\gamma\gamma$-scattering cross section,
how to produce electrons with LPA, the expected helicity specified
$\gamma$-ray yield based on the nonlinear Compton process and then
possibility to realize a realistic experiment.

\begin{figure}
\includegraphics[width=0.7\linewidth]{Fig2.eps}
\caption{ The unpolarized QED-based differential cross section
$d\sigma_{qed}/d\Omega$ at $\vartheta = \pi/2$
in units of $\mu$b as a function of
the center-of-mass system energy in photon-photon collisions,
$E_{cms}$.
}
\label{Fig1}
\end{figure}

\begin{figure}
\includegraphics[width=0.6\linewidth]{Fig3.eps}
\caption{
Helicity dependence of QED-based differential cross sections
at $E_{cms}=1.4$ MeV as a function of the photon scattering angle $\vartheta$
in CMS.
}
\label{Fig1b}
\end{figure}

\section{QED-based photon-photon interaction}
The elastic scattering of photons by photons has been thoroughly evaluated
in \cite{R5,R6,R7}.
The differential cross section per solid angle
for the unpolarized photon-photon scattering
in the lower energy limit
is approximated as
$
\frac{d\sigma_{qed}}{d\Omega} \sim \frac{\alpha^2 {r_0}^2}{4\pi^2}
\frac{139}{(90)^2}\left(\frac{k}{m}\right)^6 (3+\cos\vartheta)^2,
$
where
$m$ is electron mass in units of $\hbar=c=1$,
$k$ is the incident photon energy in the center-of-mass system,
$\alpha = 1/137$ is the dimensionless fine structure constant,
$r_0 = \alpha/m \sim 2.8 \times 10^{-13}$~cm is
the classical electron radius, and
$\vartheta$ is the polar angle of scattered photons with respect to
the colliding axis between two incident photons in the center-of-mass
system on the reaction plane.
For laser photons of $k \sim 1$~eV, the total cross section is
$\sigma \sim 10^{-42}$~b.
This is extremely small due to the steep $k^6$ dependence.
This situation is shown in Fig.\ref{Fig1} where the differential
cross section of the unpolarized photon-photon elastic scattering,
$d\sigma_{qed}/d\Omega$ at $\vartheta = \pi/2$, is plotted as a function of
the center-of-mass system energy of photon-photon collisions, $E_{cms}$.
The cross section curve is obtained from the numerical
calculation applying formulae in Ref.~\cite{R7}.
On the other hand, if $k = 0.5-1.0$~MeV is realized,
the total cross section is maximized up to $\sigma \sim 1$~$\mu$b.
Therefore, in order to detect the perturbative QED
cross section, it is reasonable to perform the scattering experiment
tuned at that energy range.
In the region $E_{cms} = 1-2$~MeV, we see a flat-top character.
This allows relatively large fluctuations on $E_{cms}$, which is preferable
for the photon-photon collider exploiting $\gamma$-rays
via laser Compton scattering off electrons from LPA having a percent-level energy spread.

Figure \ref{Fig1b} shows the helicity dependent QED-based differential cross section
$d\sigma^{\pm \rightarrow any}/d\theta$ compared to the unpolarized cross section
where the symbols $\pm \rightarrow any$ indicate the same and opposite
helicity states between incident two $\gamma$-rays
going to all the possible helicity states.
These predictions are results of numerical calculations at $E_{cms}=1.4$ MeV
based on Ref.~\cite{R7}.
We find a larger cross section in the same helicity incidence case
particularly around $\vartheta \sim \pi/2$.

\begin{figure}
\includegraphics[width=1.0\linewidth]{Fig4.eps}
\caption{
An all-optical table-top (3.4 m $\times$ 1.3 m) $\gamma\gamma$-collider:
a) top-view including two LPAs and the detector system to capture
the $\gamma\gamma\rightarrow\gamma\gamma$ scattering
and
b) collision geometry around the interaction point, IP,
where $\gamma$-rays are produced at each Compton scattering point (CP)
in head-on collisions and $D$ is the distance between IP and CP.
}
\label{Fig2}
\end{figure}

\section{Table-top $\gamma\gamma$-collider}
Figure \ref{Fig2} a) shows a table-top (3.4 m $\times$ 1.3 m) $\gamma\gamma$-collider
with a realistic detector system to capture large angle
$\gamma\gamma\rightarrow\gamma\gamma$ scattering events.
The system consists of two LPAs to generate electron beams with
which two incident $\gamma$-beams are further produced via the inverse
Compton process in head-on collisions. Initially synchronized two laser pulses
are incident from the top and bottom sides of the top view, respectively.
They are individually split into two drive pulses for LPA and scatter pulses
for the successive inverse Compton process with electrons delivered from LPA.
Figure \ref{Fig2} b) illustrates the collision geometry around the interaction point (IP)
of $\gamma\gamma$-scattering, which is located at the distance $D$
from the inverse Compton scattering point (CP).

\begin{figure}
\includegraphics[width=1.0\linewidth]{Fig5.eps}
\caption{
A focusing system for an electron beam produced by LPA consisting of
three PMQs from the ejection point (EP) of the electron beam at
the gas cell for LPA to the inverse Compton scattering point (CP).
The blue and red curves are traces of beam envelopes simulated
by TRACE3D~\cite{R30} for horizontal and vertical directions, respectively.
}
\label{Fig3}
\end{figure}

An electron beam with 210 MeV and 1.6 nC are produced from
a two-stage laser wakefield accelerator~\cite{R23} comprising a 5-mm long gas
cell filled with the mixed gas (e.g., 94 \% helium and 6 \% nitrogen) for
the injector stage and a variable-length gas cell filled with pure helium
for the accelerator stage.
Designing parameters of the laser wakefield accelerator can be carried out by
relying on the scaling law of nonlinear plasma wakefields
in the bubble regime~\cite{R24,R25,R26}.
Provided that a laser pulse with 41 TW peak power and 85 fs duration is
focused on 12 $\mu$m spot radius on the entrance of the injector cell operated
at plasma density of $3.3\times 10^{18}$ cm${}^{-3}$,
strong nonlinear wakefields can be generated so that a 1.6 nC electron bunch
could be trapped due to ionization-induced injection~\cite{R27,R28}
and accelerated up to 40 MeV, followed by boosting its energy up to 210 MeV
at the length of 2.6 cm in the accelerator cell operated at plasma density
of $1.1 \times 10^{18}$  cm${}^{-3}$. The relative energy spread and
the normalized emittance of resultant output beams are estimated to
be 4 \% in r.m.s. and 0.15 mm mrad, respectively.

The electron bunch is then focused via a set of permanent-magnet-based
quadrupoles (PMQs)~\cite{R29} consisting of three elements
-275 T/m, 770 T/m and -650 T/m with the common inner radius 3.0 mm and
outer radius 12.0 mm, respectively, over 71.2 cm
as shown in Fig.\ref{Fig3}
which displays variations of the horizontal and vertical beam envelopes
simulated by TRACE3D~\cite{R30}
for the given incident parameters of the LPA configuration above.

Based on these counter-propagating electron beams,
we evaluate the number of incident $\gamma$-rays produced by the inverse nonlinear
Compton process between a circularly polarized laser pulse and an unpolarized
electron bunch produced by LPA.
It is worth noting that once we fix the focusing geometry,
we cannot increase the effective number of $\gamma$-rays as much as we like
even if we could increase the laser pulse energy. This is because of
the nonlinear nature of Compton scattering represented by the parameter
$\eta \equiv e\sqrt{-\langle A_{\mu} A^{\mu}\rangle}/mc^2$ where $A_{\mu}$ is
the four-vector potential of the laser pulse~\cite{SLACDron}.
The $\eta$ parameter increases the effective electron mass $m^*=m\sqrt{1+\eta^2}$
in the laser field. Hence, the scattered photon energy in the single photon absorption case
is effectively lowered and the photon yield in the Compton edge energy also diminishes.
On the other hand, the larger number of laser photons increases the
luminosity factor in the Compton scattering.
Therefore, there is an optimal $\eta$ value so that the scattered
photon energy is kept within 0.5-1.0 MeV range
with the maximized $\gamma$-ray yield.

In Fig.\ref{Fig4}, we plot this situation in three different $\eta$ cases:
solid ($\eta=0.31$), dashed ($\eta=0.63$) and dash-dotted ($\eta=0.88$).
Figure \ref{Fig4} a) shows the differential cross sections
$d\sigma^{(\pm)}_{el}/d\theta$ in units of $(\mu\mbox{m})^2/rad$
as a function of photon scattering
angle measured from the incident direction of the electron beam,
where the abscissa is displayed in units of $1/\gamma^*$
with the Lorentz factor defined by
the effective electron mass in the laser field.
The symbols $(\pm)$ denote circular polarization flip $(-)$ and non-flip $(+)$
cases in the transition between the initial and final state photons.
The $(-)$ and $(+)$ cases are plotted with thicker and thiner lines, respectively.
The cross sections are numerically calculated based on the expressions
available in Ref.~\cite{SLACDron}.
In any $\eta$ cases, in order to enhance the purity of the $(-)$ case,
we must consider $\gamma$-rays produced only in $\theta < 1/\gamma^*$
as the effective number of $\gamma$-rays useful for the test of helicity specified
$\gamma\gamma$-scatterings.
Figure \ref{Fig4} b) shows $E_{cms}$ distributions by taking all possible
energy combinations between two incident $\gamma$-rays from the Compton scattering
points with the cross section weights in the $(-)$ case
within $\theta < 1/\gamma^*$ in Fig.\ref{Fig4} a),
where the area of the distributions indicate squares of the numbers of generated
$\gamma$-rays, $N_{\gamma}$, corresponding to the numerators in the luminosity factor
for the $\gamma\gamma$-collider as we discuss below.

The head-on luminosity factor for
$\gamma\gamma$-scattering is defined as
$
\frac{f N^2_{\gamma}}{4\pi r^2_b}
$
where we assume $\gamma\gamma$-collisions take place along the electron
beam axis within the effective $\gamma$-beam radius at IP, $r_b$.
Because the distance between the two CP points, $D$, should be kept as small as possible
to increase $N_{\gamma}$,
$r_b$ effectively coincides with the beam waist of the Compton seed laser,
$w_0$, if the electron beam radius is much smaller than $w_0$.
On the other hand, in order to enrich the $+ \rightarrow any$ case
in the QED-based $\gamma\gamma$-scattering,
the requirement to enhance the $(-)$ case
by limiting the effective scattering angle $\theta < 1/\gamma^*$ demands
the relation $D/\gamma^*=w_0$.
This implies that we can control the purity of the initial circular polarization states
by adjusting $D$ in experiments.

\begin{figure}
\includegraphics[width=0.8\linewidth]{Fig6.eps}
\caption{
Nonlinear effects in Compton scattering
with $\eta = 0.31$(solid), 0.63(dashed) and 0.88(dash-dotted)
:
a) differential Compton scattering cross sections
$\sigma^{(\pm)}_{el}/d\theta$ as a function of photon scattering
angle measured from the incident direction of the electron beam
in units of $1/\gamma^*$ with the Lorentz factor of the
effective electron mass in the laser field.
The symbols $(\pm)$ denote circular polarization flip $(-)$ and non-flip $(+)$
cases in the transition between the initial and final state photons.
b) $E_{cms}$ distributions by taking all possible
energy combinations between two incident $\gamma$-rays from the Compton scattering
points with the cross section weights in the case of $(-)$ in $\theta < 1/\gamma^*$,
where the distributions are normalized to squares of
the numbers of generated $\gamma$-rays.
}
\label{Fig4}
\end{figure}

\section{Design of the detection system}
We summarize optimal parameters for an all-optical photon-photon collider
in the case of $\eta = 0.63$ in Table 1.
For instance, if we apply this detection system to two
synchronized lasers operated at $f\sim10$ Hz, a data taking period over
several months will be sufficient to claim the statistical significance
of the QED-based elastic scattering events for the $+ \rightarrow any$ case.

\renewcommand{\arraystretch}{0.9}
\begin{table}[h!]
\caption{
Parameters for the all-optical photon-photon collider LPA drive laser
}
\begin{center}
\begin{tabular}{lr}
\\ \hline
LPA drive laser
\\ \hline
Wavelength [$\mu$m]	& 0.8\\
Repetition rate [Hz] & 	10\\
Pulse energy [J] & 3.5\\
Peak power [TW] & 41\\
Pulse duration [fs]&85
\\ \hline
LPA electron beam
\\ \hline
Beam energy [MeV]&210\\
Plasma density [$10^{18}$ cm${}^{-3}$]&1.1\\
Accelerator length [cm]&2.6\\
Charge per bunch [nC]&1.6\\
Bunch duration [fs]& $\sim 10$\\
Normalized emittance [mm mrad]&	$\sim 0.15$\\
Horizontal rms beam size at Compton IP [$\mu$m]&0.8\\
Vertical rms beam size at Compton IP [$\mu$m] &0.5
\\ \hline
Laser for Compton scattering
\\ \hline
Wavelength [$\mu$m]&0.8\\
Repetition rate [Hz]&10\\
Pulse energy [mJ]&89\\
Pulse duration [fs]&209\\
Spot radius at Compton IP [$\mu$m]&4.0\\
Interaction angle [degree]&0\\
Focused intensity [$10^{17}$ W/cm${}^{2}$]&8.5
\\ \hline
Compton $\gamma$-ray beam
\\ \hline
Photon energy between $\theta = 0 - 1/\gamma^*$ rad [MeV]&0.37-0.75\\
Effective photon flux [$10^9$ s${}^{-1}$]&8.65\\
Effective spot radius $r_b$ at $\gamma\gamma$-IP [$\mu$m]&4.0
\\ \hline
$\gamma\gamma$-scattering
\\ \hline
Averaged CMS photon-photon collision energy [MeV]&1.5\\
Averaged QED cross section [$\mu$b]&3.33\\
Rate of QED events [$10^{-6}$ s${}^{-1}$]&2.4\\
\hline
\end{tabular}
\end{center}
\label{Tab1}
\end{table}
\renewcommand{\arraystretch}{1.0}

Due to the short $D$, we simultaneously have
to consider electron-electron scattering after the inverse Compton scattering
occurs,
because this M{\o}ller's scattering can produce dominant background events against
the elastic $\gamma\gamma$-scattering as shown in Fig.\ref{Fig7} a) and b).
However, as far as only large angle scattering events in the range of
$\vartheta=45-135$ degree are measured, the partially integrated cross section
of M{\o}ller's scattering over that solid angle is suppressed to 17.7 $\mu$b
for the electron energy of 210 MeV, which is evaluated by the following
differential cross section with respect to solid angle
$
\frac{d\sigma_M}{d\Omega} = \frac{\alpha^2}{4E^2(E^2-m^2)^2}
\left[\frac{4(2E^2-m^2)^2}{\sin^4\vartheta}
-\frac{8E^4-4E^2m^2-m^4}{\sin^2\vartheta}+(E^2-m^2)^2 \right]
$
where $E$ is the incident electron energy in the center-of-mass system.
On the other hand, the cross section of QED $\gamma\gamma$-scattering is
relatively enhanced over the background events in the same $\vartheta$ range.
Taking the higher luminosity factor for the electron-electron scattering
than that of $\gamma\gamma$-scattering into account,
the event rate for M{\o}ller's scattering reaches 0.36 Hz for $f=10$ Hz and
the accidental rate for two types of scattering events to occur within
the same shot can be evaluated as $2.6\times10^{-6}$ Hz.
This contaminated event rate
corresponds to 36 \% of that of the QED $\gamma\gamma$-scattering.
Therefore, even if one throws such contaminated events away without
the detailed event-by-event offline analysis,
the statistical loss of the QED events is still acceptable.
In reality, the energy deposits on the detector
as well as the event topologies, whether electromagnetic showers exist or not,
are very different between QED $\gamma\gamma$-scattering and M{\o}ller's
scattering as displayed in Fig.\ref{Fig7} a) and b).
Therefore, one can readily distinguish two types of scattering events
at the offline analysis and possibly distinguish them even within
the same shot, if the two $\gamma$-clusters are sufficiently
isolated from the background electromagnetic showers.
Moreover, counting the number of M{\o}ller's scattering events is
indispensable to directly measure the electron-electron luminosity
in the actual experimental condition. This information is essentially important
to deduce the $\gamma\gamma$ luminosity in addition to
the direct measurement of the $\gamma$-ray flux.

\begin{figure}
\includegraphics[width=1.0\linewidth]{Fig7.eps}
\caption{
Comparison of event displays simulated by GEANT4 for
a) a QED-based $\gamma\gamma\rightarrow\gamma\gamma$ event and
b) a background $e^- e^- \rightarrow e^- e^-$ event, where
the blue and black trajectories denote photons and
electrons, respectively. The detector system covers the polar angle
$\vartheta =45-135$ degree and the azimuthal angle $\varphi=0-360$ degree
around the beam axis consisting of 90 scintillator crystals
(18 crystals in $\varphi \times 5$ layers in $\vartheta$)
made of Ce doped Gd${}_2$SO${}_5$ (GSO:Ce).
GSO:Ce is a well-balanced scintillator from the point of view
on the scintillation-photon yield for a sub-MeV $\gamma$-ray,
the time resolution and the radiation length $X_0$.
The radiation length $X_0=1.38$ cm is reasonably small to
suppress the lateral spreads of the electromagnetic showers produced
by 210 MeV electrons. The individual crystal has the depth of $10 X_0$.
}
\label{Fig7}
\end{figure}

\section{Conclusion}
Verification of the QED-based helicity dependent
light-by-light scattering
with the all-optical photon-photon collider is realizable
in the table-top scale.
Provided PW-class lasers such as Extreme Light Infrastructure (ELI)~\cite{R31}
we may be able to test the QED-based light-by-light scattering especially in
the same circular polarization case of incident $\gamma$-rays
with an experimentally feasible rate.
Furthermore, the rapid development of fiber-based high-power and high-repetition
rate lasers at 10~kHz~\cite{R32} would allow us to test the opposite circular polarization
case and enable the comparison between the two helicity configurations
quantitatively.
The larger statistics also opens up opportunities to perform general
resonance searches
around the MeV energy scale based on the landmark of quantitative
verification of the QED-based perturbative photon-photon scattering.
Finally we emphasize that our proposal corresponds to the first case of LPA-based
electron-electron and $\gamma\gamma$-colliders applied to
fundamental particle physics.
If succeeded, it would be a ground-breaking advancement
in experimental particle physics
even if the center of mass energy is still in the MeV range.

\begin{acknowledgments}
We cordially thank Y. Iwashita for deeply considering the possibility to
introduce the permanent-magnet-based dipole at IP in the electron transport
line to avoid M{\o}ller's scattering in advance.
K. Homma appreciates E. Milotti and C. Curceanu for discussions on
the derivation for the QED-based cross section.
He also acknowledges for the support by the Collaborative Research
Program of the Institute for Chemical Research,
Kyoto University (grants No. 2015-93).
K. Nakajima was supported by the Project Code (IBS-R012-D1).
\end{acknowledgments}


\section{Introduction}
Strong stability preserving Runge--Kutta (RK) methods were first introduced by
Shu and Osher \cite{shu1988} in the context of time integration for first-order hyperbolic
conservation laws:
\begin{align} \label{pde}
	{\mathcal U}_t + {\mathcal F}({\mathcal U})_x & = 0, &
    {\mathcal U}(x,t=0) = {\mathcal U}_0.
\end{align}
After semi-discretization, \eqref{pde} takes the form of an initial-value
ordinary differential equation system:
\begin{align}
u'(t) & = f(u) & u(0) = u_0\, , \label{ode}
\end{align}
where $f$ is a discrete approximation to $-\mathcal F$.
In the scalar case $\mathcal U$ is dissipative, and it is natural to seek a
semi-discretization that is dissipative:
\begin{align} \label{dissipative}
\frac{d}{dt}\|u\| \le 0,
\end{align}
where $\|\cdot\|$ denotes a convex functional (e.g., a norm, a semi--norm, \ldots).
This is achieved by biasing the discretization $f$ in the upwind direction.
A necessary condition for \eqref{dissipative} is monotonicity under an
explicit Euler step \cite[p. 501]{Kraaijevanger1991}:
\begin{align} \label{FEcond1}
	\| v + h f(v)\| & \le \|v\|, & \mbox{for all $v$, and for $h$ satisfying } 0 \le h \le h_0,
\end{align}
where $h_0>0$ (in general $h_0$ may depend on $v$).
Let $u_{n}, u_{n+1}$ denote approximations, computed by some numerical integrator,
to the solution at successive time steps $t_{n}$ and $t_{n+1}=t_{n}+h$.
Under the forward Euler monotonicity condition \eqref{FEcond1},
it is possible to prove
that many Runge--Kutta and linear multistep methods also give monotone
solutions; i.e., solutions that satisfy
\begin{align}
	\|u_{n+1}\| & \le \|u_{n}\| & \mbox{for $h$ satisfying } 0 \le h \le R h_0.
\end{align}
Such methods are known as  strong stability preserving (SSP) methods, and
the factor $R$ is known as the radius of absolute monotonicity or SSP coefficient
of the method.
SSP methods necessarily have non-negative coefficients, since the monotonicity
property is proved using \eqref{FEcond1} and convexity.

Monotonicity cannot be ensured using only assumption \eqref{FEcond1} in general for methods with negative coefficients \cite[Thm. 4.2]{Kraaijevanger1991}, or even for some methods
(such as the classical fourth-order RK method) with non-negative coefficients \cite[Thm. 9.6]{Kraaijevanger1991}.
In order to accommodate such methods, a second discrete approximation to $-\mathcal F$
is introduced and referred to as $\tilde{f}$.  This discretization must be monotone
under an explicit Euler step with negative step size:
\begin{align} \label{FEcond2}
	\| v - h f(v)\| & \le \|v\|, & \mbox{for all $v$, and for $h$ satisfying }
    										0 \le h \le \tilde{h}_0,
\end{align}
where $\tilde{h}_0>0$. In the context of hyperbolic problems, $\tilde{f}$ must be biased in the downwind direction,
and typically $\tilde{h}_0 = h_0$.
The downwind spatial discretization $\tilde{f}$ is to be used in place of $f$ wherever
a negative coefficient  appears in the time integration method, in order to ensure
monotonicity of the overall method.  Introduction of a downwind discretization
makes it possible to ensure monotonicity for a broader class of methods,
including the classical RK method of order four.  It also makes it possible
to ensure monotonicity for many methods under larger step sizes.

Methods that use both upwind and downwind operators can naturally be viewed as
{\em perturbed Runge--Kutta methods}. Although they are also connected to
additive RK methods (see \cite{higueras2005a,higueras2006}),
in the present work we will employ the perturbation viewpoint,
and refer to methods that use downwind discretization as perturbed RK methods.

During the last quarter century, a number of additional authors have studied
monotonicity for methods that use downwind discretization.
The main motivation for this work has been to
break the  ``order barrier'' that restricts explicit RK methods to order
four, or to find new methods with larger SSP coefficients.
In this context, numerical optimization of the SSP coefficient for RK methods with negative
coefficients was conducted for explicit methods in \cite{ruuth2004,ruuth2006,gottlieb2006}
and for implicit methods in \cite{ketcheson2011a}.
In each case, optimization was
carried out over methods with a specified order and number of stages.

The present work stems from a different motivation.  Monotonicity preservation
is not the only numerical property of interest in applications, and practitioners
may wish to use a particular integrator that has small or zero SSP coefficient.
Our goal is then to perturb the prescribed method in order to achieve
larger monotonicity-preserving timesteps.
Little work has been done in this direction, because it is not known
how to find the best perturbation for a given method.  That problem
is the main focus of this work.
Most of our results concern only explicit methods,
although one major result (Theorem \ref{thm:feasible}) pertains
also to implicit methods.

\subsection{Perturbed Runge--Kutta methods}
A Runge-Kutta method applied to the initial value problem \eqref{ode}
computes approximations $u_n \approx u(t_n)$ by
\begin{subequations} \label{Runge--Kutta}
\begin{align}
    Y & = u_n e + h K F\, ,  \\
    u_{n+1} & = Y_{s+1}.
\end{align}
\end{subequations}
Here $s$ is the number of stages, $e$ is a vector whose entries are equal to one, $Y$ is the vector containing the stage values and the numerical solution, $Y=(Y_1, \ldots, Y_s, Y_{s+1})^t$,
$[F]_i = f(Y_i)$, and $K$ is the $(s+1)\times(s+1)$ matrix of Butcher
coefficients:
$$K=\begin{pmatrix} A & 0 \\ b^t & 0 \end{pmatrix}\, .
$$
In this work we study perturbed Runge--Kutta methods:
\begin{subequations} \label{PRunge--Kutta}
\begin{align}
    Y & = u_n e + h K F + h\widetilde{\bbA}(F-\widetilde{F})\, ,  \\
    u_{n+1} & = Y_{s+1},
\end{align}
\end{subequations}
where $\widetilde{K}$ is given by
$$\widetilde{\bbA}=\begin{pmatrix} \tilde A & 0 \\ \tilde b^t & 0 \end{pmatrix}\, ,
$$
and $\widetilde{F}$ is defined analogously to $F$.
We assume that the
perturbation $\tilde A$ has the same structure
(strictly lower-triangular, lower-triangular, or full)
as the matrix $A$.

Observe that the perturbed method \eqref{PRunge--Kutta} reduces to the RK method
\eqref{Runge--Kutta} when $\widetilde{F} = F$.
Method \eqref{PRunge--Kutta} may be viewed as approximating the solution
of the perturbed problem
$$u'(t) = f(u) + (f(u) - \tilde{f}(u)),$$
where $\tilde{f}\approx f$.

We assume that $f$ and $\tilde{f}$ satisfy the explicit Euler assumptions
\eqref{FEcond1} and \eqref{FEcond2}, respectively, with $\tilde{h}_0 = h_0$.

Most previous works, including \cite{ruuth2004,ruuth2006}, have focused
on methods with the following property
\begin{quote}
	{Property C}:

    We say that a perturbation $\widetilde{K}$ to an RK method $K$ possesses
    {\em property C} if, for each value of $j$
    \begin{align}
    	\widetilde{K}_{ij}   \ne 0 \mbox{ (for some $i$)} \quad \implies  \quad
                  K_{ij}   = 0 \mbox{ (for all $i$).}
    \end{align}
\end{quote}
In words, property C means that in the $j$th column, only one of $K, \widetilde{K}$
has any nonzero entries.  Thus, only one of $f(y_j), \tilde{f}(y_j)$
need ever be evaluated, so only $s$ total function evaluations are required per step.
In \cite{gottlieb2006} it was shown that
for WENO discretizations, the cost of computing both
$f(y_j)$ and $\tilde{f}(y_j)$ is much less than twice the cost of computing $f(y_j)$
alone.  Therefore methods that without property C may also be of practical
interest.  In the present work, we do not assume property C.

\subsection{Scope and outline}
The central question of the present work is
\begin{itemize}
    \item Given a fixed RK method, what perturbation results in the largest
    		radius of absolute monotonicity for the perturbed method?
\end{itemize}
We investigate this question in the context of both linear problems
(Section \ref{sec:lin:prob}) and nonlinear problems (Section \ref{sec:nonlinear}).

For explicit methods applied to linear problems, the question above can be cast in terms of
absolute monotonicity of the (bivariate) stability polynomial.  In Section
\ref{sec:upperbound:lin} we prove a general upper bound on the radius
of absolute monotonicity of the stability polynomial of an explicit perturbed RK method
with $s$ stages and linear order $p$.
 In Section \ref{sec:upperlin:num}, we provide an algorithm for
computing tighter bounds, and tabulate some of the resulting numerical values.
Examples of optimal methods are given in \ref{sec:lin-examples}.

Section \ref{sec:nonlinear} is devoted to perturbed Runge--Kutta methods for
nonlinear problems.
Theorem \ref{thm:feasible} states that a perturbation with positive
radius of absolute monotonicity exists for every Runge--Kutta method.
Theorems \ref{prop_vpos} and \ref{Pbound} give simple upper bounds on the optimal
perturbed radius for a given explicit method.
In Section \ref{sec:algorithms} we give two algorithms for computing optimal
perturbations.  The first is provably correct but approximate, while the second is heuristic but exact and agrees with the first in all cases we have tested.  Both are applicable only to explicit methods.
These algorithms have been implemented in the software package
Nodepy \cite{nodepy}.
We conclude Section \ref{sec:nonlinear} with an application of the theorems and
algorithms to optimal perturbations of some RK methods from the literature.
Among the results is the first truly optimal perturbation for the classical
4th-order method of Kutta.

Section \ref{sec:conclusions} contains some conclusions as well as some open questions to be studied in the future.

Finally, in the Appendix we give some details on perturbations for the family of second order 2-stage methods and the classical fourth order RK method.

\section{Explicit perturbed Runge--Kutta methods for linear problems}\label{sec:lin:prob}

To study the behavior of the perturbed Runge--Kutta method \eqref{PRunge--Kutta}
for linear problems, we apply it to a linear scalar test problem,
setting $f(u)=\lambda u$ and $\tilde{f}(u)=\tilde{\lambda} u$ in \eqref{PRunge--Kutta}.
This results in the iteration $$u_{n+1} = \phi_{(\bbA, \Apert)}(z, -\tilde{z}) \, u_n,$$
where $z=h\lambda$, $\tilde{z}=h\tilde{\lambda}$ and
\begin{align}
\phi_{(\bbA, \Apert)}(z, \tilde{z})= 1 +  \left(z b^t + (z+\tilde{z}) \tilde{b}^t\right)\left(I -  z A - (z+\tilde{z}) \tilde{A} \right)^{-1} e\, . \label{bistabfunc}
\end{align}
We refer to \eqref{bistabfunc} as the {\it stability function of the perturbed
Runge--Kutta method \eqref{PRunge--Kutta}}.

We say that a function $\psi:\mathbb{R}\to \mathbb{R}$  is absolutely monotonic (a.m.) at
$\xi$ if all derivatives at $\xi$ exist and they are non-negative.
For a function $\psi:\mathbb{R}^2\to \mathbb{R}$ the concept of absolute
monotonicity can be defined in a similar way: $\psi$ is a.m.
at $(\xi, \tilde \xi)$ if all derivatives at $(\xi, \tilde \xi)$ exist and they
are non-negative.

Given a function $\psi(z, \tilde{z})$, we define the {\it radius of absolute monotonicity}  as
\begin{align}
R(\psi)=  \sup\left\{ r\in \mathbb{R} \, | \,  r=0, \hbox{ or }  r>0\, , \, \hbox{and }
\psi(z, \tilde{z}) \text{ is a.m. at } (-r, -r) \right\}\, . \label{threshold:pert}
\end{align}

For a perturbed Runge--Kutta method \eqref{PRunge--Kutta} with coefficients $(K, \widetilde{\bbA})$,
we write ${R}_{\hbox{\footnotesize Lin}}(K, \widetilde{\bbA})$ to denote the radius of absolute monotonicity of its stability
function:
$$ {R}_{\hbox{\footnotesize Lin}}(K, \widetilde{\bbA})=R(\phi_{(\bbA, \Apert)})\, .  $$
The quantity ${R}_{\hbox{\footnotesize Lin}}(K,\widetilde{\bbA})$ is referred to as the {\it threshold factor} due to its
role in the step size for monotonicity.
For a given Runge--Kutta method \eqref{Runge--Kutta} with coefficients $K$,
we are interested in determining perturbations $\widetilde{\bbA}$ that give the largest
threshold factor.  The corresponding {\it  threshold factor of the optimal perturbation} is
denoted by
\begin{align} \label{opt-thresh-fact}
  {R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)=  \sup_{\widetilde{\bbA}} {R}_{\hbox{\footnotesize Lin}}(K, \widetilde{\bbA})\, .
\end{align}
The supremum in \eqref{opt-thresh-fact} is taken over all strictly lower
triangular matrices $\widetilde{K}$, in order to preserve the explicit nature of the method.
A perturbation ${\widetilde{\bbA}}$ such that
$${R}_{\hbox{\footnotesize Lin}}(K, \widetilde{\bbA}) = {R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K) $$
will be called an {\it optimal perturbation of the method $K$ for the linear
problem}.

Taking $\widetilde{\bbA}=0$ gives a (not perturbed) Runge--Kutta method \eqref{Runge--Kutta} and a (not perturbed) stability function $\phi_K$.
In this case we denote the threshold factor ${R}_{\hbox{\footnotesize Lin}}(K, 0)$ simply by $R(\phi_K)$.
Clearly
\begin{align}R(\phi_{K}) \le {R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)\, .\label{bound:RRL}
\end{align}
In the next section, we give upper bounds on ${R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)$.

\subsection{Upper bounds on the threshold factor for optimal perturbations \label{sect:opt-dw-tf}}
For any explicit perturbed RK method $(K,\widetilde{\bbA})$ of  linear order $p$,
it can be seen that $\phi_{(K,\widetilde{\bbA})} \in \widetilde{\Pi}_{s,p}$, where
$\widetilde{\Pi}_{s,p}$, with $p\le s$, denotes the set of bivariate polynomials with
the following properties:
\begin{enumerate}
    \item $ \psi(z,\tilde{z}) = \displaystyle \sum_{j=0}^p \frac{z^j}{j!}
            + \sum_{j=p+1}^s \sigma_j z^j + (z+\tilde{z})\Psi(z,\tilde{z}) $;
    \item $\Psi$ is a polynomial of combined degree at most $s-1$.
\end{enumerate}
In this section we investigate
$$\widetilde\cR_{s,p}=\sup\left\{R(\psi) \, |\, \psi(z, \tilde{z})\in\widetilde{\Pi}_{s,p}\right\}\, ,$$
Clearly
\begin{align}{R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)\leq \widetilde\cR_{s,p}\, . \label{bound:RLRsp}
\end{align}
However, not all functions in $\widetilde{\Pi}_{s,p}$ can
be realized as the stability function of an $s$-stage perturbed Runge-Kutta
method \eqref{PRunge--Kutta}, so inequality \eqref{bound:RLRsp} is often strict (see Example \ref{ex:2eo2} below).
In case the optimal polynomial is
realizable, the corresponding method may be
of interest for the integration of linear systems.

In Subsection \ref{sec:upperbound:lin} we give an upper bound for
$\widetilde\cR_{s,p}$.  In Subsection \ref{sec:upperlin:num}, we give an algorithm to
compute, $\widetilde\cR_{s,p}$ for given $s$ and $p$, along with numerical values.

\subsubsection{Upper bound on $\boldsymbol{\widetilde\cR_{s,p}}$}
\label{sec:upperbound:lin}~\\

\begin{lemma} \label{lem:1sp}
Let $\varphi(z)$ be a polynomial satisfying
\begin{align}\label{polpert}
\varphi(z) & = 1 + \gamma_1 z +\cdots + \gamma_p z^p +\gamma_{p+1} z^{p+1} + \cdots + \gamma_s z^s\\
\gamma_j & \geq \frac{1}{j!}, \quad  j=1, \ldots, p\,  . \nonumber
\end{align}
Then the radius of absolute monotonicity of  $\varphi$ satisfies
\begin{align}R(\varphi) \le\sqrt[p]{s(s-1)\ldots (s-p+1)}\, . \label{R1sp}
\end{align}
\end{lemma}

\begin{proof}  If  $R(\varphi)=0$, inequality \eqref{R1sp}
is trivial.  Let $\varphi(z)$ satisfy \eqref{polpert}
and be absolutely monotonic at $-r$ with $r>0$.
Then it can be written as
$$\varphi(z)=\sum_{j=0}^s \alpha_j \left(1 + \frac{z}{r}\right)^j=\sum_{j=0}^s \alpha_j \left(\sum_{\ell=0}^j  \frac{z^\ell}{r^\ell} {j \choose \ell}\right)=\sum_{\ell=0}^s \left(\sum_{j=\ell}^s \alpha_j  {j \choose \ell}\right ) \frac{z^\ell}{r^\ell}\,,
$$
where $\alpha_j\geq 0$, and $\sum_{j} \alpha_j=1$. As $\varphi$ is of the form \eqref{polpert},   the coefficient of $z^p$ is larger than $1/p!$. Some computations give
$$\frac{1}{p !} \leq \left(\sum_{j=p}^s \alpha_j  {j \choose p}\right ) \frac{1}{r^p} \leq  \left(\sum_{j=p}^s \alpha_j  \right ) {s \choose p} \, \frac{1}{r^p} \leq  {s \choose p} \frac{1}{r^p}=\frac{s \, (s-1) \cdots (s-p+1)}{ p ! \, \, r^p}.
$$
Consequently,
$$r\leq \sqrt[p]{s(s-1)\cdots (s-p+1)}\, .  $$
\end{proof}

We remark that equality in \eqref{R1sp} is obtained for the polynomial
\begin{equation}\varphi(z)=\left(1+ \frac{z}{r}\right)^s\, , \label{optTihi}
\end{equation}
where $r=\sqrt[p]{s(s-1)\cdots (s-p+1)}$.

With the previous results, we can prove an upper bound on $\widetilde\cR_{s,p}$.

\begin{theorem}
\begin{align}\widetilde\cR_{s,p}\leq \sqrt[p]{s(s-1)\cdots (s-p+1)}\, .  \label{incRtsp}
\end{align}
\end{theorem}
\begin{proof}
If $R(\psi)=0$ for all $\psi \in \widetilde{\Pi}_{s,p}$, then $\widetilde\cR_{s,p}=0$ and
inequality \eqref{incRtsp} is true. Otherwise, there exists a function  $\psi
\in \widetilde{\Pi}_{s,p}$  a.m. at $(-r,-r)$ with $r>0$.
By \cite[Lemmas 2.9 and 2.10]{higueras2006}, $\psi$ is a.m.
at the points $(\xi, \xi)$, with $\xi \in [-r, 0]$.
Writing $\psi(z, \tilde{z}) = \sum \sum \mu_{jk} z^j \tilde{z}^k$
and differentiating shows that all coefficients $\mu_{jk}$ are non-negative
since $\psi$ is a.m. at $(0,0)$.  Thus  $\psi(z,z)$ (viewed as a function
of one variable) is of the form \eqref{polpert} and is a.m. at $-r$.
Application of Lemma \ref{lem:1sp} gives the desired result.
\end{proof}

\subsubsection{Numerical computation of bounds $\boldsymbol{\widetilde\cR_{s,p}}$}\label{sec:upperlin:num}
In this section we provide a means to compute tighter bounds on $\widetilde\cR_{s,p}$ for
given $s,p$, using linear programming.
The material in this section closely follows \cite[Section 4.6.2]{ketchesonphdthesis}.
The stability function \eqref{bistabfunc} of an explicit $s$-stage perturbed Runge--Kutta method
\eqref{PRunge--Kutta} with linear order $p$ can also be written in the form
\begin{align} \label{dwtaylor1}
\psi(z,\tilde{z})= \sum_{j=0}^s \sum_{\ell=0}^j \gamma_{j\ell}
  \left(1+\frac{z}{r}\right)^{j-\ell} \left(1+\frac{\tilde{z}}{r}\right)^{\ell}
  \quad \quad \mbox{with } \gamma_{j\ell} = \frac{r^j}{j!}\frac{\partial^j\psi_i}{\partial z^{j-\ell} \partial \tilde{z}^{\ell}};
\end{align}
furthermore, $\psi(z,-z) =\phi_K(z)= \exp(z)+ {\mathcal O}(z^{p+1})$.
After considerable manipulation we find that
$\psi(z,-z)=\sum_{i=0}^s C_i z^i$ where
\begin{align*}
C_i(r,\gamma) = \sum_{j=i}^s \sum_{\ell=0}^j \gamma_{j\ell} \sum_{m=\max(0,i-\ell)}^{\min(i,j-\ell)}
{j-\ell \choose m}{\ell \choose i-m} \frac{(-1)^{i-m}}{r^i}
\end{align*}

Hence we have the following problem for existence of a polynomial
\eqref{dwtaylor1} with perturbed threshold factor at least $r$ and order at least $p$:
\begin{subequations}
\label{optprob3}
\begin{align}
\mbox{Given } r \mbox{ find $\gamma$ such that} & \nonumber \\
& \gamma_{j\ell} \ge 0 & 0 \le \ell\le j \le s \\
& C_i(r,\gamma) = \frac{1}{i!} & 0 \le i \le p.
\label{eqconstr}
\end{align}
\end{subequations}
Since \eqref{eqconstr} is a system of linear equations (in $\gamma$)
then for any given value of $r$ \eqref{optprob3} represents a linear programming
feasibility problem.  Hence we can use bisection and an LP solver to find
$\widetilde\cR_{s,p}$, as was done for similar problems in \cite{ketcheson2008,ketcheson2009a}.
Table \ref{dwtable} gives computed values of $ \widetilde\cR_{s,p}$ for $s$ and $p$ up to ten.

\begin{table}
\begin{center}
\caption{$\widetilde\cR_{s,p}$: upper bounds on the threshold factors for optimal perturbations\label{dwtable}}
\begin{tabular}{|l|rrrrrrrrrr|}
\hline
s p  & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 \\ \hline
1  & 1.00  & &  &  &  &  &  &  &  &  \\
2  & 2.00  & 1.41 &  &  &  &  &  &  &  &  \\
3  & 3.00  & 2.45 & 1.60 &  &  &  &  &  &  &  \\
4  & 4.00  & 3.46 & 2.49 & 2.00 &  &  &  &  &  &  \\
5  & 5.00  & 4.47 & 3.20 & 2.94 & 2.18 &  &  &  &  &  \\
6  & 6.00  & 5.48 & 4.00 & 3.65 & 3.11 & 2.58 &  &  &  &  \\
7  & 7.00  & 6.48 & 4.86 & 4.45 & 3.88 & 3.55 & 2.76 &  &  &  \\
8  & 8.00  & 7.48 & 5.77 & 5.31 & 4.57 & 4.32 & 3.72 & 3.15 &  &  \\
9  & 9.00  & 8.49 & 6.62 & 6.22 & 5.24 & 5.02 & 4.52 & 4.14 & 3.33 &  \\
10 & 10.00 & 9.49 & 7.42 & 7.09 & 5.95 & 5.70 & 5.25 & 4.96 & 4.32 & 3.73 \\
\hline
\end{tabular}
\end{center}
\end{table}

\subsection{Examples}\label{sec:lin-examples}
We now give some examples of optimal polynomials and Runge--Kutta methods.

\subsubsection{Polynomials achieving $\boldsymbol{\widetilde\cR_{s,p}}$}
The algorithm just described also provides coefficients for an optimal
polynomial, which may or may not be realizable as the stability function of
a perturbed Runge--Kutta method.

The optimal first-order polynomial for any $s$ is just the stability polynomial
of a (not perturbed) Runge--Kutta method consisting of $s$ repeated forward Euler steps.

The optimal order-two polynomial of degree $s$ also has a simple form:
\begin{align} \label{optpoly2}
\psi_{s,2} (z, \tilde{z})& = \frac{2(s+r)-1}{2(s+r)}\left(1+\frac{z}{r}\right)^s +
\frac{1}{2(s+r)}\left(1+\frac{\tilde{z}}{r}\right)^s,
\end{align}
where $r=\widetilde\cR_{s,2}=\sqrt{s(s-1)}$. Observe that bound \eqref{incRtsp} is sharp for $p=2$.

Some of the other optimal polynomials also have rational coefficients.
Two optimal degree-four fourth order polynomials we found are
\begin{align*}
\psi_{4,4}^{1} (z, \tilde{z})& = \frac{1}{3}\left(1+\frac{z}{r}\right)^2 +
\frac{17}{48}\left(1+\frac{z}{r}\right)^4 +
\frac{14}{48}\left(1+\frac{z}{r}\right)^2\left(1+\frac{\tilde{z}}{r}\right)^2 +
\frac{1}{48}\left(1+\frac{\tilde{z}}{r}\right)^4\, ,
\end{align*}
and
$$ \psi_{4,4}^{2}(z, \tilde{z})= \frac{7}{16} \left(1+\frac{z}{r}\right)^4+\frac{3}{8}  \left(1+\frac{z}{r}\right)^2 \left(1+\frac{\tilde{z}}{r}\right)^2+\frac{1}{6}
   \left(1+\frac{z}{r}\right)^3 \left(1+\frac{\tilde{z}}{r}\right)+\frac{1}{48} \left(1+\frac{\tilde{z}}{r}\right)^4\, , $$
where $r=\widetilde\cR_{4,4}=2$.
Thus the optimal polynomial in $\widetilde{\Pi}_{s,p}$ is in general not unique.

We have also computed (not shown in this paper) the exact values of
$\widetilde\cR_{s,p}$ and polynomials $\psi_{s,p}$ for $p=3$ and $s=3,4,5,6$. In each
case,  we found an optimal polynomial of the form \eqref{dwtaylor1}
where $p$ of the coefficients $\gamma_{sj}$ are non-zero and the coefficients
$\gamma_{ij}$ with $i<s$ are all zero.

Exact values of  $\widetilde\cR_{s,p}$ can be found in a systematic way as follows.
First, a high-precision approximation can be computed using bisection and
linear programming as described above.
In practice, this yields a set of coefficients $\gamma_{jl}$ in which only $p$ values are non-zero  \cite{hadjimichael_draft}.
Setting
the remaining values to zero {\em a priori} in \eqref{eqconstr} yields a
system of $p+1$ equations in $p$ unknowns which, nevertheless, possesses
a solution.  The solution may be found using a symbolic linear algebra package.
We do not pursue this further in the present work.

\begin{remark}
As noted already, not all polynomials of the form \eqref{dwtaylor1} can be realized as
the stability function of a perturbed Runge--Kutta method \eqref{PRunge--Kutta} with $s$ stages.
For example, the polyomial \eqref{optpoly2} with $s=2$ is not the stability
function of any two-stage method (i.e. using only evaluations of
$f(u_n), \tilde{f}(u_n), f(y_1), \tilde{f}(y_1)$).
It can be realized as the stability function of a method that has three stages,
using evaluations of $f(u_n), \tilde{f}(u_n), f(y_1), \tilde{f}(y_2)$.  The difference
in cost between such methods depends on the nature of $f, \tilde{f}$; see \cite{gottlieb2006}.
For this reason, we stress that the values
in Table \ref{dwtable} are only {\em upper bounds} on what can be achieved.
We do not pursue the topic further here.   \hfill $\square$
\end{remark}

\subsubsection{Optimal threshold factors for perturbations of specified 2-stage 2nd-order methods}
We have no general method for finding ${R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)$ nor a corresponding
method.  In this section we report results of some symbolic searches.
In the case of the second-order methods, due to the small number of free
parameters, it is not difficult to prove that the results below are truly
optimal.

\begin{example}\label{ex:2eo2}
We consider explicit perturbed second-order 2-stage methods
\begin{align}\begin{tabular}{c|cc}
$0$ & $0$ & $0$ \\
 $\alpha$ & $\alpha$ & $0$ \\
\hline
$K$ \phantom{\huge I}  & $1-\frac{1}{2 \alpha}$ & $\frac{1}{2 \alpha}$
\end{tabular}
\qquad \begin{tabular}{c|cc}
  & $0$ & $0$ \\
  & $\tilde a_{21}$ & $0$ \\
\hline
$\widetilde{\bbA}$ \phantom{\huge I}  & $\tilde{b}_1$ & $\tilde{b}_2$
\end{tabular}\label{2eo2OptLin}.
\end{align}
For these methods, function \eqref{bistabfunc} can be expanded as
\begin{align}\phi_{(\bbA, \Apert)}(z,\tilde{z})= 1 + z + \frac{1}{2} z^2 + \beta_{11} z (z+\tilde{z}) + \beta_{1} (z+\tilde{z})+ \beta_2 (z+\tilde{z})^2\, , \label{2estord2}
\end{align}
where
\begin{align}\label{betas}
\beta_{11}= b^t \tilde{A} e + \tilde{b}^t A e= \tilde{b}_2 a_{21} + b_2 \tilde a_{21} \, , \qquad  \beta_1= \tilde{b}^t e= \tilde{b}_1+ \tilde{b}_2 \, , \qquad \beta_2=\tilde{b}^t \tilde{A} e=\tilde{b}_2 \tilde a_{21}
\end{align}
The polynomial \eqref{2estord2} is realizable (in the sense that it corresponds to a 2-stage Runge--Kutta method \eqref{2eo2OptLin}) if the first and last equations in \eqref{betas} can be solved for $\tilde a_{21}$ and $\tilde{b}_2$ in $\mathbb{R}$. A simple computation gives that a necessary condition is
$\beta_{11}^2-2 \beta_2\geq 0$.  Note that the stability function is independent of $\alpha$.

With the help of Maple, we have computed the largest $r$ such that \eqref{2estord2} is a.m. at $(-r, -r)$ and the polynomial is realizable.  We have obtained that the optimal perturbation, denoted by $\widetilde{\bbA}_L$, satisfies $\tilde{b}_2=\tilde a_{21}=0$ and $\tilde{b}_1=\frac{1}{3} \left(\sqrt{7}-2\right)$, and \begin{equation}{R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)=\frac{1}{3} \left(1+\sqrt{7}\right)\approx 1.21525\, . \label{RLinopt}
\end{equation}
Observe that ${R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)< \widetilde\cR_{2,2}=\sqrt{2}$.
The stability function \eqref{bistabfunc} for the  optimal perturbed method is
$$\tilde{\phi}_{(K, \widetilde{\bbA}_L)}(z,\tilde{z})= \frac{1}{9} \left(4+\sqrt{7}\right) \left(1+ \frac{z}{r}\right)^2+ \frac{1}{9} \left(5-\sqrt{7}\right) \left(1+ \frac{\tilde{z}}{r}\right) \, , \label{2estord2opt}
$$
where $r={R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)$.
\end{example}

\begin{example} \label{ex:4eo4}
 We consider perturbations of the classical four stage, order four method, of the form
\begin{equation}
\begin{tabular}{c|cccc}
$0$ & $0$ & $0$ & $0$ & $0$ \\[1ex]
$\frac{1}{2}$ & $\frac{1}{2}$ & $0$ & $0$ & $0$ \\[1ex]
$\frac{1}{2}$ & $0$ &  $\frac{1}{2}$ &$0$ & $0$ \\[1ex]
$1$ & $0$ & $0$ & $1$ & $0$ \\[1ex]
\hline
$K$ \phantom{\huge I} & $\frac{1}{6}$  & $\frac{1}{3}$& $\frac{1}{3}$ & $\frac{1}{6}$
\end{tabular} \qquad \begin{tabular}{c|cccc}
 & $0$ & $0$ & $0$ & $0$ \\[1ex]
 & $0$ & $0$ & $0$ & $0$ \\[1ex]
 & $\tilde a_{31}$ & $0$ & $0$ & $0$ \\[1ex]
& $\tilde a_{41}$ & $\tilde a_{42}$ & $0$ & $0$ \\[1ex]
\hline
$\widetilde{\bbA}$ \phantom{\huge I}  & $\tilde{b}_{1}$& $\tilde{b}_{2}$ & $0$ & $0$
\end{tabular}
\label{4eo4}
\end{equation}
We consider these perturbations because,
 in order to obtain a nonzero  SSP coefficient for nonlinear problems, the analysis done in \cite{higueras2005a} shows that only the entries $\tilde a_{31}$, $\tilde a_{41}$, $\tilde a_{42}$, $\tilde{b}_1$ and $\tilde{b}_2$ in $\widetilde{\bbA}$ need be nonzero. To study SSP coefficients for the linear case, we have to analyze the perturbed stability function, that in this case  is of the form
\begin{align}\phi_{(\bbA, \Apert)}(z,\tilde{z}) = 1 + z + \frac{1}{2} z^2 +   \frac{1}{6} z^3+  \frac{1}{24} z^4+   \beta_{1} (z+\tilde{z})   + \beta_{11} z (z+\tilde{z})+ \beta_{21} z^2 (z+\tilde{z})\,  \label{esfun:4eo4}
\end{align}
where
$$\beta_{1}= \tilde{b}_1+\tilde{b}_2\, , \qquad \beta_{11}= \frac{1}{6} \left(3 \,\tilde{b}_2+2 \,\tilde a_{31}+\tilde a_{41}+\tilde a_{42}\right)\, , \qquad \beta_{21}=  \frac{1}{12} \left(2 \,\tilde a_{31}+\tilde a_{42}\right)  \, .
$$
After some computations with \eqref{esfun:4eo4}, we obtain a coefficient
${R}_{\hbox{\footnotesize Lin}}(K,\widetilde{\bbA})\approx 1.66728$ that is the positive root of  the polynomial $15 \, x^4-4 \,x^3-12 \,x^2-24 \,x-24=0$.  The coefficients are
$$\beta_{1}=\frac{ 7 \,r^3-2\, r^2-6 \,r-12 }{12} \, , \qquad  \beta_{11}=\frac{ 5\, r^2-2 \,r-6 }{12} \, , \qquad   \beta_{21}=\frac{r-1}{6}\, .
$$
where $r={R}_{\hbox{\footnotesize Lin}}(K,\widetilde{\bbA})$. With these values, the perturbed stability function can be written as
$$\phi_{(\bbA, \Apert)}(z,\tilde{z}) =
  \gamma_{01} \left(1+\frac{\tilde{z}}{r}\right) +\gamma_{11}
   \left(1+\frac{z}{r}\right)\left(1+\frac{\tilde{z}}{r}\right)+
  \gamma_{21}
  \left(1+\frac{z}{r}\right)^2 \left(1+\frac{\tilde{z}}{r}\right) + \gamma_{40} \left(1+\frac{z}{r}\right)^4\, ,
$$
where
$$\gamma_{01}=\frac{r \left(2\, r^3-r^2-6\right)}{6} \, ,    \quad \gamma_{11}=\frac{r^2 \left(r^2+2 \,r-6\right)}{12} \, , \quad \gamma_{21}=\frac{r^3\, (r-1) }{6}  \, ,  \quad \gamma_{40}= \frac{r^4}{24} \, . $$
 This perturbed stability function can be realized with the family of perturbations
\begin{subequations}\label{coef4eo4L}
\begin{align*}
\tilde a_{31}&= \frac{1}{2} \left(2 \, r-2-\tilde a_{42}\right)\, ,\, \\ \tilde a_{41}&=\frac{1}{2} \left(5  \,r^2-6 \,
   r-2-6  \,\tilde{b}_2\right)\,, \\ \tilde{b}_1&= \frac{1}{12} \left(7  \,r^3-2  \,r^2-6  \,r-12-12  \,\tilde{b}_2\right).
\end{align*}
\end{subequations}
\end{example}

\section{Perturbed Runge--Kutta methods for nonlinear problems}\label{sec:nonlinear}

In this section we seek to answer the question posed in the
introduction: for a given method $K$, what perturbation $\widetilde{K}$ gives
the largest value of $R(K,\widetilde{K})$?  We begin by providing
some upper bounds.

It is convenient to write
scheme \eqref{Runge--Kutta} in canonical Shu-Osher form \cite{SSPbook}
\begin{align} \label{cso}
    Y & = v_r u_n + \alpha_r \left( Y + \frac{h}{r} F\right)
\end{align}
where
\begin{align} \label{canonical-coefficients}
   v_r       = (I + rK)^{-1} e\, , \qquad
     \alpha_r   = r(I + rK)^{-1} K\, .
\end{align}
Observe that matrices $K$ and $\alpha_r$ have the same structure (strictly lower triangular, lower triangular or full).

A Runge--Kutta method \eqref{Runge--Kutta} is said to be absolutely monotonic at $r$ if $(I + rK)^{-1}$
exists and all entries of $\alpha_r, v_r$ are non-negative. For a given method
$K$, the largest $r$ such that the method is absolutely monotonic at $r$ is
known as the SSP coefficient, Kraaijevanger coefficient,
or radius of absolute
monotonicity \cite{Kraaijevanger1991}, and it is denoted herein by $R(K)$:
\begin{align}R(K)=\sup\left\{r\, | \, r=0 \text{ or } r>0, (I + rK )^{-1} \hbox{ exists, and }   \alpha_r,v_r\ge 0 \right\}\, . \label{ramdef0}
\end{align}
As usual, the inequalities above should be understood component--wise.

Next, we consider perturbed Runge--Kutta methods \eqref{PRunge--Kutta}.
To study absolute monotonicity of perturbed Runge--Kutta
methods, we write method \eqref{PRunge--Kutta} also in a canonical
Shu-Osher-like form
\begin{align} \label{dwRunge--Kutta_cso}
Y & = \gamma_r u_n + \alpha^\textup{up}_r  \left(Y + \frac{h}{r}F\right) +
                    \alpha^\textup{down}_r \left(Y - \frac{h}{r}\widetilde{F}\right)\, ,
\end{align}
where
\begin{subequations} \label{pert-cso-coeff}
\begin{align}
\gamma_r & = (I + rK + 2r\widetilde{\bbA})^{-1} e\,  ,\label{pert-cso-coeffa}\\[0.5ex]
\alpha^\textup{up}_r & = r(I + rK + 2r\widetilde{\bbA})^{-1} (K + \widetilde{\bbA})\, ,  \\[0.5ex]
\alpha^\textup{down}_r & = r(I + rK + 2r\widetilde{\bbA})^{-1} \widetilde{\bbA}\, .
\end{align}
\end{subequations}
Observe that method  \eqref{dwRunge--Kutta_cso}, with $ \gamma_r=(I-\alpha^\textup{up}_r-\alpha^\textup{down}_r)e$, is a perturbed Runge--Kutta scheme with Butcher coefficients
\begin{align}K=\frac{1}{r}(I-\alpha^\textup{up}_r-\alpha^\textup{down}_r)^{-1} (\alpha^\textup{up}_r-\alpha^\textup{down}_r)\, , \qquad \widetilde{\bbA}=\frac{1}{r}(I-\alpha^\textup{up}_r-\alpha^\textup{down}_r)^{-1} \alpha^\textup{down}_r\, , \label{AApert} \end{align}
provided that $(I-\alpha^\textup{up}_r-\alpha^\textup{down}_r)^{-1}$ exists.

The {\em radius of absolute monotonicity} of a perturbed Runge--Kutta method  $(K,\widetilde{\bbA})$ is
the largest $r$ such that $\gamma_r$, $\alpha^\textup{up}_r$ and $\alpha^\textup{down}_r$ in \eqref{pert-cso-coeff} exist and are non-negative
\cite[Definition 3.1]{higueras2005a}:
\begin{align}R(K,\widetilde{\bbA})=\sup\left\{r\, | \, r=0 \text{ or } \, r>0, \,(I + rK + 2r\widetilde{\bbA})^{-1} \hbox{ exists, and }   \gamma_r,\, \alpha^\textup{up}_r,\, \alpha^\textup{down}_r\ge 0 \right\}\, . \label{ramdef}
\end{align}

\subsection{Zero-well-defined perturbations}
Regularity of $(I-\alpha^\textup{up}_r-\alpha^\textup{down}_r)$ is evidently important
in our study. Observe that from \eqref{pert-cso-coeff} we have
\begin{equation}\label{zerowell}
(I-\alpha^\textup{up}_r-\alpha^\textup{down}_r)(I+r K)=(I-2 \alpha^\textup{down}_r)\, .
\end{equation}
Consequently, if $I+r K$ is regular for some $r$, then $(I-\alpha^\textup{up}_r-\alpha^\textup{down}_r)$
is regular if and only if $(I-2 \alpha^\textup{down}_r)$ is regular.

If $I-2\alpha^\textup{down}_r$ is singular, then the stage equations
do not have a unique solution even for the trivial ODE given
by $f=0$.  Hence we say that methods for which
$I-2\alpha^\textup{down}_r$ is singular are not {\em zero-well-defined}.  See \cite[Chap. 3]{SSPbook}
for the analogous definition in the context of traditional Runge--Kutta methods.

\subsection{Optimal perturbations}
The {\em optimal perturbed SSP coefficient} of a Runge--Kutta method $K$ is
denoted by
$$R^\textup{opt}(K)=\sup_{\widetilde{\bbA}} R(K,\widetilde{\bbA})\, .
$$
For a given method $K$ that is (explicit/diagonally implicit/fully implicit),
we consider the supremum  over perturbations $\widetilde{\bbA}$ that are zero-well-defined and
correspond to the same class of methods.
A matrix $\widetilde{\bbA}$ such that $R(K,\widetilde{\bbA})=R^\textup{opt}(K)$ is called an
optimal perturbation. For some methods, the optimal perturbation is not unique.

The following result shows that every method can be perturbed so as to
give a method with strictly positive SSP coefficient.
\begin{theorem}\label{thm:feasible}
Let $K$ be a Runge-Kutta method that belongs to a specified class of methods
(explicit, diagonally implicit, or fully implicit).  Then it is always possible to find a perturbation $\widetilde{\bbA}$ within the same class such that $R(K,\widetilde{\bbA})>0$.
\end{theorem}
\begin{proof}From \cite[Proposition 3.7]{higueras2005a}, we have
$R(K,\widetilde{\bbA})>0$ if and only if the Butcher coefficients satisfy
\begin{align}\label{ineqa}
    K + \widetilde{\bbA} \geq 0 \, , \qquad \widetilde{\bbA} \geq 0\, ,
\end{align}
and the following inequalities  hold,
\begin{subequations} \label{ineqb}
\begin{eqnarray}
&& \hbox{Inc }((K+ 2 \widetilde{\bbA})  (K+ \widetilde{\bbA})) \leq
\hbox{Inc }(K+ \widetilde{\bbA})\, , \label{3.6} \\[1ex]
&& \hbox{Inc }((K+ 2 \widetilde{\bbA}) \widetilde{\bbA} ) \leq
\hbox{Inc }(\widetilde{\bbA})\, ,
\end{eqnarray}
\end{subequations}
where $\hbox{Inc } (F)$ denotes the incidence matrix of matrix $F$ defined as
$\hbox{Inc } (F)=(g_{ij})$ where $g_{ij}=1$ if $f_{ij}\neq 0$, and $g_{ij}=0$
if $f_{ij}=0$.

Consider first the implicit case.  By making all entries of $\widetilde{\bbA}$ positive
we can satisfy \eqref{ineqb}, and by making them large enough we can satisfy
\eqref{ineqa}.  For the explicit and diagonally implicit cases, note that if $K, \widetilde{\bbA}$ are (strictly) lower-triangular, then the left-hand sides
of \eqref{ineqb} are also.  Thus by making all the (strictly) lower-triangular
entries of $\widetilde{\bbA}$ positive, and by taking them large enough, we can
satisfy the above inequalities.
\end{proof}

Observe that a perturbed Runge--Kutta method $(K,\widetilde{\bbA})$ can be interpreted as an
additive Runge--Kutta method $(K+\widetilde{\bbA},\widetilde{\bbA})$ for functions $(f, \tilde{f})$,
and conditions \eqref{pert-cso-coeff} are the ones required for the absolute
monotonicity of this additive scheme at $(z, \tilde{z})=(-r, -r)$ (see
\cite{higueras2006}). From Lemma 2.8 in \cite{higueras2006}, we obtain that
the stability function $\phi_{(\bbA, \Apert)}$ defined by \eqref{bistabfunc},
is  absolutely monotonic at $(\xi, \tilde \xi)=(-r, -r)$. Consequently,
\begin{equation} R(K,\widetilde{\bbA})\leq {R}_{\hbox{\footnotesize Lin}}(K,\widetilde{\bbA})\leq {R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K) \, . \label{bound1}
\end{equation}
Furthermore, from  SSP theory and inequality \eqref{bound:RRL},  we have
\begin{equation}R(K)\leq R(\phi_K)\leq {R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)\, , \qquad  R(K)\leq  R^\textup{opt}(K)\leq {R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)\, . \label{bound2}
\end{equation}
The following example illustrates that $R(\phi_K)$ can be either larger
or smaller than $R^\textup{opt}(K)$.

\begin{example}\label{ex_2stage2order}
We consider the family of second order 2-stage Runge-Kutta methods \eqref{2eo2OptLin} for $\alpha\in\mathbb{R}$. For this family we have
\begin{align*}
v_r\geq 0 \quad & \Longleftrightarrow \qquad \alpha>0 \qquad \hbox{and}\quad 0\leq r\leq\frac{1}{\alpha}\, , \\[0.5ex]
\alpha_r\geq 0 \qquad & \Longleftrightarrow \qquad \alpha\geq\frac{1}{2} \quad \hbox{and}\quad 0\leq r\leq\frac{2\alpha-1}{\alpha}\, .
\end{align*}
Thus
\begin{equation}\label{R2stageorder2}
R(K)=\begin{cases} 0, & \hbox{if} \quad \displaystyle  \alpha\leq \frac{1}{2}\, , \\[1ex]
\displaystyle  \frac{2 \, \alpha -1}{\alpha }\, ,  & \displaystyle \hbox{if}   \quad\frac{1}{2}<\alpha \leq 1\, ,\\[1ex]
\displaystyle   \frac{1}{\alpha }\, ,  & \hbox{if}   \quad 1<\alpha\, .
\end{cases}
\end{equation}
In Figure \ref{f22} we show the threshold factor $R(\phi_K)$ (thin solid blue line) and the SSP coefficient  $R(K)$ (thick solid black line). We also show the corresponding optimal coefficients for perturbed methods, namely,  the optimal threshold factor ${R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)$  (thin dashed blue line) and the optimal SSP coefficient $R^\textup{opt}(K)$ for the perturbed method (thick dashed black line).

We see that for optimal SSP method ($\alpha=1$) it is not possible to increase the SSP coefficient by means of perturbations. However, for $\alpha=(\sqrt{7}-1)/2$ it is possible to obtain a perturbation that raises the SSP coefficient to $R^\textup{opt}(K)={R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)=\left(1+\sqrt{7}\right)\approx 1.21525$ (see \eqref{RLinopt}).

 We also see that for  $2/3< \alpha<1$ we obtain
$R(\phi_K) <R^\textup{opt}(K)
$, whereas for $0<\alpha < 2/3$  and for $1< \alpha$
we have $R^\textup{opt}(K)  < R(\phi_K) $.

Coefficients of the perturbations that give rise to these values are given in
Appendix \ref{sec:2stage-details}.
\end{example}

\begin{figure}\begin{center}\includegraphics[scale=0.45]{Figures/2stage2order}
\caption{Family of second order 2-stage methods: SSP coefficients for unperturbed methods and optimal SSP coefficients for perturbed methods.
}
\label{f22}
\end{center}
\end{figure}

\subsection{Upper bounds on the SSP coefficient for perturbed RK methods\label{sec:nonlinear-bounds}}
In this section, we explore some upper bounds on the SSP coefficient $R^\textup{opt}(K)$ where $K$ is an $s$-stage order $p$ method.
A straightforward upper bound is obtained from inequality \eqref{bound2} and Theorem \ref{incRtsp}:
\begin{align} R^\textup{opt}(K)\leq \sqrt[p]{s(s-1)\ldots (s-p+1)}\,  .\label{bound4}
\end{align}

As the next Theorem shows, the largest positive value such that vector $v_r$ in \eqref{canonical-coefficients} is non-negative is also an upper bound for $R^\textup{opt}(K)$.
\begin{theorem}\label{prop_vpos}
Consider an explicit Runge-Kutta method $K$  and let $r_e$ be the largest positive value such that vector $v_r$ in \eqref{canonical-coefficients} is non-negative. Then
\begin{equation}
R^\textup{opt}(K)\leq r_e\, .
\end{equation}
\end{theorem}
\begin{proof}
Let $r=R^\textup{opt}(K)$.  Then $\gamma_r=(I-\alpha^\textup{up}-\alpha^\textup{down})e\geq 0$, and thus
from \eqref{canonical-coefficients}  and \eqref{zerowell}   we get
\begin{equation}(I-2\alpha^\textup{down}_r)v_r   \geq 0\, . \label{vpos}
\end{equation}
As $\alpha^\textup{down}_r\geq 0$, and since we consider only explicit, zero-well-defined perturbations,
$I-2\alpha^\textup{down}_r$ is an $M$ matrix.  Thus $(I-2\alpha^\textup{down}_r)^{-1}\geq 0$. If we multiply \eqref{vpos} by  $(I-2\alpha^\textup{down}_r)^{-1}$ we obtain that $v_r\geq 0$.
\end{proof}

From Theorem \ref{prop_vpos} we obtain that
\begin{equation}
{R}(K)\leq R^\textup{opt}(K)\leq r_e\, . \label{in_re}
\end{equation}
Consequently, for those methods such that ${R}(K)=r_e$, the SSP coefficient cannot be increased by perturbation. This is the case for the family of second-order two-stage methods.
For $\alpha\geq 1$, ${R}(K)= r_e=1/\alpha$ (see Example \ref{ex_2stage2order}).

On the other hand, if ${R}(K)<r_e$ one can try to find a perturbation to increase the SSP coefficient. This is the case for the classical 4-stage order 4 method for which ${R}(K)=0$ and  $r_e\approx 1.2956$, the real root of $x^3 - 2 x^2 + 4 x - 4=0$.

Another interesting bound, for explicit methods only, can be obtained in terms of the
Butcher coefficients of the Runge-Kutta method $K$.
\begin{theorem} \label{Pbound} Consider an explicit Runge-Kutta method $(K, \widetilde{\bbA})$ with perturbed SSP coefficient $R^\textup{opt}(K)>0$. Let $K=(a_{ij})$.  Then
\begin{align}R^\textup{opt}(K)\leq  \frac{1}{\max_{ij} |a_{ij}|}\label{Rpert}
\end{align}
\end{theorem}
\begin{proof} The proof is similar to that of \cite[Lemma 3.2]{ruuth2004}.
Consider an optimal perturbation $\widetilde{\bbA}$ and set $r=R^\textup{opt} (K, \widetilde{\bbA})>0$;  consider too the canonical representation
\eqref{dwRunge--Kutta_cso}. Let
$\Lambda=\alpha^\textup{up}_r+\alpha^\textup{down}_r=(\alpha_{ij})$, $\Gamma=
\alpha^\textup{up}_r/r=(\beta_{ij})$, $\tilde \Gamma= \alpha^\textup{down}_r/r=(\tilde\beta_{ij})$;
observe that $\Lambda, \Gamma, \tilde \Gamma\geq 0$, and that $\Lambda=r
(\Gamma+\tilde \Gamma)$.  As $(I-\Lambda)e=\gamma_r\geq 0$ and $\alpha_{ik}\geq
0$, we have $\alpha_{ik}\leq 1$; as $(I-\Lambda)K=\Gamma-\tilde \Gamma$, we
have
\begin{align}a_{ik}= \beta_{ik}-\tilde \beta_{ik} + \sum_{j=k+1}^{i-1} \alpha_{ij} a_{jk}\, . \label{Acoef}
\end{align}
As $\alpha_{ik}= r (\beta_{ik} + \tilde \beta_{ik})$, then $\beta_{ik} + \tilde \beta_{ik}= \alpha_{ik}/r\leq 1/r$. In particular, from \eqref{Acoef},
$$|a_{21}|=\left|\beta_{21} - \tilde \beta_{21}\right|\leq \beta_{21} + \tilde \beta_{21}\leq \frac{1}{r}\, . $$
We proceed by induction on row $\ell$ of $K$. Assume that $|a_{ij}|\leq 1/r$, for $i= 2, \ldots, \ell$,  $j=1, \ldots , \ell-1$, and consider row $\ell+1$. Then, from  \eqref{Acoef},
\begin{align*}|a_{\ell+1, 1}|&=\left|\beta_{\ell+1, 1}-\tilde \beta_{\ell+1, 1} + \sum_{j=2}^{\ell} \alpha_{\ell+1, j} \, a_{j, 1}\right|\leq \beta_{\ell+1, 1}+\tilde \beta_{\ell+1, 1} + \sum_{j=2}^{\ell} \alpha_{\ell+1, j} |\, a_{j, 1}|\\
&\leq  \frac{1}{r} \alpha_{\ell+1, 1}+ \frac{1}{r} \sum_{j=2}^{\ell} \alpha_{\ell+1, j} \leq  \frac{1}{r} \sum_{j=1}^{\ell} \alpha_{\ell+1, j}  \leq \frac{1}{r}  \, .
\end{align*}
A similar argument can be used to show that $|a_{\ell+1,j}|\leq 1/r$, $j= 2, \ldots, \ell$. The Theorem follows by induction.
\end{proof}

Consequently, $$R(K)\leq  R^\textup{opt} (K)\leq \frac{1}{\max_{ij} |a_{ij}|}\,  .$$
For those methods such that
$R (K)= 1/{\max_{ij} |a_{ij}|}
$
it is not possible to increase the SSP coefficient by perturbing the method.
This is the case for all known optimal explicit SSP RK methods of orders one through
four, with any number of stages \cite{ketcheson2008}.

For the restricted class of perturbations considered in \cite{ruuth2004},
similar results were obtained in  \cite[Theorems 3.1, 3.4, 3.5 and
3.6]{ruuth2004}. Theorem \ref{Pbound} extends those results, showing that
no improvement in the radius of absolute monotonicity is possible for many
optimal SSP methods, even when more general perturbations are considered.

\subsection{Relations among the Butcher and canonical Shu-Osher representations}

For perturbations of explicit Runge--Kutta methods (or of other methods with
one stage equal to $u_n$) there exists a certain simple transformation that may yield a larger
value of $R(K,\widetilde{\bbA})$ -- and that never yields a smaller value.
\begin{proposition} \label{propfirstcol}
	Let an $s$-stage explicit perturbed Runge--Kutta method be given with coefficients
    $\gamma_r, \alpha^\textup{up}_r, \alpha^\textup{down}_r \ge 0,$ where $r$ is the radius of absolute
    monotonicity of the method.  Consider the perturbed method with coefficients
    \begin{subequations} \label{transform}
    \begin{align}
    	\widehat{\gamma} & = (1, 0 , \ldots, 0)^t\, ,  \\
        \widehat{\alpha}^\textup{up}_{i,1} & = \alpha^\textup{up}_{i,1} + (\gamma_r)_{i}/2 & 2 \le i \le s  \\
        \widehat{\alpha}^\textup{down}_{i,1} & = \alpha^\textup{down}_{i,1} + (\gamma_r)_{i}/2 & 2 \le i \le s .
    \end{align}
    \end{subequations}
    Then the perturbed method with coefficients $(\gamma_r, \alpha^\textup{up}_r, \alpha^\textup{down}_r)$
    and the modified perturbed method with coefficients $(\widehat{\gamma},
    \widehat{\alpha}^\textup{up}, \widehat{\alpha}^\textup{down})$
    correspond to the same RK method $K$.  The modified perturbed method has radius of absolute
    monotonicity at least equal to $r$.
\end{proposition}
\begin{proof}
It is easily seen that the modified method is equivalent to the original
one when $f=\tilde{f}$, so they correspond to the same unperturbed method.
Meanwhile, the transformation never leads to negative coefficients,
so the modified method is a.m. at $r$.
\end{proof}
\begin{remark}Proposition \ref{propfirstcol} is also valid for Runge-Kutta methods whose first row is equal to zero.
\end{remark}

In the Butcher form \eqref{PRunge--Kutta} it is obvious which
perturbed methods $(K,\widetilde{K})$ correspond to a given method $K$.
In the canonical Shu-Osher form it is less obvious.
The following lemma characterizes which methods of the form \eqref{dwRunge--Kutta_cso}
are perturbations of a given method \eqref{cso}.
\begin{lemma}\label{equivprop}
If method \eqref{dwRunge--Kutta_cso} is a perturbation of method \eqref{cso}, then
their coefficients are related as follows:
\begin{subequations} \label{alpharel1}
\begin{align}
 (I-2\alpha^\textup{down}_r)\alpha_r & = (\alpha^\textup{up}_r - \alpha^\textup{down}_r) \label{alpharel1a} \\
 (I-2\alpha^\textup{down}_r)v_r & = \gamma_r.
\end{align}
\end{subequations}
Furthermore, if \eqref{alpharel1} holds and
the perturbation is zero-well-defined,
then \eqref{dwRunge--Kutta_cso} is a perturbation of \eqref{cso}.
\end{lemma}
\begin{proof}
To prove the first part,
take $\tilde{f} = f$ in \eqref{dwRunge--Kutta_cso} to obtain:
\begin{align*}
Y & = \gamma_r u_n + (\alpha^\textup{up}_r + \alpha^\textup{down}_r) Y
                    + (\alpha^\textup{up}_r - \alpha^\textup{down}_r) \frac{h}{r}F.
\end{align*}
Subtract $2\alpha^\textup{down}_r Y$ from both sides to get
\begin{align} \label{zwd}
(I-2\alpha^\textup{down}_r) Y & = \gamma_r u_n + (\alpha^\textup{up}_r - \alpha^\textup{down}_r)
                    \left(Y +  \frac{h}{r}F\right).
\end{align}
Substituting \eqref{cso} in the above gives
\begin{align*}
(I-2\alpha^\textup{down}_r) v_r u_n + (I-2\alpha^\textup{down}_r)\alpha_r \left(Y +  \frac{h}{r}F\right)
    & = \gamma_r u_n + (\alpha^\textup{up}_r - \alpha^\textup{down}_r) \left(Y +  \frac{h}{r}F\right),
\end{align*}
Equating coefficients yields \eqref{alpharel1}.

To prove the second part, assume $I-2\alpha^\textup{down}_r$ is invertible and
write \eqref{alpharel1} as
\begin{subequations} \label{alpharel2}
\begin{align}
 \alpha_r & = (I-2\alpha^\textup{down}_r)^{-1}(\alpha^\textup{up}_r - \alpha^\textup{down}_r) \\
 v_r & = (I-2\alpha^\textup{down}_r)^{-1}\gamma_r.
\end{align}
\end{subequations}
Substitute \eqref{alpharel2} in \eqref{cso}, multiply on the left by
$(I-2\alpha^\textup{down})^{-1}$, and follow the steps above in reverse.
\end{proof}

Lemma \ref{equivprop} does {\em not} imply that the perturbation \eqref{dwRunge--Kutta_cso}
is unique for a given $r$; see Proposition \ref{propfirstcol}.

\begin{remark} \label{invertibility-remark}
The necessity of the zero-well-defined condition in the second part of
Proposition \ref{equivprop} can be seen from the following example.
We take the implicit trapezoidal Runge-Kutta method
\begin{align*}
 \begin{array}{c|cc}
0 & 0 & 0 \\
1/2 & 1/2 & 1/2 \\ \hline
    & 1/2 & 1/2
\end{array}.
\end{align*}
The canonical form \eqref{cso} is then
\begin{align*}
\alpha_r & = \begin{pmatrix}
0 & 0 & 0 \\[1ex]
\frac{r}{r+2} & \frac{r}{r+2} & 0 \\[1ex]
\frac{r}{r+2} & \frac{r}{r+2} & 0
\end{pmatrix}\, ,
& v_r & = \begin{pmatrix} 1 \\[1ex] \frac{2-r}{r+2} \\[1ex] \frac{2-r}{r+2} \end{pmatrix}.
\end{align*}
Then \eqref{alpharel1} is satisfied -- for any $r$ -- by
\begin{align*}
\alpha^\textup{up} = \alpha^\textup{down} & = \begin{pmatrix}
1/3 & 0 & 0 \\
0 & 1/2 & 0 \\
0 & 1/2 & 0
\end{pmatrix}\, , &
\gamma & = \begin{pmatrix} 1/3 \\ 0 \\ 0 \end{pmatrix}.
\end{align*}
However, this method -- which involves a perturbation that is not zero-well-defined
-- is not a perturbation of the original method. \hfill $\square$
\end{remark}

\subsection{Computing optimal perturbations\label{sec:algorithms}}
In this section we present two algorithms to symbolically
or numerically find the optimal perturbed SSP coefficient and a corresponding perturbation of a given RK method.  The first algorithm is proven to approximate the
optimal value to any accuracy, contingent on the computational solution of
linear program subproblems.  It is only valid for explicit perturbations.  The second
algorithm is analytical and exact, and valid for both explicit and implicit
methods, but it is not proven to give the optimal value.  The results of the
two algorithms coincide (to high precision) for all explicit methods on which
we have tested them.

\subsubsection{Provably optimal algorithm for explicit perturbations}
In the foregoing, we have shown that finding an optimal perturbation
consists of determining the largest $r$ such that there exists a splitting
satisfying \eqref{alpharel1} with positive coefficients.
Note that the range of
values for which a method $(K,\widetilde{K})$ is absolutely monotonic is
always the interval $[0, R(K,\widetilde{K})]$.  Therefore, one way to
find the largest $r$ is to devise a method for testing for
a given $r$ whether there exists a perturbation $\widetilde{K}$ such that $R(K,\widetilde{K})\ge r$.
For given method \eqref{Runge--Kutta} and value of $r$,
the system of equations \eqref{alpharel1} together with the inequalities
$\alpha^\textup{up}_r, \alpha^\textup{down}_r \ge 0$
constitutes a linear programming (LP) feasibility problem.
The following theorem is an immediate consequence of Lemma \ref{equivprop}.

\begin{theorem}
Let an $s$-stage RK method $K$ and a positive number $r$ be given.  There exists a
perturbation $\widetilde{\bbA}$ with $R(K,\widetilde{\bbA})\ge r$ if and only if there exists
an $(s+1)\times(s+1)$ matrix $\alpha^\textup{down}_r$ such that $(I-2\alpha^\textup{down}_r)$ is regular and
the following componentwise inequalities hold:
\begin{subequations} \label{LP}
\begin{align}
 (I-2\alpha^{down}_r)\alpha_r + \alpha^{down}_r & \ge 0 \\
 (I-2\alpha^{down}_r)v_r & \ge 0 \\
 \alpha^\textup{down}_r \ge 0.
\end{align}
\end{subequations}
\end{theorem}
The linear program \eqref{LP} can be solved by standard LP solvers.
By embedding this solution in a one-dimensional root-finding algorithm,
optimal perturbations can be found.  An algorithm based on bisection follows.
\begin{algorithm}\caption{Optimal explicit perturbation}
\label{alg:explicit-LP}
\begin{algorithmic}
\INPUT $K$
\State $r_\textup{max} := 1/\max|a_{ij}|, r_\textup{min} := 0$.
\While {$r_\textup{max} - r_\textup{min} > \epsilon$}
	\State $r=\frac{r_\textup{max} + r_\textup{min}}{2}$.
	\State Compute the coefficient matrices $\alpha_r, v_r$ using \eqref{canonical-coefficients}.
	\State Solve the LP given by \eqref{LP}.
    \If {it is feasible}
    	\State $r_\textup{min}:=r$
    \Else
    	\State $r_\textup{max}:=r$
    \EndIf
\EndWhile \\
\Return $r_\textup{min}$
\end{algorithmic}
\end{algorithm}

Assuming the solution of the LP is correct, the algorithm provably
finds an optimal explicit perturbation.  However,
for implicit perturbations the LP solver may converge to a solution
(like the method in Remark \ref{invertibility-remark} above)
for which $I-2\alpha^\textup{down}_r$ is singular.

\subsubsection{Iterated splitting algorithm}
We next investigate how to choose $\alpha^\textup{up}_r, \alpha^\textup{down}_r$ directly so as
to find a perturbation with radius of a.m. at least $r$.
The following result suggests an approach.
\begin{lemma} \label{signlemma}
Given an explicit Runge--Kutta method \eqref{cso},
let $\alpha^\textup{up}_r\ge0, \alpha^\textup{down}_r\ge0$ denote coefficients of a
zero-well-defined perturbation of \eqref{cso}.
Then there exist matrices $\alpha^+\ge 0$, $\alpha^-\ge 0$ such that
\begin{subequations} \label{alphar}
\begin{align}
    \alpha^\textup{up}_r     & = (I+2\alpha^-)^{-1} \alpha^+ \\
    \alpha^\textup{down}_r   & = (I+2\alpha^-)^{-1} \alpha^-,
\end{align}
\end{subequations}
and $\alpha_r=\alpha^+ - \alpha^-$.
\end{lemma}
\begin{proof}
Since the perturbation is zero-well-defined, we can define
\begin{subequations} \label{a1a2}
\begin{align}
    \alpha^+ & = (I-2\alpha^\textup{down})^{-1} \alpha^\textup{up}_r \\
    \alpha^- & = (I-2\alpha^\textup{down})^{-1} \alpha^\textup{down}_r.
\end{align}
\end{subequations}
Then, by \eqref{alpharel1a}, $\alpha_r=\alpha^+-\alpha^-$.  Furthermore, since $I-2\alpha^\textup{down}_r$ is an $M$-matrix,
we have $\alpha^+\ge0$ and $\alpha^-\ge0$.
Solving \eqref{a1a2} for $\alpha^\textup{up}_r,\alpha^\textup{down}_r$ gives \eqref{alphar}.
\end{proof}

In the next algorithm we use the following notation:
\begin{align*}
((x)^+)_{ij} & = \begin{cases}
x_{ij} & \mbox{ if } x_{ij}\ge0 \\
0 & \mbox{ if } x_{ij}<0.  \end{cases} &
((x)^-)_{ij} & = \begin{cases}
0 & \mbox{ if } x_{ij}\ge0 \\
-x_{ij} & \mbox{ if } x_{ij}<0\, ,   \end{cases}
\end{align*}
and thus $x=(x)^+-(x)^-$ is a sign splitting of matrix $x$, with $(x)^+\geq 0$, $(x)^-\geq 0$.

Given a perturbed Runge-Kutta method \eqref{dwRunge--Kutta_cso} with $\gamma_r=e_1$, where $e_1=(1, 0, \ldots, 0)^t$, and $\alpha^\textup{up}$ or $\alpha^\textup{down}$ containing negative values,
 we construct
\begin{subequations}\label{newpert}
\begin{align}
\tilde \gamma_r&=\left(I+2 \, (\alpha^\textup{up}_r)^-+ 2 \,(\alpha^\textup{down}_r)^-\right)^{-1}  e_1\, ,\\[0.5ex]
\tilde\alpha_r^\textup{up}&=\left(I+2 \,(\alpha^\textup{up}_r)^-+ 2 \,(\alpha^\textup{down}_r)^-\right)^{-1} \left((\alpha^\textup{up}_r)^++ (\alpha^\textup{down}_r)^-\right)\, ,   \\[0.5ex]
\tilde\alpha_r^\textup{down}&=\left(I+2 \,(\alpha^\textup{up}_r)^-+ 2 \,(\alpha^\textup{down}_r)^-\right)^{-1}  \left((\alpha^\textup{up}_r)^-+(\alpha^\textup{down}_r)^+\right)\, ,
\end{align}
\end{subequations}
where
$\alpha^\textup{up}=(\alpha^\textup{up}_r)^+-(\alpha^\textup{up}_r)^-$, $\alpha^\textup{down}_r=(\alpha^\textup{down}_r)^+-(\alpha^\textup{down}_r)^-$, provided that $I+2 (\alpha^\textup{up}_r)^-+ 2 (\alpha^\textup{down}_r)^-$ exists. Using Lemma \ref{equivprop}, it is straightforward to prove that, if method $\alpha^\textup{up}_r$, $\alpha^\textup{down}_r$ is a perturbation of   \eqref{cso}, then \eqref{newpert} is also perturbation of \eqref{cso}.
Next, for explicit methods, we perform transformation \eqref{transform}.
In this way, \eqref{newpert} followed by transformation \eqref{transform} give a perturbation of the form \eqref{dwRunge--Kutta_cso} with $\gamma_r=e_1$,  that we denote by  $\hat \alpha_r^\textup{up}$, $\hat \alpha_r^\textup{down}$. If $\hat \alpha_r^\textup{up}\geq 0$,  $\hat \alpha_r^\textup{down}\geq 0$, then $r$ is an SSP coefficient; otherwise, we can repeat the above process.

The following lemma studies the sign of $(\hat \alpha_r^\textup{up})_{ij}$, $(\hat \alpha_r^\textup{down})_{ij}$ when   $(\alpha^\textup{up}_r)_{ij}<0$ or $(\alpha^\textup{down}_r)_{ij}<0$. For the sake of clarity, we drop the index $r$.

\begin{lemma}\label{lemma_row} We consider a perturbed explicit Runge-Kutta method  with coefficients $\gamma=e_1$, $\alpha^\textup{up}$, $\alpha^\textup{down}$, and
the perturbation   $\hat \alpha^\textup{up}$, $\hat \alpha^\textup{down}$ obtained by computing \eqref{newpert} followed by transformation \eqref{transform}.  Assume that  $j_0\geq 2$ is the first row with negative terms in $\alpha^\textup{up}$ or $\alpha^\textup{down}$. Let $m_0$ be the largest index $m_0\geq 1$ such that  $\alpha^\textup{up}_{j_0,m_0}<0$ or $\alpha^\textup{down}_{j_0,m_0} <0$. Then
\begin{enumerate}
\item For first to $(j_0-1)$-th row, we have: $\hat \alpha^\textup{up}_{i,j}=\alpha^\textup{up}_{i,j}$ and $\hat \alpha^\textup{down}_{i,j}=0$ for $1\leq i\leq j_0-1$, $1\leq j\leq j_0-2$.
\item For the $j_0$-th row, we have:
\begin{enumerate}
\item If  $m_0=1$, then $\hat \alpha^\textup{up}_{j_0,1}<0$ or $\hat \alpha^\textup{down}_{j_0,1}<0$.
\item If $m_0\geq 2$, then,  $\hat \alpha^\textup{up}_{j_0,m_0}\geq 0$  and $\hat \alpha^\textup{down}_{j_0,m_0}\geq 0$.
\item For $1\leq m_0\leq j_0-2$, we have $\hat \alpha^\textup{up}_{j_0,\ell}\geq 0$  and $\hat \alpha^\textup{down}_{j_0,\ell}\geq 0$ for $\ell=m_0+1, \ldots, j_0-1$.
\end{enumerate}
\end{enumerate}

 \end{lemma}

 \begin{proof}  If $j_0$ is the first row with negative terms in $\alpha^\textup{up}$ or $\alpha^\textup{down}$, straightforward computations give that $\hat \alpha^\textup{up}_{i,j}=\alpha^\textup{up}_{i,j}$ and $\hat \alpha^\textup{down}_{i,j}=0$ for $1\leq i\leq j_0-1$, $1\leq j\leq j_0-2$, and
\begin{subequations} \begin{align}
\hat \alpha^\textup{up}_{j_0,1}&=  \alpha^\textup{up}_{j_0,1} -2  \sum_{i=2}^{j_0-1} \left((\alpha^\textup{up}_{j_0,i})^-+ (\alpha^\textup{down}_{j_0,i})^-\right) \,  \alpha^\textup{up}_{i,1}\, , \label{d1} \\\quad
\hat \alpha^\textup{up}_{j_0,\ell}&=  (\alpha^\textup{up}_{j_0,\ell})^++ (\alpha^\textup{down}_{j_0,\ell})^- -2  \sum_{i=\ell+1}^{j_0-1} \left((\alpha^\textup{up}_{j_0,i})^-+ (\alpha^\textup{down}_{j_0,i})^-\right) \, \alpha^\textup{up}_{i,\ell} \, , \quad \ell=2, \ldots, j_0-1\, .\label{dm}
\end{align}
\end{subequations}
and
\begin{subequations}
\begin{align}
\hat \alpha^\textup{down}_{j_0,1}&=  \alpha^\textup{down}_{j_0,1} -2  \sum_{i=2}^{j_0-1} \left((\alpha^\textup{up}_{j_0,i})^-+ (\alpha^\textup{down}_{j_0,i})^-\right) \, \alpha^\textup{down}_{i,1}\, , \label{d11} \\\quad
\hat \alpha^\textup{down}_{j_0,\ell}&=  (\alpha^\textup{up}_{j_0,\ell})^-+ (\alpha^\textup{down}_{j_0,\ell})^+  -2  \sum_{i=\ell+1}^{j_0-1} \left( (\alpha^\textup{up}_{j_0,i})^-+ (\alpha^\textup{down}_{j_0,i})^-\right) \, \alpha^\textup{down}_{i,\ell} \, , \quad \ell=2, \ldots, j_0-1\, .\label{dm1}
\end{align}
\end{subequations}
Let $m_0$ be the largest index $m_0\geq 1$ such that  $\alpha^\textup{up}_{j_0,m_0}<0$ or $\alpha^\textup{down}_{j_0,m_0} <0$. In this case,
$\alpha^\textup{up}_{j_0,i}\geq 0$, $\alpha^\textup{down}_{j_0, i}\geq 0$ for $i=m_0+1, \ldots, j_0-1$, and  thus
$$(\alpha^\textup{up}_{j_0,i})^+=\alpha^\textup{up}_{j_0,i}\, , \quad (\alpha^\textup{down}_{j_0,i})^+=\alpha^\textup{down}_{j_0,i}\, , \quad (\alpha^\textup{up}_{j_0,i})^-=(\alpha^\textup{down}_{j_0,i})^-=0\, ,\qquad i=m_0+1, \ldots, j_0-1\, . $$
 If $m_0=1$, from \eqref{d1} and \eqref{d11} we get  $\hat \alpha^\textup{up}_{j_0,1} = \alpha^\textup{up}_{j_0,1}$ and $\hat \alpha^\textup{down}_{j_0,1}= \alpha^\textup{down}_{j_0,1}$,  and thus $\hat \alpha^\textup{up}_{j_0,1}<0$ or $\hat \alpha^\textup{down}_{j_0,1}<0$.  If $m_0\ge  2$, from \eqref{dm} and \eqref{dm1} we get
 $$\hat \alpha^\textup{up}_{j_0,m_0} = (\alpha^\textup{up}_{j_0,m_0})^++ (\alpha^\textup{down}_{j_0,m_0})^- \geq 0\, , \qquad \hat \alpha^\textup{down}_{j_0,m_0} =  (\alpha^\textup{up}_{j_0,m_0})^-+(\alpha^\textup{down}_{j_0,m_0})^+ \geq 0\, .
 $$
Finally, for $1\leq m_0\leq j_0-2$, from  \eqref{dm} and \eqref{dm1} we get that, for $\ell=m_0+1, \ldots, j_0-1$, we have
$$\hat \alpha^\textup{up}_{j_0,\ell} = (\alpha^\textup{up}_{j_0,\ell})^+ \geq 0\, ,  \qquad \hat \alpha^\textup{down}_{j_0,\ell} =  (\alpha^\textup{down}_{j_0,\ell})^+ \geq 0\,.
$$
 \end{proof}

Consequently, if matrices $\alpha^\textup{up}_r$ and $\alpha^\textup{down}_r$  contain negative elements in the second or later columns, an iterated construction of perturbations $\hat \alpha_r^\textup{up}$, $\hat \alpha_r^\textup{down}$ removes these negative values obtaining a perturbation with non-negative elements from second column on. However, if in a row $j_0$ we have:
\begin{equation}\label{stop}
  \alpha^\textup{up}_{j_0,1}<0 \, \qquad \hbox{and} \qquad \alpha^\textup{up}_{j_0,\ell}\geq 0 \qquad \ell=2, \ldots, j_0-1\, ,
\end{equation}
or
\begin{equation}\label{stop2}
  \alpha^\textup{down}_{j_0,1}<0 \, \qquad \hbox{and} \qquad \alpha^\textup{down}_{j_0,\ell}\geq 0 \qquad \ell=2, \ldots, j_0-1\, ,
\end{equation}
 the new perturbation $\hat \alpha_r^\textup{up}$, $\hat \alpha_r^\textup{down}$ will also contain negative elements in the first column.

We now give Algorithm \ref{alg:2} to determine whether there exists a perturbation
with a.m. radius $r$ for a given method.
\begin{algorithm}\caption{Existence of a perturbation with radius $r$}
\label{alg:2}
\begin{algorithmic}
\INPUT $r, K$
    \State Compute the coefficient matrices $\alpha_r, v_r$ using \eqref{canonical-coefficients}.

    \State Set $\alpha^\textup{up} = \alpha_r$ and $\alpha^\textup{down} = 0$.
	\While {$\alpha^\textup{up}$ or $\alpha^\textup{down}$ has any negative entries}
	    \State If $K$ has a zero row, perform the transformation \eqref{transform}.
	    \State If $\alpha^\textup{up}, \alpha^\textup{down} \ge 0$, stop.  This is a feasible perturbation.
        \State If condition \eqref{stop} or \eqref{stop2} hold, stop. A feasible perturbation cannot be found.
        \State Set $\alpha^- = (\alpha^\textup{up})^- + (\alpha^\textup{down})^+$ and
    		  $\alpha^+ = (\alpha^\textup{up})^+ + (\alpha^\textup{down})^-$
        \State Compute a new splitting:
    		\begin{align*}
            	\alpha^\textup{up} & = \left(I + 2 ((\alpha^\textup{up})^- + (\alpha^\textup{down})^-)\right)^{-1}\alpha^+ \\
                \alpha^\textup{down} & =\left(I + 2((\alpha^\textup{up})^- + (\alpha^\textup{down})^-)\right)^{-1}\alpha^-
            \end{align*}
\EndWhile \\
\end{algorithmic}
\end{algorithm}

\begin{remark}
This approach seems to lead to optimal splittings for all the explicit
methods on which we have tested it.  However, for all implicit methods
we have tested, it fails to increase the radius of absolute monotonicity at all.
Even for explicit methods, we have no proof that it's optimal because one could use
$(\alpha_r)^+ + \delta, (\alpha_r)^- + \delta,$
in place of $(\alpha_r)^+, (\alpha_r)^-$, where $\delta$ is any non-negative matrix.  \hfill $\square$
\end{remark}

\subsection{Examples}\label{sec:examples}
In this section we compute optimal perturbations of some existing methods,
using the algorithms described in the last section.

\begin{table}
\begin{tabular}{ll|l|llllll}
\toprule
Order & Stages & Method & $R(K)$ & $R^\textup{opt}(K)$ & Bound \eqref{Rpert} & Bound \eqref{bound4} & Property C \\
\midrule
1 & 1 & Forward Euler & 1 &          1 &      1 &       1 &                  True \\
\midrule
2 & 2 & Midpoint         & 0 &      0.732 &      1 &   1.414 &                   True \\
  & 2 & Minimal trunc. error & 0.5 &          1 &  1.333 &   1.414 &                   True \\
  & 2 & SSP22 \cite{shu1988}          & 1 &          1 &      1 &   1.414 & True \\
  & 2 & SSP22* \cite{gottlieb2006} & 0.784 &      1.215 &  1.215 &   1.414 & True \\
\midrule
3 & 3 & Heun33 \cite{heun1900neue} & 0 &      0.776 &  1.333 &   1.817 &  False \\
  & 3 & SSP33 \cite{shu1988}   & 1 &          1 &      1 &   1.817 & True \\
\midrule
4 & 4 & RK44 (Kutta)          & 0 &      0.685 &      1 &   2.213 &                  False \\
  & 5 & Merson \cite{merson1957operational} & 0 &      0.242 &    0.5 &   3.309 & False \\
  & 10 &SSP104 \cite{ketcheson2008}       & 6 &          6 &      6 &   8.425 & False \\
\midrule
5 & 6 & Fehlberg \cite{fehlberg1969klassische}     & 0 & 0.057 &  0.125 & 3.727 & False \\
  & 7 & Dormand-Prince \cite{Dormand1980} & 0 &      0.040 &  0.086 &   4.789 & False \\
  & 8 & Bogacki  \cite{Bogacki1996} & 0 &      0.313 &  0.859 &   5.827 & False \\
  & 7 & SSP75 \cite{ruuth2004} & 0 &      1.396 &  1.792 &   4.789 & False \\
  & 8 & SSP85 \cite{ruuth2004} & 0 &      1.875 &  1.919 &   5.827 & True \\
  & 9 & SSP95 \cite{ruuth2004} & 0 &      2.738 &  3.198 &   6.853 & False \\
\midrule
6 & 9 & Calvo \cite{calvo1990} & 0 & 0.021 &  0.059 &   6.265 & False \\
\midrule
8 & 13 & Prince-Dormand \cite{Prince1981} & 0 &      0.013 &  0.059 &   9.212 & False \\
\bottomrule
\end{tabular}
\caption{Properties of some RK methods and their optimal perturbations.  The optimal perturbed
radius of absolute monotonicity was computed by both the linear programming algorithm
and the iterated splitting algorithm; in every case they gave identical results (up to
roundoff errors).  Decimal values have been truncated to the number of digits shown.
\label{tbl:explicit-methods}}
\end{table}

We have computed optimal perturbations for several known explicit methods using the two
algorithms described above.  In all cases, the two algorithms gave the same
values.  It thus seems possible that Algorithm 2
also gives truly optimal results in general, but we do not have a proof.  Properties
of the methods studied are given in Table \ref{tbl:explicit-methods}.  Several interesting
facts are evident:
\begin{itemize}
	\item For all optimal SSP methods (up to order four), perturbation cannot yield a
          larger coefficient.  This is evident already from the bound \eqref{Rpert}.
          For all other methods, some improvement is achieved.
	\item Consistent with Theorem \ref{Pbound}, for every method considered,
    	  it is possible to achieve $R^\textup{opt}>0$ by some perturbation.
    \item The simple bound \eqref{Rpert} predicts the optimal coefficient to within a factor of three in
    		every case.
    \item The methods SSP75, SSP85, and SSP95 are optimal methods found in \cite{ruuth2004},
    		with property C.  By considering methods without property C, we obtain slightly
            larger coefficients for perturbations of SSP75 and SSP95.  On the other hand,
            relaxing the column assumption gives no benefit in the case of the SSP85 method.
    \item The values found have been truncated to three decimal places but are known to greater precistion.  For the 4-stage, order-four method of Kutta, the three-digit value of $R^\textup{opt}(K)$ given in the table matches the value found by Shu and Osher.  However, the exact (irrational) value is slightly larger and is given in the appendix.
\end{itemize}

\section{Conclusions}\label{sec:conclusions}
In this work we have studied SSP coefficients for perturbations of a given explicit Runge-Kutta method. We have considered both the linear and the nonlinear case, and have obtained useful bounds on the threshold factor and on the radius of absolute monotonicity
for perturbed Runge--Kutta methods.
 We have also provided an algorithm for computing
optimal perturbations of explicit Runge-Kutta methods, and given optimal perturbations for many methods from the literature.

This work seems to provide a complete picture for the case of most interest: explicit
methods applied to nonlinear problems.
Nevertheless, some other interesting issues remain unsolved.  These include:
\begin{itemize}
	\item A method to compute optimal perturbations for linear problems.
	\item An algorithm for obtaining optimal splittings of implicit methods.
\end{itemize}
These may be a starting point for future work.

\section{Appendix}

In this section we give additional details on SSP coefficients and optimal perturbations of
second order 2-stage Runge--Kutta methods and the classical 4-stage fourth order Runge--Kutta method.

\subsection{Second order 2-stage methods} \label{sec:2stage-details}

We consider the family of 2-stage second order methods \eqref{2eo2OptLin}. In example \ref{ex:2eo2}   we studied perturbations that increase the SSP coefficient for the linear case. For nonlinear problems, in example \ref{ex_2stage2order},  figure \ref{f22} shows the values of $R^\textup{opt}(K)$ for $\alpha\in[-3,3]$.

In this section,  for each $\alpha$,  we give the expressions for $R^\textup{opt}(K)$ and we show optimal perturbations $\widetilde{\bbA}_{NL}$ such that $R(K, \widetilde{\bbA}_{NL})=R^\textup{opt}(K)$.  It is important to point out the convenience of choosing
$\widetilde{\bbA}_{NL}= \widetilde{\bbA}_{L}$, where  $\widetilde{\bbA}_{L}$ denotes the optimal perturbation for the linear case. In this case, we have not only $R(K, \widetilde{\bbA}_{L})=R^\textup{opt}(K)$ but also ${R}_{\hbox{\footnotesize Lin}}(K,  \widetilde{\bbA}_{NL})={R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)$. The computations required to obtain the results in this section have been done with the symbolic computation program \emph{Mathematica}.

If we denote by $r=R^\textup{opt}(K)$,  we have that
\begin{equation}
r= \begin{cases}
\displaystyle \frac{1}{|\alpha|}\, ,  &  \qquad \displaystyle  \hbox{if} \quad  \alpha \in\left(-\infty, -\frac{1}{2} \left(1+\sqrt{7}\right)\right]\bigcup  \left[\frac{1}{2} \left(-1+\sqrt{7}\right), \infty\right)\,,  \\[3ex]
\displaystyle   \frac{-1+\alpha+\sqrt{ 3 \alpha ^2-2 \alpha +1}
   }{|\alpha| }\, ,  &   \qquad \displaystyle  \hbox{if} \quad  \alpha\in\left(-\frac{1}{2} \left(1+\sqrt{7}\right), 0\right) \bigcup \left(0, \frac{1}{2} \left(-1+\sqrt{7}\right)\right)\, . \end{cases}\label{radiopert22}
\end{equation}
Next we give optimal perturbations $\widetilde{\bbA}_{NL}$.

For $\alpha<0$, we obtain that it is not possible to obtain a perturbation of the form \eqref{2eo2OptLin}  with $\tilde b_2=0$ and $\tilde a_{21}=0$. Consequently, $\widetilde{\bbA}_{NL}\neq \widetilde{\bbA}_{L}$ and we always have that ${R}_{\hbox{\footnotesize Lin}}(K,  \widetilde{\bbA}_{NL})<{R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)$. Optimal perturbations of the form \eqref{2eo2OptLin} for different values of $\alpha<0$ must satisfy the following conditions.
\begin{itemize}
\item For $-\frac{1}{2} \left( 1+\sqrt{7}\right)\leq\alpha <0$, the coefficients $\tilde a_{21}$, $\tilde b_1$ and $\tilde b_2$ in $\widetilde{\bbA}_{NL}$ must satisfy
   $$-\alpha \leq \tilde a_{21}\leq \frac{1-r \, \alpha}{2 \, r}\, , \qquad \tilde b_1=-\frac{r\, \tilde a_{21}}{2\, \alpha } \, , \qquad \tilde b_2=-\frac{1}{2\, \alpha}\, ,
   $$
   where $r=R^\textup{opt}(K)$.
\item For $\alpha\leq-\frac{1}{2} \left( 1+\sqrt{7}\right)$,
   we should have
   $$\tilde a_{21}=-\alpha,\qquad  -\frac{1}{2 \,\alpha }\leq\tilde b_1\leq \frac{-2 \,\alpha ^2-2
 \,  \alpha +1}{4\, \alpha }\, , \qquad  -\frac{1}{2 \,\alpha }\leq
  \tilde b_2\leq \frac{2\, \alpha  \,\tilde b_1-1}{4\, \alpha
   }\, . $$
\end{itemize}
For  $\alpha>0$ we can find optimal perturbations with $\tilde b_2=0$ and $\tilde a_{21}=0$. Coefficient $\tilde{b}_1$ must satisfy the following conditions.
\begin{itemize}
\item For $0<\alpha \leq
   \left(-1+\sqrt{7}\right)/2$,  we have that
   \begin{align}  \qquad \tilde{b}_1=\frac{\sqrt{3 \alpha ^2-2
   \alpha +1}-\alpha}{2 \alpha }  \, . \label{tb1NL}
   \end{align}
Thus there is a unique $\widetilde{\bbA}_{NL}$ of the form \eqref{2eo2OptLin}. In this case, we   have  $R(K)<R(K, \widetilde{\bbA}_{NL})=R^\textup{opt}(K)$.

   \item For $\left(-1+\sqrt{7}\right)/2<\alpha< 1$, we also get $R(K)<R^\textup{opt}(K)$,  but in this case the optimal perturbation  $\widetilde{\bbA}_{NL}$ is not unique. All the perturbations with $\tilde{b}_1$ satisfying
$$   \frac{1-\alpha }{\alpha } \leq \tilde{b}_1\leq
   \frac{2 \alpha ^2-2 \alpha +1}{4 \alpha } \, ,
$$
are optimal. In particular, we can take $\widetilde{\bbA}_{NL}=\widetilde{\bbA}_{L}$. With this choice, $R(K, \widetilde{\bbA}_{L})=R^\textup{opt}(K)=1/\alpha$ and ${R}_{\hbox{\footnotesize Lin}}(K,  \widetilde{\bbA}_{L})={R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)\approx 1.22$. Furthermore,   $\alpha=\left(-1+\sqrt{7}\right)/2$  provides the largest  SSP coefficient within the family of 2-stage second order method (see figure \ref{f22}).

\item For $1\leq \alpha$, we have $R(K)=R^\textup{opt}(K)=1/\alpha$ and the optimal perturbation $\widetilde{\bbA}_{NL}$  is not unique. All the values
$$0\leq \tilde{b}_1\leq
   \frac{2 \alpha ^2-2 \alpha +1}{4 \alpha }
$$
give  optimal perturbations. We can take $\widetilde{\bbA}_{NL}=0$,   but in this case ${R}_{\hbox{\footnotesize Lin}}(K, 0)< {R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)$. A better choice is $\widetilde{\bbA}_{NL}=\widetilde{\bbA}_{L}$. Observe that, for $\alpha=1$, we get the optimal SSP coefficient $R(K)=1$ that cannot be increased by perturbations.
  \end{itemize}

Next, we consider some concrete values of $\alpha$ to show the the expressions of the perturbations. For each value, we give the Butcher tableau of the perturbation and matrices  $\alpha^\textup{up}$ and $\alpha^\textup{down}$ in \eqref{dwRunge--Kutta_cso}.
\begin{itemize}
\item For $\alpha=1/2$ we get method RK2a in \cite{HuKoVa} with $R(K)=0$. With perturbation
\begin{align*}
\widetilde{\bbA}= \begin{pmatrix}
0 & 0 & 0 \\
0 & 0 & 0\\
\tilde{b}_1 & 0 & 0
\end{pmatrix} , \, \alpha^\textup{up}= \begin{pmatrix}
0 & 0 & 0 \\
\tilde{b}_1 & 0 & 0\\
0 & 2 \tilde{b}_1  & 0
\end{pmatrix}
  , \,  \alpha^\textup{down} = \begin{pmatrix}
0 & 0 & 0 \\
0 & 0 & 0\\
1-2\tilde{b}_1  & 0 & 0
\end{pmatrix}
 , \, \gamma = \begin{pmatrix}
1\\
1-\tilde{b}_1\\
0
\end{pmatrix},
\end{align*}
where $\tilde{b}_1=\frac{1}{2} \left(\sqrt{3}-1\right)$, we get  $R(K,\widetilde{\bbA}) ={R}_{\hbox{\footnotesize Lin}}(K, \widetilde{\bbA})= \sqrt{3}-1$.
\item For $\alpha=2/3$,  we have a nontrivial SSP coefficient $R^\textup{opt}(K)=1/2$, but we can increase this value to $R(K,\widetilde{\bbA}_1)=R^\textup{opt}(K)=1$ with perturbation
\begin{align*}
\widetilde{\bbA}_1= \begin{pmatrix}
0 & 0 & 0 \\
0 & 0 & 0\\
1/4 & 0 & 0
\end{pmatrix}
\, , \, \alpha^\textup{up}= \begin{pmatrix}
0 & 0 & 0 \\
2/3 & 0 & 0\\
0 & 3/4 & 0
\end{pmatrix}
\, , \,  \alpha^\textup{down} = \begin{pmatrix}
0 & 0 & 0 \\
0 & 0 & 0\\
1/4 & 0 & 0
\end{pmatrix}
\, , \, \gamma  = \begin{pmatrix}
1\\
1/3\\
0
\end{pmatrix}\, .
\end{align*}
For this perturbation, $R(\phi_K)={R}_{\hbox{\footnotesize Lin}}(K, \widetilde{\bbA}_1)=1$. We can take $\gamma=(1,0,0)^t$ by modifying the first column of $\alpha^\textup{up}$ and $\alpha^\textup{down}$ according to \eqref{transform},
\begin{align*}
\widetilde{\bbA}_2= \begin{pmatrix}
0 & 0 & 0 \\
1/6 & 0 & 0\\
3/8 & 0 & 0
\end{pmatrix}\, , \, \alpha^\textup{up}= \begin{pmatrix}
0 & 0 & 0 \\
5/6 & 0 & 0\\
0 & 3/4 & 0
\end{pmatrix}
\, , \,  \alpha^\textup{down} = \begin{pmatrix}
0 & 0 & 0 \\
1/6 & 0 & 0\\
1/4 & 0 & 0
\end{pmatrix}
\, , \, \gamma = \begin{pmatrix}
1\\
0\\
0
\end{pmatrix}.
\end{align*}

\item As it has been pointed out above, the largest value in the $\alpha$-family of 2-stage second order schemes is $R^\textup{opt}(K)=(1+\sqrt{7})/3$ and it is obtained for $\alpha=(\sqrt{7}-1)/2$. The perturbation is  of the form \eqref{2eo2OptLin} with $\tilde{b}_1= \left(\sqrt{7}-2\right)/2$, and
\begin{align*}
\alpha^\textup{up}= \begin{pmatrix}
0 & 0 & 0 \\
1  & 0 & 0\\
0 & \frac{1}{9} \left(4+\sqrt{7}\right)  & 0
\end{pmatrix}
\, , \,  \alpha^\textup{down} = \begin{pmatrix}
0 & 0 & 0 \\
0 & 0 & 0\\
\frac{1}{9} \left(5-\sqrt{7}\right)  & 0 & 0
\end{pmatrix}
\, , \, \gamma_r = \begin{pmatrix}
1\\
0 \\
0
\end{pmatrix}
\end{align*}
This is the perturbation obtained in \cite[Table V]{gottlieb2006} by numerical search in the class of perturbations considered in \cite{gottlieb2006}.
\end{itemize}

\subsection{Classical fourth order 4-stage method}\label{4eo4L}

For nonlinear problems, applying the analysis above, we find that the optimal
perturbation of the classical method has SSP coefficient given by
the real root of $x^3+2 x^2+4 x-4=0$, which is approximately
$R^\textup{opt}(K)\approx 0.685016$.  The corresponding perturbation is
{\em not unique.}  For instance, we can take $\gamma_r=(1,0,0,0,0)$,
and all entries of $\alpha^\textup{down}_r$ equal to zero except
\begin{align}
	(\alpha^\textup{down}_r)_{31} & = \frac{r^2}{4}   & (\alpha^\textup{down}_r)_{42} & = \frac{r^2}{2},
\end{align}
where $r=R^\textup{opt}(K)$.
However, there exist other optimal perturbations with additionally $(\alpha^\textup{down}_r)_{42}=\epsilon$ where $0\le \epsilon \le 0.782$.

We remark that nearly-optimal perturbations for this method are given
in \cite[p. 448]{shu1988} and \cite{Higueras2010}.
Interestingly, these different perturbed methods have different values of
${R}_{\hbox{\footnotesize Lin}}(K, \widetilde{\bbA})$.