\section{Introduction}

Time-harmonic solutions to the wave equation can be computed in the
frequency or in the time domain. In the frequency domain, a discrete
version of the Helmholtz equation is solved. In the time domain, the
periodic solutions to a discrete wave equation with time-harmonic
forcing term are sought. Frequency domain methods involve less degrees
of freedom. However, the indefinite linear systems resulting from
discretizing the Helmholtz equation are often difficult to
solve. Time-domain methods can be attractive because they require
relatively little memory and are easy to implement if a time-domain
solver is available. The introduction of
\cite{appelo2019waveholtz} contains a recent overview of time-harmonic
wave equation solvers.

Time-domain methods are based on the correspondence between Helmholtz
and wave equations. In this paper we will use as example the damped
wave equation
\begin{equation} \label{eq:damped_wave_equation}
  \frac{1}{c^2} \pdpd{^2 u}{t^2} + R \pdpd{u}{t} - \Delta u  = f ,
\end{equation}
where $u = u(t,x)$ is the wave field, $f = f(t,x)$ is the forcing
term, $\Delta$ is the Laplacian, $c(x)$ is the spatially dependent
wavespeed, $R$ is a spatially dependent damping coefficient, and
$x$ is in some domain $\Omega$ with Dirichlet and/or Neumann
boundary conditions. Let $u(x,t), f(x,t)$ and $U(x),
F(x)$ be related by
\begin{equation} \label{eq:time_harmonic_u_and_f}
  u(x,t) = e^{i \omega t} U(x) , \qquad \text{ and }
  f(x,t) = e^{i \omega t} F(x) .
\end{equation}
Then $u$ satisfies a wave equation with forcing term $f$
if and only if $U$ satisfies the following Helmholtz equation with
forcing term $F$
\begin{equation} \label{eq:damped_Helmholtz_equation}
  - \Delta U  - \frac{\omega^2}{c^2} U + i \omega R U = F .
\end{equation}
This equation is supplemented with Dirichlet and/or Neumann boundary
conditions that carry over from those for
(\ref{eq:damped_wave_equation}).

We briefly review some time-domain approaches.
The most basic time-domain method is derived from the
limiting-amplitude principle. This principle states that solutions
$u(x,t)$ to (\ref{eq:damped_wave_equation}) with zero initial
conditions and forcing term
\begin{equation} \label{eq:basic_idea_forcing_term}
  f(x,t) = e^{i \omega t} F(x)
\end{equation}
satisfy
\begin{equation} \label{eq:limiting-amplitude_principle}
  u(x,t)  = e^{i \omega t} U(x) + O(1) , \qquad t \to \infty
\end{equation}
under certain conditions on the problem, see
\cite{EncyclopediaOfMath_LimitingAmplitudePrinciple} and references
therein.  Thus, if $u(x,t)$ is a (numerical) solution to this initial
boundary-value problem, and $T$ is some large time, measured in
periods, then an approximate solution to the Helmholtz equation is
given by
\begin{equation} \label{eq:limiting_amplitude_approx_sol}
  e^{- i 2 \pi T} u(x, 2\pi \omega^{-1} T ) .
\end{equation}
We will call this the limiting-amplitude approximate solution for
time $T$.
A more advanced method is the
exact controllability method \cite{bristeau1998controllability}. In
this method the periodicity of the solutions is enforced using
optimization. The starting value for the optimization procedure is
typically some partially converged limiting-amplitude solution.
Recently more insights in and some improvements to this method were
obtained \cite{grote2019controllability}, and its parallel
implementation was studied \cite{grote2020parallel}.
In \cite{appelo2019waveholtz} an optimization approach called
WaveHoltz was introduced. This method again uses a form of
optimization but with a different optimization functional.

An important feature of the methods just described is that they approximate
periodic solutions of a given discrete wave equation. This has
two consequences that are in general not desirable. First, the results
will be negatively affected by both spatial and temporal
discretization errors, while solutions to discrete Helmholtz equations
only have spatial discretization errors.  Secondly, the method is
limited to situations in which a time-domain scheme is available, and
not (directly) applicable if one only has a discretization of
(\ref{eq:damped_Helmholtz_equation}) available, or perhaps only a
linear system with similar properties.

In this paper, we will address both of these shortcomings by
developing a new time-domain solver for discrete
Helmholtz equations. In fact the method can be applied to any linear
system
\begin{equation} \label{eq:discrete_damped_Helmholtz}
  H U = F ,
\end{equation}
where $H$ is a complex $N \times N$ matrix, such that
\begin{equation} \label{eq:Re_H_requirement}
  \text{there is $c>0$ such that $\Re H + c I$ is symmetric positive semidefinite} ,
\end{equation}
 and
\begin{equation} \label{eq:Im_H_requirement}
  \text{$\Im H$ is symmetric positive semidefinite}
\end{equation}
In equation (\ref{eq:discrete_damped_Helmholtz})
$F$ is a vector in ${\mathbb{C}}^N$ and
$U$ is the unknown, also in ${\mathbb{C}}^N$.
We next describe the steps involved in the construction of  the new
time-domain solver.

First an $N \times N$ system of second order ODE's
\begin{equation} \label{eq:general_semi-discrete_equation}
  \pdpd{u^2}{t^2} + A u + B \pdpd{u}{t}  = f
\end{equation}
and a frequency parameter $\omega$ are constructed.  We will look for
time-harmonic solutions with frequency $\omega$ of the system
(\ref{eq:general_semi-discrete_equation}). Note that $\omega$ is in
general {\em not} the physical frequency parameter used to derive
(\ref{eq:discrete_damped_Helmholtz}).  It is a computational
parameter, chosen together with the matrices $A$ and $B$, and it depends
on $H$ in a way to be specified.  Time-harmonic functions
$u = e^{i\omega t} U$ and $f = e^{i \omega t} F$, with
$U, F \in {\mathbb{C}}^N$, satisfy (\ref{eq:general_semi-discrete_equation}) if
and only if
\begin{equation} \label{eq:time-harmonic_from_semi-discrete}
  ( - \omega^2 I + i \omega B  + A) U = F .
\end{equation}
Therefore we will choose $A,B$ and $\omega$ such that
\begin{equation} \label{eq:choice_of_ABomega_intro}
  H = -\omega^2 I + i \omega B + A .
\end{equation}
The system (\ref{eq:general_semi-discrete_equation}) plays
the role of a semi-discrete wave equation.

Secondly, this system is time-discretized.
This is done in such a way that time-harmonic solutions of the
discrete-time system are exactly those of the continuous-time system.
I.e.\ if $u_n, f_n$ are related to $U,F$ by $u_n = e^{i\omega n \Delta
  t} U$ and
$f_n = e^{i \omega n \Delta t} F$ then $u_n,f_n$ are solutions to the
discrete-time system if and only if $U,F$ satisfy
(\ref{eq:time-harmonic_from_semi-discrete}).
For this purpose we present two modified leapfrog methods.

So by using non-standard finite-difference derivatives,
there are no time-discretization errors for the time-harmonic
solutions.  This is related to ideas from the papers
\cite{holberg1987computational,tam1993dispersion} where it was shown
that errors in finite-difference discretizations can be analysed
precisely in the Fourier domain, and that finite-difference
approximations of derivatives can be designed to minimize such errors
over some range of wavenumbers.
See also \cite{berland2006low} for an optimized time-stepping method.

Having a time-discretization of
(\ref{eq:general_semi-discrete_equation}) at hand, the third step is
to define a map from a right-hand side $F$ in
(\ref{eq:discrete_damped_Helmholtz}) to an approximation for the
solution $U$.  For this we follow the idea of equation
(\ref{eq:limiting_amplitude_approx_sol}) with one modification, which
is the inclusion of a smooth window function in the time-harmonic
forcing term.

It is known that the convergence of limiting-amplitude approximate
solutions can be slow in some cases, e.g.\ in case of resonant wave
cavities. Therefore we will not use the approximate solution operator
directly. Instead we propose to use it as a preconditioner in an
iterative solution method for (\ref{eq:discrete_damped_Helmholtz}),
such as GMRES or BiCGSTAB. This is the fourth and last step of our
construction. This use as a preconditioner can also be called Krylov
acceleration, as the method itself already approximates the true solution.
It can be seen as an alternative to the exact controllability and
WaveHoltz methods.
Because it is used this way, we will call the new approximate solution
operator a {\em time-domain preconditioner}.
The above distinction between frequency- and time-domain methods
appears no longer satisfactory
for this method: It is based on time-domain methodology, but solves a
frequency-domain discrete wave equation.

The term time-domain preconditioner was used before in
\cite{zschiedrich2006advanced}. In that work a related problem was
solved, but the resulting method was substantially different, as the
Schr\"odinger equation was used instead of the wave equation and the
resulting time-domain scheme was implicit, and not exactly
solved.

The behavior of the proposed method will be explained using theory and
examples. It is established theoretically under which conditions the
new time-discretizations are stable and that the approximate solutions
converge to the exact solution of (\ref{eq:discrete_damped_Helmholtz})
in the limit of large time parameter $T$. The numerical examples
confirm the converge of the approximate solutions and show that Krylov
acceleration is indeed useful to accelerate the solution
process. Krylov acceleration makes a difference in particular in the
case that a resonant low-velocity zone is present. Parameters of the
method include the large time parameter $T$ used in the
preconditioner, and a parameter for the window
function. Results depend weakly on these parameters, i.e.\ there is a
large set of suitable parameter choices.

A few situations that do not fit in the classical time-domain setting,
but can be handled with the new solver are as follows. First one can
use it with discretizations that have been designed specifically for
the Helmholtz equation. For example finite-difference discretizations
that minimize dispersion errors such as those from
\cite{BabuskaEtAl1995} for the 2-D case and
\cite{stolk2016dispersion,sutmann2007compact,TurkelEtAl2013} for the
3-D case. Our examples are about this application. One can
also imagine situations where the physical time-domain model is
complicated to simulate, but reduces to a relatively simple Helmholtz
equation in the frequency domain.  A third use case is in the context
of a multigrid method. In multigrid methods the coarsest level system
still has to be solved by another (non-multigrid) method. It is based
on the linear system one started with and on the choice of
multigrid method. In Helmholtz equations the convergence can be very
sensititve to the choice of the coarse level system and it is
recommended to use certain prescribed discretizations
\cite{StolkEtAl2014,stolk2016dispersion}.

However, the usefulness of the method is not restricted to these
cases, and some of the ideas could also be incorporated into other
time-domain methods.

The contents of the remainder of the paper is as follows. In
section~\ref{sec:method}, the construction of the new solver is
described. In section~\ref{sec:examples} examples of this construction
are given in case $H$ results from certain finite-difference
discretizations. Section~\ref{sec:analysis} contains the theoretical
results. After that, section~\ref{sec:numerical_examples} contains the
numerical examples. We conclude the main text with a discussion
section. Appendices contains some material on Fourier analysis used
in section~\ref{sec:analysis}, and some remarks on the compact-stencil
finite-difference discretization described in section~\ref{sec:examples}.

\section{Methods}
\label{sec:method}

In this section we describe in detail the construction of a
time-domain preconditioner for a matrix $H$ satisfying
(\ref{eq:Re_H_requirement}) and (\ref{eq:Im_H_requirement}).
We recall from the introduction
that there are three main steps:
(i) the definition of a suitable second order system of ODE's of the
form (\ref{eq:general_semi-discrete_equation});
(ii) the definition of a suitable time-integration method for this system of
ODE's;
(iii) the definition of a linear map that produces approximate
solutions based on the limiting-amplitude principle.
Step (ii), the time-integration method, will be discussed first, since
the properties of the time-integration method affect the choice of the
system of ODE's (\ref{eq:general_semi-discrete_equation}).
Then steps (i) and (iii) and the
application of the method as a preconditioner are discussed.
Occasionally we will point forward to section~\ref{sec:analysis}, where
some properties of the time-discretizations and the preconditioner are
proven.

\subsection{Frequency-adapted time discretizations of
  (\ref{eq:general_semi-discrete_equation})}
\label{subsec:time-discretization}

The leapfrog or basic Verlet method is a standard method to integrate
equations of the form (\ref{eq:general_semi-discrete_equation}) in
case that $B = 0$. It is obtained, basically, by discretizing the second order
time derivative using standard second order finite differences.
To allow for nonzero $B$, the damping term has to be discretized as well.
A standard way to do this is with central differences
\cite{brunger1984stochastic, schlick2010molecular,
 sandvik2018numerical}. This yields the equation
\begin{equation} \label{eq:central_differences_damped_leapfrog_verlet}
  \frac{1}{\Delta t^2} \left( u_{n+1} - 2 u_n + u_{n-1}
  \right)
  + \frac{1}{2 \Delta t} B \left( u_{n+1} - u_{n-1} \right)
  + A u_n = f_n ,
\end{equation}
from which $u_{n+1}$ can be solved.
Note that $u_n$ denotes the discrete approximation to $u(n \Delta t)$,
and that $u(t) \in {\mathbb{C}}^N$.
This only leads to an explicit method if $B$ is diagonal. To obtain an
explicit method in case $B$ is non-diagonal, the damping term can be
discretized by backward differences, which yields the equation
\begin{equation} \label{eq:backward_differences_damped_leapfrog_verlet}
  \frac{1}{\Delta t^2} \left( u_{n+1} - 2 u_n + u_{n-1}
  \right)
  + \frac{B}{\Delta t} \left( u_n - u_{n-1} \right)
  + A u_n = f_n .
\end{equation}
Since this method is an order less accurate and has stricter CFL
conditions it is only proposed for the case that $B$ is non-diagonal,
in the other case (\ref{eq:central_differences_damped_leapfrog_verlet}) is preferred.

We will formally define the time-integrators resulting from
(\ref{eq:central_differences_damped_leapfrog_verlet}) and
(\ref{eq:backward_differences_damped_leapfrog_verlet}).

\begin{definition} \label{def:Icd_Ibd}
  Let
  \begin{equation} \label{eq:K_L_central_differences}
  \begin{aligned}
    K = {}& \Delta t^2 A
    \qquad {}&
    L = {}& \Delta t B
    \qquad {}&
    g_n = \Delta t^2 f_n .
  \end{aligned}
\end{equation}
  {\em Central differences damped leapfrog} will be defined as the time
  integrator given by
  \begin{equation} \label{eq:cd_integrator}
    u_{n+1} = I_{\rm cd}(u_n,u_{n-1},f_n)
    := \left( I + \tfrac{1}{2}L \right)^{-1}
    \left( (2 - K) u_n - (I - \tfrac{1}{2}L) u_{n-1} + g_n \right) .
  \end{equation}
  {\em Backward differences damped leapfrog} will be defined as the time
  integrator given by
  \begin{equation} \label{eq:bd_integrator}
    u_{n+1}
    = I_{\rm bd}(u_n,u_{n-1},f_n)
    := (2 - K - L) u_n - (I - L) u_{n-1} + g_n  .
  \end{equation}
\end{definition}

The stability of these methods is studied in
section~\ref{sec:analysis}. According to Theorem~\ref{th:stability_cd},
$I_{\rm cd}$ is stable if
\begin{equation} \label{eq:stability_cond_K_L}
  \text{$K$ and $L$ are positive semidefinite},
\end{equation}
and
\begin{equation} \label{eq:stability_cond_cd_upper}
  \text{$4I - K$ is positive definite} ,
\end{equation}
while $I_{\rm bd}$ is stable if
(\ref{eq:stability_cond_K_L}) is satisfied and
\begin{equation} \label{eq:stability_cond_bd_upper}
  \text{$4I - K - 2L$ is positive definite} .
\end{equation}
Conditions (\ref{eq:stability_cond_cd_upper}) and
(\ref{eq:stability_cond_bd_upper}) lead to CFL bounds. This will
be discussed below.

\medskip

As mentioned in the introduction, we look for discretizations such
that the time-harmonic solutions of frequency $\omega$ of the
discrete-time system are exactly those of the continuous time system
(\ref{eq:general_semi-discrete_equation}). Due to discretization
errors this is not the case for
the integrators $I_{\rm cd}$ and $I_{\rm bd}$.  In the next
proposition we will show that the time-harmonic solutions to
(\ref{eq:general_semi-discrete_equation}) satisfy recursions of the
same form as (\ref{eq:central_differences_damped_leapfrog_verlet}) and
(\ref{eq:backward_differences_damped_leapfrog_verlet}), but with
different choices of $A,B$. From these recursions modified schemes can derived
that have the desired property.

\begin{proposition} \label{prop:reproduce_symbol_at_omega}
  Let $u_n$ and $f_n$ be related to $U,F \in {\mathbb{C}}^N$ by
  \begin{equation}
    u_n = e^{i \omega n \Delta t} U , \qquad
    f_n = e^{i \omega n \Delta t} F
  \end{equation}
  and let
  \begin{equation} \label{eq:define_alpha_beta}
    \begin{aligned}
      \alpha = {}& \frac{(\Delta t \, \omega)^2}{2 - 2 \cos (\omega \Delta t)}
      = \frac{(\Delta t \, \omega)^2}{4 \sin( \frac{\omega \Delta
          t}{2})^2}  ,
      & \qquad \text{ and } \qquad
      \beta = {}& \frac{\omega \Delta t}{ \sin( \omega \Delta t)} .
    \end{aligned}
  \end{equation}
  Then $U,F$ satisfy (\ref{eq:time-harmonic_from_semi-discrete}) if
  and only if $u_n$, $f_n$ satisfy
  \begin{equation} \label{eq:mod_central_differences_damped_leapfrog_verlet}
    \frac{1}{\Delta t^2} \left( u_{n+1} - 2 u_n + u_{n-1}
    \right)
    + \frac{1}{2 \Delta t} \tilde{B} \left( u_{n+1} - u_{n-1} \right)
    + \tilde{A} u_n = \alpha^{-1} f_n ,
  \end{equation}
  where
  \begin{equation} \label{eq:redefineABf_cd}
    \tilde{A} = \alpha^{-1} A  , \qquad \text{ and } \qquad
    \tilde{B} = \alpha^{-1} \beta B ,
  \end{equation}
  and if and only if $u_n$ and $f_n$ satisfy
  \begin{equation} \label{eq:mod_backward_differences_damped_leapfrog_verlet}
    \frac{1}{\Delta t^2} \left( u_{n+1} - 2 u_n + u_{n-1}
    \right)
    + \frac{\hat{B}}{\Delta t} \left( u_n - u_{n-1} \right)
    + \hat{A} u_n = \alpha^{-1} f_n ,
  \end{equation}
  where
  \begin{equation} \label{eq:redefineABf_bd}
    \hat{A}  =  \alpha^{-1} A
    - \frac{\beta (1-\cos(\omega \Delta t))}{\alpha \Delta t} B ,
    \qquad \text{ and } \qquad
    \hat{B}   = \alpha^{-1} \beta B .
  \end{equation}

\end{proposition}

\begin{proof}
  To prove the first claim,
  $\tilde{A}$, $\tilde{B}$ and $\tilde{c}$
  will be constructed such that
  \begin{equation} \label{eq:recursion_tilde}
  \frac{1}{\Delta t^2} \left( u_{n+1} - 2 u_n + u_{n-1}
  \right)
  + \frac{1}{2 \Delta t} \tilde{B} \left( u_{n+1} - u_{n-1} \right)
  + \tilde{A} u_n = \tilde{c} f_n ,
  \end{equation}
  if and only if (\ref{eq:time-harmonic_from_semi-discrete}).
  Inserting $u_n = U e^{i n \omega \Delta t}$ into
  (\ref{eq:recursion_tilde}), results in
  \begin{equation}
    \left[ \frac{2 \cos(\omega \Delta t) - 2}{\Delta t^2}
      + \frac{i \, \sin( \Delta t \omega)}{\Delta t} B
      + \tilde{A}
    \right] U e^{i n \omega \Delta t} = \tilde{c} e^{i n \omega \Delta t} F .
  \end{equation}
  Using the definitions of $\alpha$ and $\beta$ to rewrite the
  left-hand side, this is equivalent to
  \begin{equation}
    \left[ - \frac{\omega^2}{\alpha} + i \frac{\omega}{\beta} \tilde{B}
      + \tilde{A} \right] U = \tilde{c} F .
  \end{equation}
  Multiplying by $\alpha$ results in
  \begin{equation}
    \left[ - \omega^2 + i \omega \frac{\alpha}{\beta} \tilde{B}
      + \alpha \tilde{A} \right] U = \alpha \tilde{c} F .
  \end{equation}
  This is equivalent to (\ref{eq:time-harmonic_from_semi-discrete})
  if $\tilde{c} = \alpha^{-1}$ and $\tilde{A}$ and $\tilde{B}$ are defined
  as in (\ref{eq:redefineABf_cd}).

  To prove the second claim,
  $\hat{A}$, $\hat{B}$ and $\hat{c}$
  will be constructed such that
  \begin{equation} \label{eq:recursion_hat}
  \frac{1}{\Delta t^2} \left( u_{n+1} - 2 u_n + u_{n-1} \right)
  + \frac{1}{\Delta t} \hat{B} \left( u_n - u_{n-1} \right)
  + \hat{A} u_n = \hat{c} f_n ,
  \end{equation}
  if and only if (\ref{eq:time-harmonic_from_semi-discrete}).
  Inserting $u_n = U e^{i n \omega \Delta t}$ into
  (\ref{eq:recursion_hat}), results in
  \begin{equation}
    \left[ \frac{2 \cos(\omega \Delta t) - 2}{\Delta t^2}
      + \frac{i \, \sin( \Delta t \omega)}{\Delta t} \hat{B}
      + \frac{1 - \cos( \Delta t \omega)}{\Delta t} \hat{B}
      + \hat{A}
    \right] U e^{i n \omega \Delta t} = \hat{c} e^{i n \omega \Delta t} F .
  \end{equation}
  Using the definitions of $\alpha$ and $\beta$ and multiplying by
  $\alpha$ results in the equivalent equation
  \begin{equation}
    \left[ - \omega^2
      + i \omega \frac{\alpha}{\beta} \hat{B}
      + \frac{\alpha(1 - \cos( \Delta t \omega))}{\Delta t} \hat{B}
      + \alpha \hat{A} \right] U = \alpha \hat{c} F .
  \end{equation}
  This is equivalent to (\ref{eq:time-harmonic_from_semi-discrete})
  if $\hat{c} = \alpha^{-1}$ and $\hat{A}$ and $\hat{B}$ are defined
  as in (\ref{eq:redefineABf_bd}).
\end{proof}

Based on the proposition we define the following time integrators.
The time-harmonic solutions with frequency $\omega$ of these
integrators correspond exactly to time-harmonic solutions of
(\ref{eq:general_semi-discrete_equation}).

\begin{definition} \label{def:Iacd_Iabd}
  {\em Frequency adapted central differences damped leapfrog}
  will be defined as the time integrator given by
  \begin{equation} \label{eq:acd_integrator}
    u_{n+1} = I_{\rm acd}(u_n,u_{n-1},f_n)
    := \left( I + \tfrac{1}{2}L \right)^{-1}
    \left( (2 - K) u_n - (I - \tfrac{1}{2}L) u_{n-1} + g_n \right) ,
  \end{equation}
  where $K$, $L$ and $g_n$ are given by
  \begin{equation} \label{eq:K_L_frequency_adapted_central_differences}
    \begin{aligned}
      K = {}& \frac{\Delta t^2}{\alpha} A ,
      & \qquad
      L = {}& \frac{\beta \Delta t}{\alpha} B ,
      & \qquad
      g_n = {}& \frac{\Delta t^2}{\alpha} f_n
    \end{aligned}
  \end{equation}
  {\em Frequency adapted backward differences damped leapfrog} will be
  defined as the time integrator given by
  \begin{equation} \label{eq:abd_integrator}
    u_{n+1}
    = I_{\rm abd}(u_n,u_{n-1},f_n)
    := (2 - K - L) u_n - (I - L) u_{n-1} + g_n  .
  \end{equation}
  where $K$, $L$ and $g_n$ are given by
  \begin{equation} \label{eq:K_L_frequency_adapted_backward_differences}
    \begin{aligned}
      K  = {}& \frac{\Delta t^2}{\alpha} A
      - \frac{\Delta t \,
        \beta (1-\cos(\omega \Delta t))}{\alpha} B \qquad
      \\

      L = {}& \frac{\beta \Delta t}{\alpha} B ,
      \\

      g_n = {}& \frac{\Delta t^2}{\alpha} f_n
    \end{aligned}
  \end{equation}
\end{definition}

The time integrators $I_{\rm cd}$ and $I_{\rm acd}$ are of the same
form with different choices for $K$ and $L$. Therefore for
$I_{\rm acd}$ the stability conditions are again
(\ref{eq:stability_cond_K_L}) and (\ref{eq:stability_cond_cd_upper}),
but now with $K,L$ as in
(\ref{eq:K_L_frequency_adapted_central_differences}). Similarly, for
$I_{\rm abd}$, the stability conditions are
(\ref{eq:stability_cond_K_L}) and (\ref{eq:stability_cond_bd_upper})
with $K,L$ as in
(\ref{eq:K_L_frequency_adapted_backward_differences}).

\subsection{Choice of the semi-discrete system and the parameters
  $\omega$ and $\Delta t$}
\label{subsec:choose_semidiscrete}

We now look for a second order system of ODE's of the form
(\ref{eq:general_semi-discrete_equation}), and a parameter $\omega$
such that the time-harmonic solutions to
(\ref{eq:general_semi-discrete_equation}) satisfy $H U = F$. Recall
that $\omega$ is the computational frequency parameter, that is chosen
in the construction of the algorithm, and not the physical frequency
used in the underlying Helmholtz equation. We will
also discuss the choice of time integrator ($I_{\rm acd}$ or $I_{\rm
  abd}$) and the choice of the parameter $\Delta t$.
The requirements are that an explicit and stable scheme is obtained.

Because of (\ref{eq:choice_of_ABomega_intro}), we set
\begin{align} \label{eq:A_B_from_H}
  A = {}& \Re H + \omega^2 I , \qquad \text{and}
  \\
  \label{eq:relate_B_ImH}
  B = {}& \omega^{-1} \Im H .
\end{align}
where $\omega$ is still to be determined.

It is required that the time-integrator leads to explicit scheme for
the given matrices $A,B$. By (\ref{eq:relate_B_ImH}) the matrix $B$ is
diagonal if and only if
\begin{equation} \label{eq:ImH_diagonal}
  \text{$\Im H$ is diagonal} .
\end{equation}
Following the remarks below
(\ref{eq:central_differences_damped_leapfrog_verlet})
the frequency adapted central differences method is used if $\Im H$ is
diagonal, and the frequency adapted backward differences method
otherwise.

Next the stability conditions must be considered.
In case $\Im H$ is diagonal these are given in
(\ref{eq:stability_cond_K_L}) and (\ref{eq:stability_cond_cd_upper}).
Because $\Im H$ is positive semidefinite, $L$ is also positive
semidefinite. To ensure that $K$ is positive semidefinite, we set
\begin{equation} \label{eq:define_omega_ImH_diagonal}
  \omega = \sqrt{ - \lambda_{\rm min}(\Re H)} .
\end{equation}
(or to a lower bound for $\sqrt{ - \lambda_{\rm min}(\Re H)}$ if this
value is not exactly known).
The other stability condition (\ref{eq:stability_cond_cd_upper}) is a form of the
well-known CFL condition. It implies that the eigenvalues
$K$ must be less than $4$. For central differences
damped leapfrog integration (not frequency adapted) it implies the condition
\begin{equation}
  \Delta t < \frac{2}{\sqrt{\lambda_{\rm max}(A)}} .
\end{equation}
For the frequency adapted variant it implies, by
(\ref{eq:K_L_frequency_adapted_central_differences}) and
(\ref{eq:define_alpha_beta}) the following condition on $\Delta t$
\begin{equation}
  \sin( \frac{\omega \Delta t}{2} )^2
  < \frac{\omega^2}{\lambda_{\rm max}(A)} ,
\end{equation}
hence
\begin{equation} \label{eq:define_Deltat_ImH_diagonal}
  \Delta t
  <
  \frac{2}{\omega} \arcsin( \frac{\omega}{\sqrt{ \lambda_{\rm max}(A) } } ) .
\end{equation}

The choices and requirements (\ref{eq:A_B_from_H}),
(\ref{eq:define_omega_ImH_diagonal}) and
(\ref{eq:define_Deltat_ImH_diagonal}) define a suitable choice of
parameters of the systems of ODEs and the $I_{\rm acd}$ integrator in
case $\Im H$ is diagonal.

In case that $\Im H$ has nonzero entries outside the diagonal and the $I_{\rm
  abd}$ integrator is used, the stability conditions are
(\ref{eq:stability_cond_K_L}) and (\ref{eq:stability_cond_bd_upper})
with $K,L$ as in
(\ref{eq:K_L_frequency_adapted_backward_differences}).
This results in two conditions that both involve $\omega$ and
$\Delta t$.
In this case it is convenient to use $\omega \Delta t$ and $\Delta t$
as parameters instead of $\omega$ and $\Delta t$. The parameter
$\omega \Delta t$ should be between 0 and $\pi$. Given
$\omega \Delta t$ the following expression for $\omega^2$ can be
derived from the condition that $K$ is positive semidefinite and
equation (\ref{eq:K_L_frequency_adapted_backward_differences})
\begin{equation}
  \omega^2 = - \lambda_{\rm min}(\Re H)
  + \frac{\beta}{\omega \Delta t} ( 1 - \cos(\omega \Delta t))
  \lambda_{\rm max}(\Im H)
\end{equation}
(instead of $\lambda_{\rm min}(\Re H)$ and
$\lambda_{\rm max}(\Im H) $ lower and upper bounds can be used respectively).
From the condition that $4I - K - 2L$ is positive definite we then get
the following scalar condition
\begin{equation} \label{eq:Deltat_bd_full_inequality}
  \frac{ (\omega \Delta t)^2} {\alpha}
  \left( \frac{\lambda_{\rm max}(\Re H)}{\omega^2} + 1 \right)
  + \frac{\beta}{\alpha} (1 + \cos(\omega \Delta t))
  \frac{\lambda_{\rm max}(\Im H)}{ \omega^2} < 4 .
\end{equation}
The following is a stronger inequality than
(\ref{eq:Deltat_bd_full_inequality})
\begin{equation}
\left( \frac{\lambda_{\rm max}(\Re H)}{-\lambda_{\rm min}(\Re H)} + 1
\right) (\omega \Delta t)^2
+ 2 \frac{\lambda_{\rm max}(\Im H)}{-\lambda_{\rm min}(\Re H)}
\omega \Delta t < 4 .
\end{equation}
From here a value of $\omega \Delta t$ can be obtained that satisfies
the conditions by solving a simple quadratic equation.
If a larger value of $\omega \Delta t$ is desired, one can look numerically
for a value as large as possible for which
(\ref{eq:Deltat_bd_full_inequality})
is still satisfied.

\subsection{Time-domain approximate solution operators}
\label{subsec:time-domain_approx_sol}

In this section the time-domain approximate solution operator will be
defined. We also show that a complex approximate solution can be computed
by solving a real time-domain wave equation. This appears to be a standard trick
in the field.
The main novelty is that the formula for the time-harmonic forcing
term (\ref{eq:basic_idea_forcing_term}) is modified so that the
forcing is turned on gradually.

The time-domain approximate solution operators will be denoted by
$S^{\rm acd}_T$ and $S^{\rm abd}_T$, depending on which of the
integrators $I_{\rm acd}$ and $I_{\rm abd}$ is used. They are defined
as follows.

\begin{definition}
  Let $P = {\rm acd}$ or $P = {\rm abd}$,
  and assume that $I_P$ is an integrator for the system
  (\ref{eq:general_semi-discrete_equation})
  where $A,B$ and $\omega$ satisfy (\ref{eq:A_B_from_H}).
  Let $\chi$ be a $C^\infty$ function that is non-decreasing, equal to
  zero for $s \le 0$, and equal to one for $s \ge 1$
  and let $T$ be a positive real constant, such that
  $n_{\rm steps} := 2\pi \omega^{-1} T / \Delta t$ is an integer.
  For $F \in {\mathbb{C}}^N$, let
  \begin{equation}
    f_n = f(n \Delta t) , \qquad
    f(t) = \chi(\frac{t}{2\pi \omega^{-1} T} ) e^{i \omega t} F .
  \end{equation}
  The {\em time-domain approximate solution operator for $H$} associated with the
  integrator $I_P$ is the linear map $S_T^P : {\mathbb{C}}^N \to {\mathbb{C}}^N$ defined by
  \begin{equation}
    S_T^P F = e^{- i 2\pi T} u_{n_{\rm steps} } ,
  \end{equation}
  where $u_n$, $n=0,1, \ldots, n_{\rm steps} $ is given by
  \begin{equation}
    u_{n+1} = I_P (u_n,u_{n-1}, f_n) , \qquad u_0 = 0 .
  \end{equation}
\end{definition}

In section~\ref{sec:analysis} the convergence of $S_T^{\rm acd} F$ and
$S_T^{\rm abd} F$ to $H^{-1} F$ will be established under the
assumption that requirements for stability of the the time integrators
discussed in subsection~\ref{subsec:choose_semidiscrete} are
satisfied. Section~\ref{sec:numerical_examples} contains numerical
examples for $S_T^{\rm acd}$.

There is still a lot of freedom to choose the window function $\chi$.
In the numerical examples we will introduce therefore an additional
parameter $\rho$, $0 < \rho \le 1$, such that the window function
is 1 on $[\rho 2\pi \omega^{-1} T , 2\pi \omega^{-1} T]$ and positive but strictly less than one on
$(0,\rho 2\pi \omega^{-1} T )$.  We will fix a function $\chi_1$ and set
\begin{equation}
  \chi(s) = \chi_\rho(s) = \chi_1(\rho^{-1} s) ,
\end{equation}
The function $\chi_1$ that is used is chosen as a simple sine square
window
\begin{equation}
  \chi_1(s) =
  \left\{
    \begin{array}{ll}
      0 & \text{if $s \le 0$}
      \\
      \sin( \pi s /2)^2 & \text{if $0<s<1$}
      \\
      1 & \text{if $s \ge 1$}.
    \end{array}
  \right.
\end{equation}
The parameter $\rho$ will be called the window parameter.

We next show that it is sufficient to solve a real time-domain wave
problem to compute the limiting-amplitude approximate solution. The
argument is given in the continuous case, but is applicable
equally well in the discrete case.

Let $F$ be a complex right-hand side for the Helmholtz equation and
$u(t,x)$ be the solution to (\ref{eq:damped_wave_equation})
with right-hand side $f(t,x) = F(x) e^{i \omega t}$.  Assuming the
limiting-amplitude principle holds, cf.\
(\ref{eq:limiting-amplitude_principle}),
an approximate solution to the Helmholtz equation is given by
\begin{equation}
  U(x) \approx e^{-i\omega t} u(t,x) , \qquad \text{for some large
    $t$} .
\end{equation}
The field $\Re u(t,x)$ can be determined by solving the real wave
equation with real forcing term, i.e.\ with forcing term
\begin{equation}
  \Re F(x) e^{i \omega t}
  = \cos(\omega t) \Re F(x) -\sin(\omega t) \Im F(x) .
\end{equation}
From (\ref{eq:limiting-amplitude_principle}) it follows that
\begin{equation}
  \Re u(t,x) = \Re U(x) \cos(\omega t) - \Im U(x) \sin(\omega t) +
  o(1) , \qquad t \to \infty .
\end{equation}
Approximations to $\Re U(x)$ and $\Im U(x)$ can hence be obtained from $\Re
u(t,x)$ by
\begin{equation} \label{eq:real_evaluation}
  \begin{aligned}
    \Re U(x) \approx {}& \Re u(t,x) ,
    \qquad t = T \frac{2 \pi}{\omega}
    \\
    \Im U(x) \approx {}& \Re u(t,x) ,
    \qquad t = \left( T - \frac{1}{4} \right) \frac{2 \pi}{\omega}
  \end{aligned}
\end{equation}
with $T$ a large integer. Therefore real time-domain simulation is
sufficient.

In the examples $\Delta t$ is chosen such that the period
$2 \pi \omega^{-1}$ is an integer multiple of $ 4 \Delta t$.
The approximations for $\Re U$ and $\Im U$ in
(\ref{eq:real_evaluation}) are then easy to compute.

\subsection{Time-domain preconditioned GMRES}

The approximate Helmholtz solver can be used as a preconditioner for
iterative methods like GMRES of BiCGSTAB.
Without preconditioning, the system to be solved is  $H U = F$.
Applying left-preconditioning means that instead the system
\begin{equation}
  P H U = P F
\end{equation}
is solved, where $P = S_T^{\rm acd}$ or $S_T^{\rm abd}$.
The right-preconditioned system is
\begin{equation}
  H P V = F .
\end{equation}
The vector $V$ is solved from this system and the solution to the
original problem is then given by $U = P V$.  We propose a solution method
where GMRES is applied to the left-preconditioned system.

\section{Examples}
\label{sec:examples}

The examples we consider are finite-difference discretizations
of the damped Helmholtz equation
\begin{equation} \label{eq:damped_Helmholtz_example}
  - \Delta U - k^2 \left( 1 - i \frac{ \tilde{R} }{\pi} \right) U = F
  .
\end{equation}
Here $\tilde{R}$ is the spatially dependent damping in units of
damping per cycle. We start with standard second order differences and
then apply the method to the discretization from
\cite{stolk2016dispersion}.

The method is not limited to finite-difference
discretizations. Finite-element discretizations can also lead to
linear systems $H U = F$ with $H$ satisfying
(\ref{eq:Re_H_requirement}) and (\ref{eq:Im_H_requirement}).

We will see that in our method, the discrete system is a
discretization of a modified PDE, see (\ref{eq:modified_PDE}) below,
and that the CFL bound for this modified PDE is independent of the
velocity.

\subsection{Second order finite differences}
It is instructive to start with a simple second order
finite-difference discretization. For readability we describe the
two-dimensional case. The degrees of freedom will be denoted by
$U^{(i,j)}$ (in two dimensions), where $i,j$ are in some rectangular
domain $D \subset {\mathbb{Z}}^2$. The matrix $H$ is defined by the equation
\begin{equation}
  ( H U )^{(i,j)} =
  h^{-2}  \bigg( 4 U^{(i,j)} - U^{(i-1,j)} - U^{(i+1,j)} - U^{(i,j-1)} -
  U^{(i,j+1)} \bigg) - (k^{(i,j)})^2  \left( 1 - i \frac{ \tilde{R}^{(i,j)} }{\pi} \right) U^{(i,j)}
\end{equation}
where $h$ denotes the grid spacing, and Dirichlet boundary conditions
are assumed, i.e.\ $U^{(i,j)} = 0$ if $(i,j) \notin D$.

In this case $\Im H$ is diagonal, and $\omega$ and $\Delta t$ are
chosen based on the values
\begin{equation}
  \begin{aligned}
  \lambda_{\rm min}(\Re H) = {}& - k_{\rm max}^2
\\
\lambda_{\rm max}(\Re H) = {}& - k_{\rm min}^2 + 8 h^{-2} .
\end{aligned}
\end{equation}
where
\begin{equation}
  k_{\rm max} := \max_{(i,j) \in D} k^{(i,j)} , \qquad \text{ and } \qquad
  k_{\rm min} := \min_{(i,j) \in D} k^{(i,j)} .
\end{equation}
We find that $\omega = k_{\rm max}$
while $\Delta t$ is chosen from (\ref{eq:define_Deltat_ImH_diagonal}) using that
$\lambda_{\rm max}(A) = 8 h^{-2} + k_{\rm max}^2 - k_{\rm min}^2$.
The scheme in the computational time domain, can be written as
\begin{equation} \label{eq:discrete_scheme_2ndOrderFd}
  \begin{aligned}
  {}& \frac{\alpha}{\Delta t^2} (u_{n+1}^{(i,j)} - 2 u_n^{(i,j)} + u_{n-1}^{(i,j)} )
  + \frac{\beta}{2 \Delta t} \frac{(k^{(i,j)})^2 \tilde{R}^{(i,j)}}{\pi k_{\rm max}}
  ( u_{n+1}^{(i,j)} - u_{n-1}^{(i,j)} )
  \\
  {}& + h^{-2}  \bigg( 4 u_n^{(i,j)} - u_n^{(i-1,j)} - u_n^{(i+1,j)} - u_n^{(i,j-1)} -
  u_n^{(i,j+1)} \bigg) + ( k_{\rm max}^2 -(k^{(i,j)})^2 )  u_n^{(i,j)} =
  f_n^{(i,j)} .
\end{aligned}
\end{equation}

The scheme (\ref{eq:discrete_scheme_2ndOrderFd}) differs substantially
from the standard second order FDTD scheme. In fact, it is a
discretization of the PDE
\begin{equation} \label{eq:modified_PDE}
  \pdpd{^2 u}{t^2}  - \Delta u
  + \frac{k^2 \tilde{R}}{\pi k_{\rm max}} \pdpd{u}{t}
  + \left( k_{\rm max}^2- k^2 \right) u = f .
\end{equation}
This is not necessarily the case. If the operator $H$ would
be rescaled, multiplying from the left and the right by a diagonal
matrix with entries $c^{(i,j)}$ on the diagonal, then the resulting
computational time domain scheme would be close to a standard FDTD
scheme.

This leads to a different behavior of the CFL bound and the minimum
number of timesteps per period. In standard schemes the CFL bound is
chosen based on the maximum velocity, but the grid spacing based on
the minimum velocity, this results in a rough equality.
\begin{equation}
  (\text{\# timesteps per period})
  \lesssim \frac{c_{\rm max}}{c_{\rm min}}
  (\text{min \# gridpoints per wavelength})
\end{equation}
For the new scheme (\ref{eq:discrete_scheme_2ndOrderFd}) we have
\begin{equation} \label{eq:timesteps_per_period_approx_bound_newscheme}
  (\text{\# timesteps per period})
  \lesssim (\text{min \# gridpoints per wavelength}) .
\end{equation}
For standard schemes, the CFL bound is locally suboptimal, where the
velocity is below the maximum. For the
modified PDE (\ref{eq:modified_PDE}) on the other hand, the CFL bound is
locally optimal or near-optimal throughout the domain.

\subsection{An optimized finite-difference method}
\label{subsec:optimized-fd}
Next we consider a discretization with a 27 point cubic
stencil (in 3-D) that was described in
\cite{stolk2016dispersion}.
There exists different variants of such compact stencil methods, see
among others \cite{BabuskaEtAl1995,sutmann2007compact,TurkelEtAl2013}.
For each gridpoint, the number of neighbors that the gridpoint
interacts with is relatively small (when compared e.g.\ to
higher order finite elements). Even so, the dispersion errors of these
schemes are also relatively small, as shown in \cite{stolk2016dispersion},
so that these methods can be used with relatively coarse meshes%
\footnote{In some applications relatively coarse meshes are used to
  save on computation time. Reduction of dispersion errors is
  important in this case. In general, other discretization
  errors can also be present, e.g.\ related to the presence
  of variable coefficients.}.

The equation discretized in \cite{stolk2016dispersion} was the real
part of (\ref{eq:damped_Helmholtz_example}).
We will discuss here the situation with constant $k$, for variable $k$
we refer to appendix~\ref{app:opt_fd}.
The stencil is the 27 point cube that we will denote with
$\{ -1,0,1\}^3$.
For constant $k$, due to symmetry there are four different
matrix coefficients. The values of these coefficients are $h^{-2} f_s(
\frac{h k}{2\pi} )$, for $s=0,1,2,3$ respectively, where the dimensionless
function $f_s$ describes the coefficient as a function of the
dimensionless quantity $\frac{h k}{2\pi} $. To be precise,
\begin{equation}
(\Re H)^{(i,j,l; p,q,r)} =
\left\{
  \begin{array}{ll}
    \frac{1}{h^2} f_{|i-p|+|j-q|+|l-r|}( \frac{h k}{2 \pi} )
    & \text{if $(i-p,j-q,l-r) \in \{-1,0,1\}^3$}
    \\
    0 & \text{otherwise}
  \end{array}
\right.
\end{equation}
In \cite{stolk2016dispersion} the $f_s$ are precomputed functions, in
a class of piecewise polynomial functions, obtained by optimization to
minimize phase errors.
For the 2-D case, in \cite{BabuskaEtAl1995} explicit
expressions for the optimal functions $f_s$, $s=0,1,2$ were derived.
The imaginary part was discretized as above, i.e.\ it was diagonal
with entries
$\frac{(k^{(i,j,l)})^2 \tilde{R}^{(i,j,l)}}{\pi}$.

Because $\Im H$ is diagonal, the integrator $I_{\rm acd}$ is used.
The discretization in \cite{stolk2016dispersion} has in addition the
property that it was second order and such that
\begin{equation}
  \text{$\Re H + k^2 I$ is symmetric positive definite} .
\end{equation}
Thus in case of constant $k$
\begin{equation} \label{eq:omega_k}
  \omega = k .
\end{equation}
The upperbound for $\Re H$ that determines the maximum for $\Delta t$
is given by
\begin{equation} \label{eq:upperbound_HoptimFd}
  \lambda_{\rm max}(\Re H) = f_0( \frac{k h}{2\pi} )
  - 6 f_1( \frac{k h}{2\pi} )
  + 12 f_2( \frac{k h}{2\pi} )
  - 8 f_3( \frac{k h}{2\pi} ) .
\end{equation}
This bound is typically somewhat less than the bound associated with
second order finite differences with the same parameters.
It is straightforward to follow the recipe of
subsection~\ref{subsec:choose_semidiscrete}. For this scheme we have
again (\ref{eq:timesteps_per_period_approx_bound_newscheme}), and it
can be considered a discretization of (\ref{eq:modified_PDE}) rather
than of (\ref{eq:damped_wave_equation}).

\section{Analysis}
\label{sec:analysis}

Here theoretical results will be obtained concerning the stability of
the time integrators defined in
subsection~\ref{subsec:time-discretization}, the solutions of these
time-integration schemes and and the convergence of
the approximate solution operators defined in
subsection~\ref{subsec:time-domain_approx_sol}.

The consequences of the stability analysis for the choice of $\Delta
t$ were already discussed in subsection~\ref{subsec:choose_semidiscrete}.

\subsection{Stability}
\label{subsec:stability}
To study the conditions under which (\ref{eq:cd_integrator})
and (\ref{eq:acd_integrator})
are stable, we study the growth of the solutions to the homogeneous recursion
\begin{equation} \label{eq:homogeneous_recursion_cd}
  ( 1 + \tfrac{1}{2} L ) u_{n+1} + (- 2+K) u_n + (I - \tfrac{1}{2}  L) u_{n-1} = 0 .
\end{equation}
Throughout we will assume that $K$ and $L$ are symmetric matrices and
that $L$ is positive semidefinite.  The basic idea is to show that the
following energy functions stays bounded
\begin{equation} \label{eq:define_Ecd}
  E_{\rm cd}(n-1/2) = \langle u_n - u_{n-1} , (4I - K)  (u_n - u_{n-1}  ) \rangle
  + \langle u_n + u_{n-1} , K (u_n + u_{n-1}  ) \rangle .
\end{equation}
Here $\langle \cdot , \cdot \rangle$ denotes the standard inner
product. If the energy function is coercive, than the solution also
stays bounded. The following theorem, states the stability conditions
that follow from such an analysis. The conditions under (\ref{it:item_ii_stability_cd})
in the theorem were already mentioned above, see equations
(\ref{eq:stability_cond_K_L}), (\ref{eq:stability_cond_cd_upper}).

\begin{theorem} \label{th:stability_cd}
  Let $u_n$, $n \in {\mathbb{Z}}$ be a solution to
  the homogeneous recursion (\ref{eq:homogeneous_recursion_cd}),
  where $K$ and $L$ are real valued symmetric matrices and $L$ is
  positive semidefinite.
  \begin{enumerate}[(i)]
  \item
    If $K$ and $4I - K$ are positive definite then $E_{\rm cd}$ is equivalent
    to a norm on ${\mathbb{R}}^{2N}$. If $L = 0$ then $E_{\rm cd}$ is conserved and
    solutions remain bounded if $t \to \pm \infty$. If $L$ is positive
    semidefinite then $\Delta E_{\rm cd}(n) \le 0$
    and solutions remain bounded if $t \to \infty$.
  \item \label{it:item_ii_stability_cd}
    If $K$ is positive semidefinite and $4I - K$ is positive definite
    and $K$ has one or more zero eigenvalues then $E_{\rm cd}$ is not equivalent to a
    norm. If $L = 0$ then $E_{\rm cd}$ is
    conserved and solutions grow at most linearly if $t \to \pm \infty$.
    If $L$ is positive semidefinite, then $\Delta E_{\rm cd}(n) \le 0$ and
    solutions grow at most linearly if $t \to \infty$.
  \end{enumerate}
\end{theorem}

\begin{proof}
  To prove the result, a bound is derived for
  \begin{equation}
    \Delta E_{\rm cd}(n) = E_{\rm cd}(n+1/2) - E_{\rm cd}(n-1/2) .
  \end{equation}
  Using that $K$ is symmetric in combination with basis rules for
  standard products results in
\begin{equation}
  \begin{aligned}
  \Delta E_{\rm cd}(n) =  {}&
  \langle u_{n+1} - 2 u_n + u_{n-1} , (4I - K) ( u_{n+1}
  - u_{n-1} ) \rangle
  \\
  {}& +
  \langle u_{n+1} - u_{n-1} , K ( u_{n+1} + 2 u_n
  + u_{n-1} ) \rangle .
\end{aligned}
\end{equation}
By using the recursion equation in two places results one obtains the estimate
\begin{equation}
  \begin{aligned}
    \Delta E_{\rm cd}(n) = {}&
  \langle -K u_n - \tfrac{1}{2} L(u_{n+1} - u_{n-1} ) , (4I - K) ( u_{n+1}
  - u_{n-1} ) \rangle
  \\
  {}& +
  \langle u_{n+1} - u_{n-1} ,
  4 K u_n + K ( -K u_n - \tfrac{1}{2}  L(u_{n+1} - u_{n-1} ) ) \rangle
  \\
  = {}&
  - 2 \langle u_{n+1} - u_{n-1} ,
  L  ( u_{n+1} - u_{n-1} ) \rangle
  \\
  \le {}& 0 ,
\end{aligned}
\end{equation}
where the last inequality holds because by assumption $L$ is positive semidefinite.
The theorem follows from this result.
\end{proof}

\medskip

To establish stability of backward differences damped leapfrog, we study
the growth of solutions to the recursion
\begin{equation} \label{eq:homogeneous_recursion_bd}
  u_{n+1} + (- 2+K + L) u_n + (I - L) u_{n-1} = 0
\end{equation}
using the energy function
\begin{equation}
  \begin{aligned}
    E_{\rm bd}(n-1/2) = {}& \langle u_n - u_{n-1} , (4I - K)  (u_n - u_{n-1}  ) \rangle
    + \langle u_n + u_{n-1} , K (u_n + u_{n-1}  ) \rangle
    \\
    {}& - 2 \langle u_n - u_{n-1}, L (u_n - u_{n-1}) \rangle .
  \end{aligned}
\end{equation}

\begin{theorem}
  Let $u_n$, $n \in {\mathbb{Z}}$ be a solution to
  the homogeneous recursion (\ref{eq:homogeneous_recursion_bd}),
  where $K$ and $L$ are real valued symmetric matrices.
  \begin{enumerate}[(i)]
  \item
    If $K$ and $4I - K - 2L$ are positive definite, then $E_{\rm bd}$ is equivalent
    to a norm on ${\mathbb{R}}^{2N}$. If $L = 0$ then $E_{\rm bd}$ is conserved and
    solutions remain bounded if $t \to \pm \infty$. If $L$ is positive
    semidefinite then $\Delta E_{\rm bd}(n) \le 0$
    and solutions remain bounded if $t \to \infty$.
  \item
    If instead $K$ is positive semidefinite with one or more zero eigenvalues
    and $4I - K - 2L$ is positive definite, then $E_{\rm bd}$ is not equivalent to a
    norm. If $L = 0$ then $E_{\rm bd}$ is
    conserved and solutions grow at most linearly if $t \to \pm \infty$,
    If $L$ is positive semidefinite, then $\Delta E_{\rm bd}(n) \le 0$ and
    solutions grow at most linearly if $t \to \infty$.
  \end{enumerate}
\end{theorem}

\begin{proof}
  Similarly as in (\ref{eq:define_Ecd}) define
  \begin{equation} \label{eq:define_Ebd}
    \Delta E_{\rm bd}(n) =  E_{\rm bd}(n+1/2) - E_{\rm bd}(n-1/2) .
  \end{equation}
  Using that $K$ is symmetric in combination with basis rules for
  standard products results in
  \begin{equation}
    \begin{aligned}
      \Delta E_{\rm bd}(n) =  {}&
      \langle u_{n+1} - 2 u_n + u_{n-1} , (4I - K) ( u_{n+1}
      - u_{n-1} ) \rangle
      \\
      {}& +
      \langle u_{n+1} - u_{n-1} , K ( u_{n+1} + 2 u_n
      + u_{n-1} ) \rangle
    \\
      {}& - 2 \langle u_{n+1} - u_n, L (u_{n+1} - u_n) \rangle
    + 2 \langle u_n - u_{n-1}, L (u_n - u_{n-1}) \rangle .
    \end{aligned}
  \end{equation}
By using the recursion equation in two places results one obtains the estimate
\begin{equation}
  \begin{aligned}
    \Delta E_{\rm bd}(n) = {}&
  \langle -K u_n -  L(u_n - u_{n-1} ) , (4I - K) ( u_{n+1}
  - u_{n-1} ) \rangle
  \\
  {}& +
  \langle u_{n+1} - u_{n-1} ,
  4 K u_n + K ( -K u_n - L(u_n - u_{n-1} ) ) \rangle
    \\
      {}& - 2 \langle u_{n+1} - u_n, L (u_{n+1} - u_n) \rangle
    + 2 \langle u_n - u_{n-1}, L (u_n - u_{n-1}) \rangle
  \\
  = {}&
  - 4 \langle u_{n+1} - u_{n-1} ,
  L  ( u_n - u_{n-1} ) \rangle
    \\
      {}& - 2 \langle u_{n+1} - u_n, L (u_{n+1} - u_n) \rangle
    + 2 \langle u_n - u_{n-1}, L (u_n - u_{n-1}) \rangle
  \\
  ={}&
  -2 \langle u_{n+1} - u_{n-1}, L ( u_{n+1} - u_{n-1} ) \rangle
  \\
  \le {}& 0 ,
\end{aligned}
\end{equation}
where the last inequality holds because by assumption $L$ is positive semidefinite.
The theorem follows from this result.
\end{proof}

\subsection{Eigenvalue analysis}

An analysis based on the eigenvalues of $K$ gives insight
in the question whether these conditions are necessary.
Suppose $v$ is an eigenvector of $K$ with eigenvalue
$k$ and simultaneously an eigenvector of $L$ with
eigenvalue $\ell$. One can think of the situation that
$L = \ell I$. We start by looking for solutions to equation
(\ref{eq:homogeneous_recursion_cd})
of the form $u_n = \Lambda^n v$.
The constant $\Lambda$ must then satisfy
\begin{equation}  \label{eq:scalar_recursion_cd}
  (1 + \ell /2) \Lambda^2 + (-2+k) \Lambda + (1- \ell /2) = 0 .
\end{equation}
The solution to this equation are
\begin{equation} \label{eq:Lambda_pm_stability}
  \Lambda_\pm = \frac{-1+k/2}{1+\ell/2}
  \pm \frac{1}{1+\ell/2} \sqrt{ -k + k^2/4 + \ell^2/4}
\end{equation}
We consider the growth of solutions, first if $\ell = 0$ (no damping).
Then if $0 < k < 4$ both solutions are complex with
$| \Lambda_\pm| = 1$. If $k=0$ then
(\ref{eq:homogeneous_recursion_cd}) has solutions
$u_n = C_1 v + C_2 n v$.  If $k<0$ or $k>4$ at least one of
$\Lambda_\pm$ has modulus larger than 1.  Next consider the case
$\ell > 0$ (damping present).  If $0 < k < 4$ then $\Lambda_\pm$ can
be complex (small damping) or real (large damping), with
$| \Lambda_\pm | < 1$ in both cases.  If $k=0$ or $k=4$ one
$\Lambda _\pm$ is equal to 1 and the other smaller than 1. If $k<0$ or
$k>4$, exponentially growing solutions exist.
For stability, $K$ and $4I - K$ need to be positive semidefinite, but
we will require $4I - K$ to be strictly positive definite (to avoid
the solution with $\Lambda  = -1$ while $k=4$, $\ell>0$).

\medskip

Next we study the undamped solutions in case the eigenvectors of $L$
are not those of $K$. In this case, the recursion
(\ref{eq:homogeneous_recursion_cd}) can be written in first order
form. Let
\begin{equation}
  y_n = \begin{bmatrix} u_{n-1} \\ u_n \end{bmatrix} , \qquad
  g_n = \begin{bmatrix} 0 \\ f_n \end{bmatrix} ,
\end{equation}
and
\begin{equation}
  \Xi =
  \begin{bmatrix}
    0 & I \\
    - (I+L)^{-1} (I-L)
    & (I+L)^{-1} (2I-K)
  \end{bmatrix}
\end{equation}
then the integrators $I_{\rm cd}$ and $I_{\rm acd}$ from
definitions~\ref{def:Icd_Ibd} and~\ref{def:Iacd_Iabd}
 can be written as
\begin{equation} \label{eq:system_dim2N}
  y_{n+1} = \Xi y_n + g_n .
\end{equation}
The properties of the solutions are directly related to the
eigenvalues, (generalized) eigenvectors and the Jordan blocks of
$\Xi$. We next list some properties of $\Xi$ that relate to
eigenvalues $\xi$ of $\Xi$ with $| \xi | = 1$. These correspond to
undamped solutions. The basic idea is to
show that there are not more of these than there should be.

\begin{theorem}
  Assume (\ref{eq:stability_cond_K_L}),
  (\ref{eq:stability_cond_cd_upper}). The following properties hold
  for the eigenvalues, (generalized) eigenvectors and Jordan blocks of
  $\Xi$:
\begin{enumerate}[(i)]
\item \label{it:i}
  all eigenvectors are of the form $[u,\xi u]$, with $u \in {\mathbb{C}}^N$;
\item \label{it:ii}
  there are no eigenvalues $\xi = -1$;
\item \label{it:iii}
  for $|\xi| = 1$, $\xi \neq 1$ there are no Jordan blocks of size $>1$;
\item \label{it:iv}
  for $|\xi| = 1$, $\xi \neq 1$ the $u$ in $[u,\xi u]$ is in $\ker
L$, and is an eigenvector of $K$ and
$\xi$ follows from (\ref{eq:Lambda_pm_stability});
\item \label{it:v}
  for $\xi = 1$, there are no Jordan blocks of size $>2$,
an eigenvector is of the form
$v = [ u, u]$ with $u \in \ker K$ and if $w$ is a generalized
eigenvector such that $\Xi w = w
+ v$, $v$ an eigenvector, then $v = [u,u]$ with $u$ also in $\ker L$, and
$w = [ w_1 , w_1 + u ]$ with $w_1 \in \ker K$.
\end{enumerate}
\end{theorem}

\begin{proof}
  Claim (\ref{it:i}) is obvious.

  Regarding (\ref{it:ii}), the equations $\Xi [ u ,\xi u]^T = \xi [ u ,\xi u]^T$
  with $\xi = -1$ gives $(4I - K ) u = 0$ which is impossible because
  $4I - K$ is positive definite by assumption.

  Claim (\ref{it:iii}) follows from the energy growth equations.
  Solution to second order recursion would be of the form
  $u_n = \xi^n w + n \xi^{n-1} v$, and the dominant term in the energy in the
  limit $n \to \infty$ would be
  $\langle n (1-\xi) v , (4I - K) n (1- \xi) v \rangle$. This term would
  be unbounded as $n \to \infty$ which is impossible, hence any Jordan
  block of size $>1$ must be associated with an eigenvalue $\xi = 1$.

  In the situation of claim (\ref{it:iv}), solutions are of the form $u_n = \xi^n
  v$ with $v \in {\mathbb{C}}^N$. Energy must be conserved hence
  $\langle u_{n+1} - u_{n-1} , L (u_{n+1} - u_{n-1}) \rangle = 0$ hence
  hence $(\xi - \xi^{-1} ) u \in \ker L$, hence $u \in \ker L$. But then
  it follows that $K u = ( - \xi^{-1} + 2 - \xi) u$, so that $u$ is a
  simultaneous eigenvector of $K$, and $L$.

  The first part of part (\ref{it:v}) follows from the bound on the
  energy. Jordan blocks of size $>2$ would lead to solutions of the form
  \begin{equation}
    u_n = w^{(2)} + n w^{(1)} + n (n-1) v ,
  \end{equation}
  The dominant contribution to the term
  $\langle u_n-u_{n-1} , (4I-K) (u_n-u_{n-1}) \rangle$ would be
  $\langle  2(n-1) v , (4I-K) 2 (n-1) v \rangle$, which would be
  unbounded, which is not possible.
  The second part of (\ref{it:v}) follows from the form of $\Xi$.
  If eigenvector is $V = [ v_1, v_2]^T$, $v_j \in {\mathbb{C}}^N$, then $v_1 = v_2$ and
  \begin{equation}
    -(1-L) v_1 + (2I-K) v_1 = (1+L) v_1
  \end{equation}
  hence $K v_1 = 0$.  The third part of (\ref{it:v}) follows from the
  form of $\Xi$.  Working out the first line of equation $\Xi w = w + v$
  yields $w_2 = w_1 + u$, working out the second line gives $L u = 0$.
\end{proof}

\subsection{Convergence of the approximate solution operators $S_T^{\rm
  acd}$ and $S_T^{\rm abd}$}

In this subsection we will write $u( n \Delta t)$ instead of
$u_n$ and $f( n \Delta t)$ instead of $f_n$. Instead of the two-sided
sequences $(u_n)$ and $(f_n)$ we write $u(t)$, $f(t)$, with $t \in \Delta t {\mathbb{Z}}$. The Fourier
transform from functions on $\Delta t {\mathbb{Z}}$ to functions on
${\mathbb{T}}_{2\pi/\Delta t}$ is as defined in
appendix~\ref{app:Fourier_analysis_grid}.

Associated with the integrator $I_{\rm acd}$ is a second order
difference operator $C_{\rm acd}$. We define this operator as follows
\begin{equation} \label{eq:operator_Cacd}
  C_{\rm acd}(u)(t) = \frac{\alpha}{\Delta t^2}
\left[ ( 1 + \tfrac{1}{2} L ) u(t+\Delta t) + (- 2+K) u(t) + (I -
  \tfrac{1}{2} L) u(t-\Delta t) \right]
\end{equation}
where $K, L$ are as defined in (\ref{eq:K_L_frequency_adapted_central_differences}).
Solutions to the time-integration satisfy
\begin{equation} \label{eq:difference_equation_acd}
  C_{\rm acd}(u)(t) = f(t) .
\end{equation}
The left-hand side can be see as the convolution of $u$ with a kernel
$\Gamma_{\rm acd}$ that is a matrix valued function on $\Delta t {\mathbb{Z}}$.
(The explicit expression for
$\Gamma_{\rm acd}$ is easily derived from (\ref{eq:operator_Cacd}).)
In the Fourier domain the convolution becomes a multiplication and we have
\begin{equation}
  \widehat{\Gamma}_{\rm acd}(\tau) \widehat{u}(\tau) =
  \widehat{f}(\tau) .
\end{equation}
From
proposition~\ref{prop:reproduce_symbol_at_omega} it follows that
\begin{equation} \label{eq:symbol_equality_property_acd}
  \widehat{\Gamma}_{\rm acd}(\omega) = - \omega^2 I + i \omega B + A .
\end{equation}
In addition we know that $\widehat{\Gamma}$ is $C^\infty$, because
$\Gamma$ is nonzero only for $t \in \{ -\Delta t, 0 , \Delta t \}$.

By $\Phi_{\rm acd}(t)$, $t \in \Delta t {\mathbb{Z}}$ we will denote the causal
Green's function associated with
(\ref{eq:difference_equation_acd}). It is a function
$\Delta t {\mathbb{Z}} \to {\mathbb{C}}^{N \times N}$
defined as the solution to
\begin{equation}
  C_{\rm acd} \Phi_{\rm acd} (t) = I \delta(t) ,
  \qquad \text{and} \qquad
  \Phi_{\rm acd}(t) = 0 \text{ if } t \le  0 ,
\end{equation}
where $\delta(t)$ is as defined in
appendix~\ref{app:Fourier_analysis_grid}.
Assuming that the conditions in
(\ref{eq:stability_cond_K_L}), (\ref{eq:stability_cond_cd_upper})
hold, $\Phi_{\rm acd}$ grows at most polynomially if $t \to \infty$
(by Theorem~\ref{th:stability_cd}), and
its Fourier transform $\widehat{\Phi}_{\rm acd}$ is well defined as a matrix valued
distribution on the circle ${\mathbb{T}}_{2\pi/\Delta t}$.

We will first show that the approximate solution operator
$S_T^{\rm acd}$ can be expressed as a weighted mean of
$\widehat{\Phi}(\tau)$ around $\tau = \omega$.  For this purpose we
define a helper function $\phi_\epsilon$, depending on $\chi$, as
follows
\begin{equation}
  \psi_T(t) = \chi(1 - \left| \frac{t}{2 \pi \omega^{-1} T} \right| ) , \qquad \text{and} \qquad
  \phi_\epsilon(\tau) = \frac{1}{2 \pi} \mathcal{F} \psi_{1/\epsilon}(\tau) .
\end{equation}
Note that $\phi_\epsilon$ is an approximation to the identity. We have
the following result

\begin{theorem} \label{th:convolution_formula_S_T_acd}
Assume $\phi_\epsilon$ and $\Phi_{\rm acd}$ are as just defined, then
\begin{equation} \label{eq:convolved_Fourier_transform_acd}
  S_{T}^{\rm acd} =  \phi_{1/T} \ast \widehat{\Phi}_{\rm acd}(\omega) .
\end{equation}
\end{theorem}

In words, $s_T^{\rm acd}$ equals the convolution product of
$\phi_{1/T}$ and $\widehat{\Phi}_{\rm acd}$, evaluated at $\tau =
\omega$. The convolution product is taken on the torus
${\mathbb{T}}_{2\pi/\Delta t}$.

\begin{proof}
Let $\tilde{T} = 2\pi \omega^{-1} T$.
By linearity, the time-domain approximate solution $S_T(F)$ satisfies
\begin{equation}
  S_T^{\rm acd}(F) = e^{- i \omega \tilde{T}}
  \Delta t \sum_{t \in \Delta t {\mathbb{Z}}, 0 \le t < \tilde{T}}  \Phi(\tilde{T}-t)
  \chi(t/\tilde{T}) e^{i \omega t} F
\end{equation}
Inserting the definition of $\psi_{T}$ this becomes
\begin{equation}
  S_T^{\rm acd}(F) = \left[
  \Delta t \sum_{t \in \Delta t {\mathbb{Z}}}  \Phi(\tilde{T}-t) \psi_T(\tilde{T}-t)
  \chi(t/\tilde{T}) e^{-i \omega (\tilde{T}-t)}  \right]  F
= \left[
  \mathcal{F} ( \psi_T \Phi )(\omega)
  \right] F .
\end{equation}
The equality (\ref{eq:convolved_Fourier_transform_acd}) follows because
the Fourier transform of the product is $\frac{1}{2\pi}$ times
the convolution product of the Fourier transforms,
see (\ref{eq:FourierTransformFG}).
\end{proof}

Concerning the convergence of the time-domain approximate solutions to the
correct solutions we have the following result.
\begin{theorem} \label{th:convergence_acd}
  If
  \begin{equation} \label{eq:H_non-singular}
    \text{$H = -\omega^2 I + i \omega B + A$ is non-singular}
  \end{equation}
  and $A,B,\omega$ and $\Delta t$ are such that $K$, $L$ defined in
  (\ref{eq:K_L_frequency_adapted_central_differences}) satisfy
  the stability conditions (\ref{eq:stability_cond_K_L}) and (\ref{eq:stability_cond_cd_upper})
  then
  \begin{equation}
    \lim_{T \to \infty} S_T^{\rm acd} = H^{-1} .
  \end{equation}
\end{theorem}

\begin{proof}
  The Fourier transform $\widehat{\Phi}_{\rm acd}$ satisfies
  \begin{equation}
    \widehat{\Gamma}_{\rm acd}(\tau) \widehat{\Phi}_{\rm acd}(\tau) = I .
  \end{equation}
  Note that
  $\widehat{\Gamma}_{\rm acd}(\tau)$ is a $C^\infty$ function and
  $\widehat{\Phi}_{\rm acd}$ is a distribution and that this equality
  is in the sense of distributions.

  By the assumption (\ref{eq:H_non-singular}) and the property
  (\ref{eq:symbol_equality_property_acd}), the matrix valued function
  $\widehat{\Gamma}_{\rm acd}(\tau)$ is non-singular on a neighborhood
  of $\tau = \omega$. It follows that
  $\widehat{\Phi}_{\rm acd}(\tau)$ is continuous on a neighborhood of
  $\tau = \omega$ and that
  \begin{equation}
    \widehat{\Phi}_{\rm acd}(\omega)
    = \widehat{\Gamma}_{\rm acd}(\omega)^{-1}
    = H^{-1} .
  \end{equation}
  Because $\phi_\epsilon$ is an approximation to the identity and
  $\widehat{\Phi}(\tau)$ is continuous at $\tau = \omega$,
  it follows from theorem~\ref{th:convolution_formula_S_T_acd} and
  standard results in integration theory and distribution theory that
  \begin{equation}
    \lim_{T \to \infty} S_T^{\rm acd} = \widehat{\Phi}(\omega)     = H^{-1} . \qedhere
  \end{equation}
\end{proof}

\medskip

A similar analysis can be performed for the approximate solution
operator $S_T^{\rm abd}$. The equivalent of (\ref{eq:operator_Cacd})
is
\begin{equation} \label{eq:operator_Cabd}
  C_{\rm abd}(u)(t) = \frac{\alpha}{\Delta t^2}
  \left[ u(t+\Delta t) + (- 2+K+L) u(t)
    + (I - L) u(t-\Delta t) \right]
\end{equation}
where $K,L$ are as defined in (\ref{eq:K_L_frequency_adapted_backward_differences}).
A causal Green's function $\Phi_{\rm abd}$ is defined satisfying
\begin{equation}
  C_{\rm abd} \Phi_{\rm abd} (t) = I \delta(t) ,
  \qquad \text{and} \qquad
  \Phi_{\rm abd}(t) = 0 \text{ if } t \le  0 ,
\end{equation}
The following theorems are proved in the same ways
as theorems~\ref{th:convolution_formula_S_T_acd} and
\ref{th:convergence_acd}

\begin{theorem} \label{th:convolution_formula_S_T_abd}
Assume $\phi_\epsilon$ and $\Phi_{\rm abd}$ are as just defined, then
\begin{equation} \label{eq:convolved_Fourier_transform_abd}
  S_{T}^{\rm abd} =  \phi_{1/T} \ast \widehat{\Phi}_{\rm abd}(\omega) .
\end{equation}
\end{theorem}

\begin{theorem} \label{th:convergence_abd}
  If
  \begin{equation} \label{eq:H_non-singular_abd}
    \text{$H = -\omega^2 I + i \omega B + A$ is non-singular}
  \end{equation}
  and $A,B,\omega$ and $\Delta t$ are such that $K$, $L$ defined in
  (\ref{eq:K_L_frequency_adapted_backward_differences}) satisfy
  the stability conditions (\ref{eq:stability_cond_K_L}) and (\ref{eq:stability_cond_bd_upper})
  then
  \begin{equation}
    \lim_{T \to \infty} S_T^{\rm abd} = H^{-1} .
  \end{equation}
\end{theorem}

\section{Numerical examples}
\label{sec:numerical_examples}

In this section we study some examples. The examples all involve an
optimized finite-difference discretization described in
subsection~\ref{subsec:optimized-fd} in two and three dimensions.

To study the 2-D method, a simple implementation in Julia was made.
For the 3-D case a mixed Julia/C implementation was made. In this case
the time stepping was done in C, the computation of coefficients and
the GMRES iteration were done in Julia. In the C implementation, for
each grid point and timestep, 14 real coefficients had to be read, of
which one also had to be written back to. The C implementation used AVX
extensions but was otherwise straightforward. In most computations
double precision numbers were used. The exception was the
time-stepping performed in C. We found that the time-domain preconditioner
could also be run in single precision, since the GMRES solver would
automatically take care of rounding errors in subsequent iterations.

\subsection{2-D examples}

We considered three 2-D examples: A constant velocity model and two
piecewise constant models defined on the unit square. The first of the
two piecewise constant models allows for resonances in the circular
region.  On the exterior of the model damping layers of thickness 32
gridpoints were added to simulate an unbounded domain.
The simulations in these examples
were done using a minimum of 6 points per wavelength. In each case
this resulted in 8 timesteps per period.
See
figure~\ref{fig:model_and_sols2D} for non-constant velocity models and
some solutions in these models.

The first objective was to the study the convergence of the
approximate solution $S_T^{\rm acd} F$ to $H^{-1} F$. For this we took models of size
$320\times320$ and $640 \times 640$ (excluding damping layers), and
let the right-hand side be a point source.
The ``taper'' parameter $\rho$ was varied, we considered values in
$\{ 0, 0.25, 0.5, 0.75 \}$. The value $0$ was not exactly zero, in
this case the time-harmonic right-hand side grew from 0 to 1 over one
period.
The convergence (relative difference between approximate and true solutions)
is given in Figure~\ref{fig:large_time_convergence}. Note the
differing axes in the Figure~\ref{fig:large_time_convergence}(b).
The main conclusions from these figures are that the convergence in
the resonant model is relatively poor, and that the choice of $\rho = 0$ (no
tapering) leads to poor convergence. The detailed amount of tapering
is less important. The lower value ($\rho = 0.25$) appears to
perform best unless extremely small relative errors are required.
\newcommand{\addsubfigurelabel}[3]{%
  \begin{minipage}{#2}
    \begin{center}
      #1\\
      #3
    \end{center}
  \end{minipage}}

\begin{figure}
  \begin{center}
    \addsubfigurelabel{(a)}{6cm}{%
      \includegraphics[width=6cm]{ vel2.pdf}}
    \addsubfigurelabel{(b)}{6cm}{
     \includegraphics[width=6cm]{ vel3.pdf}}
    \addsubfigurelabel{(c)}{6cm}{%
      \includegraphics[width=6cm]{ sol_vel2_1.pdf}}
    \addsubfigurelabel{(d)}{6cm}{
     \includegraphics[width=6cm]{ sol_vel3_1.pdf}}
  \end{center}
  \caption{Two velocity models (a,b) and solutions in these models
    (c,d).
  A relatively small example of $240 \times 240$ was used so that the
  wavefield oscillations are still visible.
}\label{fig:model_and_sols2D}
\end{figure}

\begin{figure}
  \begin{center}
    \addsubfigurelabel{(a)}{65mm}{%
    \includegraphics[width=65mm]{ largeTime_velocity1_loglog.pdf}} \hspace*{5mm}
    \addsubfigurelabel{(b)}{65mm}{%
    \includegraphics[width=65mm]{ largeTime_velocity2_loglog.pdf}}\\
    \addsubfigurelabel{(c)}{65mm}{%
    \includegraphics[width=65mm]{ largeTime_velocity3_loglog.pdf}}\\
  \end{center}
  \caption{Convergence behavior of large time approximate
    solutions. Note the different scales in figure (b).}
  \label{fig:large_time_convergence}
\end{figure}

GMRES convergence for the three velocity models of size 640 gridpoints
(excluding damping layers) are given in
Figure~\ref{fig:GMRES_convergence}. In each case the number of periods
of the preconditioner was varied.
In Figures~\ref{fig:GMRES_convergence}(a-c) the taper parameter
was 0.25
and both the error indicator from GMRES and the true error are
plotted. In Figure~\ref{fig:GMRES_convergence}(d) the GMRES error was
plotted and the taper parameter was varied.
On the x axis, the number of iterations time the number of periods per
iteration was displayed. This is roughly, but not quite proportional
to the cost, the main difference being that one more preconditioner
application is needed for the right-hand side in the preconditioned
system.
In velocity models 1 and 3, there was no speedup compared to the
direct application of the time-domain solver. Neither was the method
much slower. However, in the resonant velocity model (velocity model
2), the performance using GMRES was substantially better than when the
time-domain solver was used directly.

In the examples about 50 to 100 periods for the preconditioner worked
best. This is approximately the size of the example in wavelenghts.
\begin{figure}
  \begin{center}
    \addsubfigurelabel{(a)}{65mm}{%
    \includegraphics[width=65mm]{ krylov_velocity1_semilogy.pdf}} \hspace*{5mm}
    \addsubfigurelabel{(b)}{65mm}{%
    \includegraphics[width=65mm]{ krylov_velocity2_semilogy.pdf}}\\
    \addsubfigurelabel{(c)}{65mm}{%
    \includegraphics[width=65mm]{ krylov_velocity3_semilogy.pdf}} \hspace*{5mm}
    \addsubfigurelabel{(d)}{65mm}{%
    \includegraphics[width=65mm]{ krylov_velocity2_semilogy_tap.pdf}}
  \end{center}
\caption{Preconditioned GMRES convergence for different models and different
  parameter choices} \label{fig:GMRES_convergence}
\end{figure}

\subsection{Examples in 3-D}

We also studied the performance of the algorithm in three
dimensions. In this case, the velocity model was the SEG/EAGE Salt
Model. This is a standard, non-resonant synthetic model in the geophysics community.
In this case our goal was to obtain a first estimate of the
computational cost of the method. A minimum of 6 gridpoints per
wavelength was used. The Julia/C code was run on a 2019 MacBook Pro with
a 2.6 GHz 6-core Intel i7 processor and 16 GB main memory. A GMRES
error reduction with a factor  $10^{-5}$ was required.
The results are in Table~\ref{tab:seg_eage_salt_example}.

\begin{table} \label{tab:seg_eage_salt_example}
\begin{center}
  \begin{tabular}{| l | c | c | c |}
    \hline
    physical frequency & 2.5 Hz & 4 Hz & 6 Hz\\
    problem size   & $200 \times 200 \times 106$
                                & $280 \times 280 \times 130$
    & $392 \times 392 \times 165$\\
    degrees of freedom         & 4.24 e6 & 1.02 e7 & 2.54 e7\\
    time steps/period   & 8 & 8 & 8 \\
    periods/iteration     & 25 & 40 & 60 \\
    taper parameter          & 0.25 & 0.25 & 0.25\\
    iterations                 & 6 & 6 & 6 \\
    computation time & 45 s & 166 s & 601 s\\ \hline
  \end{tabular}
\end{center}
\caption{Computational results for the SEG/EAGE Salt model}
\end{table}

\section{Generalizations}

As outlined in the introduction, the new time-domain preconditioned
GMRES solver is based on four innovative ideas: A modified
semi-discrete system derived from the discrete Helmholtz equation
(equation (\ref{eq:general_semi-discrete_equation})), a
new time-integration scheme, introduction of a smooth window in the
computation of limiting-amplitude solutions and the use of the
time-domain approximate solution operator as a preconditioner.

Some of these ideas can be used separately. The new semi-discrete
system and the new integration method could be used with the exact
controllability or WaveHoltz methods.
In this section we discuss the
application of Krylov subspace method in combination with a different
semi-discrete equation and time integrator.

We assume the
semi-discrete system is still of the form
(\ref{eq:general_semi-discrete_equation}). Denoting
\begin{equation}
  v = \begin{bmatrix} u \\ \pdpd{u}{t}
  \end{bmatrix} ,
  \qquad \;\;
  g = \begin{bmatrix} 0 \\ f
  \end{bmatrix} ,
  \qquad \;\;
  \Xi = \begin{bmatrix}
    0 & I \\
    -A & -B
  \end{bmatrix}
\end{equation}
this can be written in first order form
\begin{equation} \label{eq:general_semi-discrete_first_order}
  \pdpd{v}{t} = \Xi v + g .
\end{equation}
The first order system can be discretized for example using a standard
fourth order Runge-Kutta method, which yields
\begin{equation}
  \begin{aligned}
  v_{n+1} = {}& (I + h \Xi + \frac{h^2}{2} \Xi^2 + \frac{h^3}{6} \Xi^3
  + \frac{h^4}{24} \Xi^4) v_n
  \\
  {}& + \frac{h}{24} ( 4 + 4 h \Xi + 2 h^2 \Xi^2 + h^3 \Xi^3 ) g_n
  + \frac{h}{12} (  8 + 4 h \Xi  + h^2 \Xi^2)  g_{n+1/2}
  + \frac{h}{6} g_{n+1}
\end{aligned}
\end{equation}
Requiring that $v_n = e^{i n \omega \Delta t} V$ and
$g_n = e^{i n \omega \Delta t} G$ leads to the following equation for
$V$ and $G$
\begin{equation}
  \begin{aligned}
    {}& \bigg[   e^{i\omega \Delta t} I
  -  I - h \Xi - \frac{h^2}{2} \Xi^2 - \frac{h^3}{6} \Xi^3
  - \frac{h^4}{24} \Xi^4 \bigg]  V =
  \\
  {}& \bigg[ \frac{h}{24} ( 4 I + 4 h \Xi + 2 h^2 \Xi^2 + h^3 \Xi^3 )
  + e^{i \omega \Delta /2} \frac{h}{12} (  8 I + 4 h \Xi  + h^2 \Xi^2)
  + e^{i \omega \Delta t} \frac{h}{6} I \bigg] G  .
\end{aligned}
\end{equation}
This is the exact discrete Helmholtz equation satisfied by the
periodic solutions of the discrete-time-and-space system.
Based on the experience in this paper, this system could be solved by
GMRES preconditioned by a time-domain preconditioner
defined as in subsection~\ref{subsec:time-domain_approx_sol} but applied to
(\ref{eq:general_semi-discrete_first_order}) solved by RK4.
For different Runge-Kutta integrators, different polynomials in front
of $V$ and $G$ result, but the principle is the same.

\section{Concluding remarks}

In this paper we constructed a time-domain preconditioner for certain
indefinite linear systems, and studied its properties and behavior
analytically and with numerical examples.  It should be emphasized
that the method does not compute in the physical time domain. We call
it a time-domain method, because of the similarity with classical
time-domain methods for time-harmonic waves.

For the practical application it would be useful to have further
examples and a comparison with alternatives, and more insight into how
the method behaves for different choices of $H$, e.g.\ for different
discretizations. We will make a few brief remarks in this direction.

First it is clear that the method requires relatively little
memory, compared to
alternatives that use LU (or LDL${}^T$) decompositions, such as
domain-decomposition and direct methods, see
\cite{stolk2017improved,poulson2013parallel} and  references in
\cite{gander2019class}.
In the 3-D implementation here, with GMRES with restart $n_{\rm r}
= 10$, most memory was used for the approximately
$n_{\rm r} + 2$ complex double
precision vectors needed for GMRES.
This was in part because the timestepping was done in single precision, and
using real fields. Memory use could be somewhat reduced by setting
$n_{\rm r} = 5$ or using a different iterative method like BiCGSTAB.

For finite element discretizations the cost in general will be
different. These methods often have stricter CFL bounds, and there may
be more nonzero matrix elements, depending on the order of the finite
elements, which leads to larger cost. In case of unstructured meshes
computations in general are less efficient, because there is more
indirect memory access and cache misses. On the other hand, the size
of the grid cells can be adapted to the local velocity which means the
number of grid cells can be smaller.

It is difficult to compare performance of different methods, as
methods are run on different computer systems, with different examples
and implementations are optimized to different degrees. We refer to
\cite{appelo2019waveholtz,CalandraEtAl2013,liu2018solving,
  poulson2013parallel,stolk2017improved,taus2019sweeps,WangDeHoopXia2011}
for some alternative methods. We believe that computation times of the
method as outlined are modest, and there is potential for further
improvements, by using GPUs or by reorganizing the code to better
make use of cached data.  It remains challenging to get the most out
of modern computer hardware in time-domain finite-difference and
finite-element simulations. We hope that recent developments in this
area, cf.\ \cite{kukreja2016devito}, will lead to further improvements.

\bibliographystyle{abbrv}