# Superstep wavefield propagation

Tamas Nemeth1, Kurt Nihei1, Alex Loddoch1, Anusha Sekar1, Ken Bube2,
John Washbourne1, Luke Decker1, Sam Kaplan1, Chunling Wu1,
Andrey Shabelansky1, Milad Bader1, Ovidiu Cristea1 and Ziyi Yin3

(1 Chevron Technical Center, a division of Chevron U.S.A. Inc.
2 University of Washington
3 Georgia Institute of Technology
May 31, 2024)

###### Abstract

This paper describes how to propagate wavefields for arbitrary numbers of
traditional time steps in a single step, called a superstep. We show how to
construct operators that accomplish this task for finite-difference time
domain schemes, including temporal first-order schemes in isotropic,
anisotropic and elastic media, as well as temporal second-order schemes for
acoustic media. This task is achieved by implementing a computational tradeoff
differing from traditional single step wavefield propagators by precomputing
propagator matrices for each model location for $k$ timesteps (a superstep)
and using these propagator matrices to advance the wavefield $k$ time steps at
once. This tradeoff separates the physics of the propagator matrix computation
from the computer science of wavefield propagation and allows each discipline
to provide their optimal modular solutions.

##  1 Introduction

Traditional wavefield propagation using finite-difference time domain (FD)
schemes advances the wavefield one timestep at a time by performing a matrix-
vector multiplication of the current wavefields with the combination of
material properties (such as velocity or density) in a pattern described by
the FD kernel. The accuracy, stability and dispersion of such FD schemes have
been studied extensively (Hairer et al. [1], Chapman [2], Durran [3], Bube et
al. [4], Bube et al. [5], Bube et al. [6], Washbourne et al. [7]) and are well
understood. Today, most modeling and imaging applications are built on linear
or linearized FD schemes from acoustic forward modeling to vertical and tilted
transverse isotropy (VTI, TTI) and elastic modeling.

Traditional recursive time steps with such a FD scheme operationally
correspond to a repeated matrix-vector multiplication, although the matrices
and vectors may not be explicitly formed. Such schemes can be viewed as a
series of matrices (one for each time step) applied to a vector representing
the current wavefield state. It is widely accepted, especially for large
matrices, that computationally it is better to carry out this computation
sequentially from “right” to “left”; start by multiplying the initial
wavefield vector on the right with the matrix for the first time step
resulting in a new vector; applying the next timestep matrix to the new
vector, until all matrix-vector multiplications are executed.

Performing operations in this manner reduces computational cost, as these
traditional schemes map well onto current hardware (CPU, GPU) where the
wavefield propagation problem is considered as a mostly compute-bound
computational problem. This established process performs well when most of
these jobs are run on high-end hardware nodes and storage systems in the data
centers owned by the organizations themselves. With the availability of broad
classes of compute nodes and storage options in the cloud, hardware choices
have become more diverse and, correspondingly, the associated cost structure
has become more complex but visible. For wavefield propagation using the cloud
infrastructure, the overall runtime nowadays mostly depends on the
sophistication of managing the storage and IO of the wavefields and earth
model material parameters both in the memory of compute nodes and online
storage, and on the agility of the scheduling systems for virtual machine
compute nodes. These developments necessitate the re-examination of the
traditional computational paradigm for wavefield propagation.

This paper presents numerical formulations for wavefield propagation where the
tradeoff for compute and storage is different from the traditional formulation
and the algorithmic components are more general. These formulations enable us
to recast the process of wavefield propagation as two distinct steps of
“Precompute” and “Compute”. During Precompute, the common reusable components,
the propagator matrices (PM) are calculated and stored. This Precompute stage
is separate from the subsequent Compute step and can therefore be executed on
separate systems possibly with a different hardware architecture. During the
actual Compute stage, the precomputed PMs are used and re-used to achieve
wavefield propagation.

There were several earlier efforts to use Precompute for part of the wavefield
propagation. Notable is the checkpointing scheme of Symes [8] where wavefields
are systematically saved to match the forward and adjoint wavefields. Kole [9]
used a matrix exponential approach to derive unconditionally stable algorithms
to approximate time evolution (numerical integration of the matrix exponential
power series) by operator splitting and one-step Chebyshev polynomials. Araujo
and Pestana [10] also used matrix exponential approaches to model the time
evolution of the first-order linear acoustic and elastic wave equations. In a
different domain, stencil computations were approximated with Gaussian fitting
and then these Gaussians were used for further propagation (Ahmad et al.
[11]). Jun [12] and Vamaraju et al. [13] formulated seismic imaging as a
repeated application of system matrices. Vdovina et al. [14], Korostyshevskaya
et al. [15] and Gibson et al. [16] developed an operator upscaling method to
capture the finer subgrid effects on a coarser domain without updating the
wavefields on a finer scale.

We have found that the optimal Precompute is to calculate the superstepping
propagator matrices associated with each gridpoint in the spatial
discretization of the wave propagation domain as these objects are the most
re-usable and can be applied to wavefields at any timesteps. The corresponding
Compute is the integration of the PMs with the wavefields at each grid point
and the output is the updated wavefield.

The key observation to superstepping is to re-arrange the FD schemes in such a
way that current and previous wavefields are related to future and subsequent
wavefields via a square system matrix. This way the square matrix can be
raised to the $k$th power corresponding to superstepping of $k$ traditional
time steps. When the initial wavefield is a delta function for each spatial
location, we can precompute the Green’s functions at each location to the
$k$th traditional time steps. These precomputed PMs become impulse response
filters at these locations. Then we apply these filters to the actual
wavefields to advance the wavefield by $k$ steps (a superstep).

Similar ideas have appeared in different contexts. Seismic interferometry uses
Green’s functions to reconstruct wavefields from other wavefields (Wapenaar et
al. [17], Schuster [18], Wapenaar et al. [19]). Superstepping can be
considered as a special case of seismic interferometry. Similarly, Marchenko
methods (Slob et al. [20], Lomas and Curtis [21], Wapenaar [22]) reconstruct
subsurface Green’s functions and use them for imaging. Wapenaar [22]
especially, shows the interrelation and differences between Green’s functions,
propagator matrices and Marchenko focusing functions. Since a propagator
matrix depends only on the medium parameters between the two time levels where
defined in the medium, we will use the terminology of propagator matrices
thereon in accordance with Wapenaar [22].

Recently, many novel machine learning techniques were developed that train to
predict system states continuously from current time to up to a given future
time. These methods, including Fourier Neural Operators [23], [24] use the
inherent nonlinearity in the neural network to provide a model for inference.
Although both machine learning methods and superstepping can predict future
states of the system, superstepping is based on the linear properties of the
system and is not an approximation to the traditional time stepping schemes.
Some neural network methods provide models to some tasks used in traditional
wavefield propagation such as boundary condition prediction (Kaur et al.
[25]). In this work we show that these tasks can also be precomputed as they
are part of the linear system matrix.

As seen in later sections, the Precompute stage can be identified in
incorporating physics into the Green’s functions and the Compute stage can be
viewed as a computer science effort to match computational tasks to the given
computing resources. For the heterogeneous cloud-based compute and storage
environment, this match may be best accomplished via efficient domain-specific
languages (DSL). Earlier, Charlesworth et al. [26] developed a DSL for finite-
difference modeling across CPU, FPGA and GPU hardware. Recently the Devito
framework [27] is showing promise as a DSL for FD stencil computations. We
hope that the superstepping scheme will provide yet another DSL abstraction
with the benefit that the Compute stage involves only computer science
considerations.

The manuscript is organized as follows. In the Theory section we give a
concise description of superstepping for both first- and second-order systems
of wave propagation and show that simply an outcome of precomputing the $k$th
power of the square is the propagator matrix. In the System matrix
representations as propagator matrices section we show that the propagator
matrices are composed of PMs for each row of the matrix and this greatly
reduces the non-zero elements of these matrices. In the Superstepping
computational model section we describe the numerical aspects of both
Precompute and Compute stages. In the Synthetic examples section we
demonstrate superstepping on a Marmousi 2D example and show the treatment of
boundary reflections in the Green’s functions. The Discussion section
emphasizes one of the main results, that the Compute stage is is the stage for
computational and algorithmic implementations rather than physics. Finally, in
the Appendices section we give detailed examples of formulating the first and
second-order acoustic, VTI anisotropic and elastic wave equations in terms of
superstepping.

##  2 Theory

Many physical phenomena can be expressed via partial differential equations
(PDE) where state variables change over time and space. An important class of
these partial differential equations is shown in the next set of equations
where state variables $\mathbf{u}$ and $\mathbf{v}$ are related to each
other’s temporal changes and typically indicate the velocity and stress
relations in the medium. The state variables $\mathbf{u}$ and $\mathbf{v}$ in
the physical domain can be both multi-dimensional and multi-component
quantities.

###  2.1 First-order partial differential equations

A typical first order system of equations with 2 state variables for wave
propagation can be written in a matrix equation form as

$\displaystyle\left[\begin{array}[]{c}\mathbf{\dot{u}}\\\
\mathbf{\dot{v}}\end{array}\right]=\begin{bmatrix}0&A\\\
B&0\end{bmatrix}\left[\begin{array}[]{c}\mathbf{u}\\\
\mathbf{v}\end{array}\right]~{},$ (5)

where $\mathbf{\dot{u}}$ and $\mathbf{\dot{v}}$ are the temporal derivatives
of state variables $\mathbf{u}$ and $\mathbf{v}$, typically wavefields of
particle velocities and stress in our case. Matrices $A$ and $B$ incorporate
material coefficients and partial derivative operators such as $\nabla$ and
$\nabla\cdot$ operators. Note that $A$ and $B$ are off the main diagonal of
the matrix indicating that the spatial variations of one quantity affects the
temporal variations of the other one. This is typical for many PDEs in wave
propagation and we will provide specific examples later. There are more
complex systems of first order equations for other physical formulations and
we will discuss later as well.

Equation 5 can be discretized using the symplectic Euler method:

$\displaystyle\mathbf{u}_{n+1}$ $\displaystyle=$
$\displaystyle\mathbf{u}_{n}+A\,\mathbf{v}_{n}~{},$
$\displaystyle\mathbf{v}_{n+1}$ $\displaystyle=$
$\displaystyle\mathbf{v}_{n}+B\,\mathbf{u}_{n+1}~{},$ (6)

where $\mathbf{\dot{u}}=(\mathbf{u}_{n+1}-\mathbf{u}_{n})/\Delta t$,
$\mathbf{\dot{v}}=(\mathbf{v}_{n+1}-\mathbf{v}_{n})/\Delta t$, and $\Delta t$
is absorbed into $A$ and $B$ for simplicity. The index $n$ denotes the
discretized time step. Sometimes the $\mathbf{u}$ variable is staggered
compared to the $\mathbf{v}$ variable and $\mathbf{\dot{u}}$ is instead equal
to $(\mathbf{u}_{n+1/2}-\mathbf{u}_{n-1/2})/\Delta t$. The staggered nature of
time steps for the variables is implied when needed. Equation 2.1 can be
expressed in a matrix form as

$\displaystyle\left[\begin{array}[]{c}\mathbf{u}_{n+1}\\\
\mathbf{v}_{n+1}\end{array}\right]=\begin{bmatrix}I&0\\\
B&I\end{bmatrix}\begin{bmatrix}I&A\\\
0&I\end{bmatrix}\left[\begin{array}[]{c}\mathbf{u}_{n}\\\
\mathbf{v}_{n}\end{array}\right]=\begin{bmatrix}I&A\\\
B&I+BA\end{bmatrix}\left[\begin{array}[]{c}\mathbf{u}_{n}\\\
\mathbf{v}_{n}\end{array}\right]~{},$ (13)

where the Lie-Trotter split-step update matrices $\begin{bmatrix}I&0\\\
B&I\end{bmatrix}$ and $\begin{bmatrix}I&A\\\ 0&I\end{bmatrix}$ describe the
update from $\mathbf{u}_{n}$ to $\mathbf{u}_{n+1}$ and from $\mathbf{v}_{n}$
to $\mathbf{v}_{n+1}$. These split-step matrices are invertible
$\begin{bmatrix}I&0\\\ B&I\end{bmatrix}^{-1}=\begin{bmatrix}I&0\\\
-B&I\end{bmatrix}$ , $\begin{bmatrix}I&A\\\
0&I\end{bmatrix}^{-1}=\begin{bmatrix}I&-A\\\ 0&I\end{bmatrix}$, and the
inverse split-step update matrices describe the update from $\mathbf{u}_{n}$
to $\mathbf{u}_{n-1}$ and from $\mathbf{v}_{n}$ to $\mathbf{v}_{n-1}$
indicating forward and backward time symmetry.

The update is a 2-step process for this system with 2 classes of variables.
The progression of the update in terms of the resulting vectors can be shown
as

$\displaystyle\left[\begin{array}[]{c}\mathbf{u}_{n}\\\
\mathbf{v}_{n}\end{array}\right]\Rightarrow\left[\begin{array}[]{c}\mathbf{u}_{n+1}\\\
\mathbf{v}_{n}\end{array}\right]\Rightarrow\left[\begin{array}[]{c}\mathbf{u}_{n+1}\\\
\mathbf{v}_{n+1}\end{array}\right]\Rightarrow\left[\begin{array}[]{c}\mathbf{u}_{n+2}\\\
\mathbf{v}_{n+1}\end{array}\right]\Rightarrow\left[\begin{array}[]{c}\mathbf{u}_{n+2}\\\
\mathbf{v}_{n+2}\end{array}\right]$ (24)

From equation 24 it is obvious that repeated applications of these matrices
constitute time propagation of these state variables. This can be expressed by
raising the propagator matrices to the $k$th power for $k$ time step
advancement, as shown in

$\displaystyle\left[\begin{array}[]{c}\mathbf{u}_{n+k}\\\
\mathbf{v}_{n+k}\end{array}\right]=\begin{bmatrix}I&A\\\
B&I+BA\end{bmatrix}^{k}\left[\begin{array}[]{c}\mathbf{u}_{n}\\\
\mathbf{v}_{n}\end{array}\right]~{}.$ (29)

Equation 29 implies that we can calculate the
$\mathbf{u}_{n+k},\mathbf{v}_{n+k}$ wavefields from the
$\mathbf{u}_{n},\mathbf{v}_{n}$ by either (a) sequentially applying the
transform matrix $k$ times as in traditional timestepping or (b) calculating
the effects of the transform matrix to the $k$th power and applying it to
$\mathbf{u}_{n}$ and $\mathbf{v}_{n}$.

###  2.2 Second-order partial differential equations

Second-order system of equations typically arise from first-order system of
equations. In this section we extend the equations for second-order systems
based on the learnings from the previous section. The corresponding second-
order system (in time) derived from equation 5 is

$\displaystyle\left[\begin{array}[]{c}\mathbf{\ddot{u}}\\\
\mathbf{\ddot{v}}\end{array}\right]=\begin{bmatrix}0&A\\\
B&0\end{bmatrix}\left[\begin{array}[]{c}\mathbf{\dot{u}}\\\
\mathbf{\dot{v}}\end{array}\right]=\begin{bmatrix}0&A\\\
B&0\end{bmatrix}^{2}\left[\begin{array}[]{c}\mathbf{u}\\\
\mathbf{v}\end{array}\right]=\begin{bmatrix}AB&0\\\
0&BA\end{bmatrix}\left[\begin{array}[]{c}\mathbf{u}\\\
\mathbf{v}\end{array}\right]~{},$ (38)

recalling that $A$ and $B$ are time invariant. Note that the $A$ and $B$
block-matrices are in the main diagonal of the system matrix, while they were
in off-diagonal position for the first-order system. Using equation 13, we can
propose the finite difference solution of equation 38 as follows:

$\displaystyle\left[\begin{array}[]{c}\mathbf{u}_{n+1}\\\
\mathbf{v}_{n+1}\end{array}\right]=\begin{bmatrix}I&A\\\
B&(I+BA)\end{bmatrix}^{2}\left[\begin{array}[]{c}\mathbf{u}_{n-1}\\\
\mathbf{v}_{n-1}\end{array}\right]~{}.$ (43)

We show next that equation 43 can be transformed to the traditional second-
order finite-difference schemes. Specifically we need to show that
$\mathbf{v}_{n+1}$ in equation 43 can be written only as a function of
$\mathbf{v}_{n}$ and $\mathbf{v}_{n-1}$ with no $\mathbf{u}_{n}$ and
$\mathbf{u}_{n-1}$ contributions. The bottom line of equation 13,

$\mathbf{v}_{n+1}=B\mathbf{u}_{n}+(I+BA)\,\mathbf{v}_{n}~{},$ (44)

can be brought to the following form by substituting the second line in
equation 2.1 to yield

$\mathbf{v}_{n+1}=2\,\mathbf{v}_{n}-\mathbf{v}_{n-1}+BA\,\mathbf{v}_{n}~{}.$
(45)

In conventional FD schemes, this equation would have been used to directly
discretize equation 38. Note that the same derivation could have been shown
for $\mathbf{u}_{n+1}$ using equations 2.1 and 13 as shown in the derivation
below

$\mathbf{u}_{n+1}=\mathbf{u}_{n}+A\mathbf{v}_{n}=\mathbf{u}_{n}+A(\mathbf{v}_{n-1}+B\mathbf{u}_{n})=\mathbf{u}_{n}+\mathbf{u}_{n}-\mathbf{u}_{n-1}+AB\,\mathbf{u}_{n}~{}.$
(46)

The second-order discretization scheme contains 3 time steps,
$\mathbf{v}_{n-1}$, $\mathbf{v}_{n}$ and $\mathbf{v}_{n+1}$. The following
system advances the time steps by one step and contains a square matrix

$\displaystyle\left[\begin{array}[]{c}\mathbf{v}_{n}\\\
\mathbf{v}_{n+1}\end{array}\right]=\begin{bmatrix}0&I\\\
-I&2I+BA\end{bmatrix}\left[\begin{array}[]{c}\mathbf{v}_{n-1}\\\
\mathbf{v}_{n}\end{array}\right]=\begin{bmatrix}0&I\\\
-I&G\end{bmatrix}\left[\begin{array}[]{c}\mathbf{v}_{n-1}\\\
\mathbf{v}_{n}\end{array}\right]~{},$ (53)

where for simplicity of notation we have designated $G=2I+BA$, where
$G=G(k=1)$ is the block matrix containing a single-step Green’s function in
each row as seen later. Squaring the system matrix advances the time steps by
two, from $\mathbf{v}_{n-1}$, $\mathbf{v}_{n}$ steps to $\mathbf{v}_{n+1}$,
$\mathbf{v}_{n+2}$ steps

$\displaystyle\left[\begin{array}[]{c}\mathbf{v}_{n+1}\\\
\mathbf{v}_{n+2}\end{array}\right]=\begin{bmatrix}0&I\\\
-I&G\end{bmatrix}^{2}\left[\begin{array}[]{c}\mathbf{v}_{n-1}\\\
\mathbf{v}_{n}\end{array}\right]=\begin{bmatrix}-I&G\\\
-G&G^{2}-I\end{bmatrix}\left[\begin{array}[]{c}\mathbf{v}_{n-1}\\\
\mathbf{v}_{n}\end{array}\right]~{},$ (60)

where the squared system matrix is symplectic (Hairer [1]). We can generalize
the system to $k$ arbitrary time steps by bringing the system matrix to the
$k$th power:

$\displaystyle\left[\begin{array}[]{c}\mathbf{v}_{n+k-1}\\\
\mathbf{v}_{n+k}\end{array}\right]=\begin{bmatrix}0&I\\\
-I&G\end{bmatrix}^{k}\left[\begin{array}[]{c}\mathbf{v}_{n-1}\\\
\mathbf{v}_{n}\end{array}\right]~{},$ (65)

Note that the structure of the second-order system matrix corresponds to the
structure of the system matrix for the first-order system. Comparing equation
60 with equation 29, the differences are: (i) in the second-order system the
main diagonal has a $-I$ component rather than $I$ as it corresponds to the
second to last time step, (ii) the off-diagonal terms are skew symmetric. As
in the case of the first-order system, we can calculate the
$\mathbf{v}_{n+k-1}$ and $\mathbf{v}_{n+k}$ steps by either (a) sequentially
applying the transform matrix k times as in traditional timestepping or (b)
calculating the effects of the transform matrix to the k-th power and applying
it to $\mathbf{v}_{n-1}$ and $\mathbf{v}_{n}$.

Similarly to the first-order systems, these transform matrices are invertible
indicating forward and backward time symmetry $\begin{bmatrix}0&I\\\
-I&G\end{bmatrix}^{-1}=\begin{bmatrix}G&-I\\\ I&0\end{bmatrix}$ .

##  3 System matrix representations as propagator matrices

In this section we relate the system matrices to physical PMs and describe the
wavefield propagation in terms of dot products applications of PMs to current
state variables (wavefields).

Denote the resulting first- and second order systems from equations 29 and 65
as

$\displaystyle\left[\begin{array}[]{c}\mathbf{u}_{n+k}\\\
\mathbf{v}_{n+k}\end{array}\right]=\begin{bmatrix}I&A\\\
B&I+BA\end{bmatrix}^{k}\left[\begin{array}[]{c}\mathbf{u}_{n}\\\
\mathbf{v}_{n}\end{array}\right]=\begin{bmatrix}S_{11}(k)&S_{12}(k)\\\
S_{21}(k)&S_{22}(k)\end{bmatrix}\left[\begin{array}[]{c}\mathbf{u}_{n}\\\
\mathbf{v}_{n}\end{array}\right]~{},$ (72)

and

$\displaystyle\left[\begin{array}[]{c}\mathbf{v}_{n+k-1}\\\
\mathbf{v}_{n+k}\end{array}\right]=\begin{bmatrix}0&I\\\
-I&G\end{bmatrix}^{k}\left[\begin{array}[]{c}\mathbf{v}_{n-1}\\\
\mathbf{v}_{n}\end{array}\right]=\begin{bmatrix}S_{11}(k)&S_{12}(k)\\\
S_{21}(k)&S_{22}(k)\end{bmatrix}\left[\begin{array}[]{c}\mathbf{v}_{n-1}\\\
\mathbf{v}_{n}\end{array}\right]~{},$ (79)

where in either case $\bm{S}(k)$ denotes the $k^{th}$ superstep propagator
matrix with the pre-calculated block matrices

$\displaystyle S(k)=\begin{bmatrix}S_{11}(k)&S_{12}(k)\\\
S_{21}(k)&S_{22}(k)\end{bmatrix}~{}.$ (80)

Each block matrix $S_{ij}(k)$ corresponds to the $i^{th}$ output state
variable vector and the $j^{th}$ input state variable vector. The size of
block matrix $S_{ij}(k)$ is $m\times q$ where $m$ is the length of the
$i^{th}$ output state variable vector and $q$ is the length of the $j^{th}$
input state variable vector. To better understand $S(k)$ we need to look at
the structure of the individual block matrix products such as

$\displaystyle\mathbf{\hat{u}}_{n+k}=S_{11}(k)\;\mathbf{u}_{n}$ (81)

where $\mathbf{\hat{u}}_{n+k}$ is the top-left product in equation 72.
Equation 81 is a discretized form of linear integral kernel operators. In wave
equations these kernel operators correspond to PMs and the functions these
kernels operate on are the time source functions. Therefore, $S_{11}(k)$ is a
matrix of PMs, each PM forming one row of the matrix and $\mathbf{u}_{n}$ is
the time source function.

A sketch of equation 81 is shown indicating contributions to the output

$\displaystyle\begin{bmatrix}\hat{u}_{0}\\\ \hat{u}_{1}\\\ \hat{u}_{2}\\\
\vdots\end{bmatrix}_{n+k}=\begin{bmatrix}[l]\leftarrow
G_{0}(k)\rightarrow\;~{}\;&\cdots\\\ \;\;\;\;\;\;\;\;\;\;\leftarrow
G_{1}(k)\rightarrow\;&\\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\leftarrow
G_{2}(k)\rightarrow&\\\
\;\;\;\;\;\;\;\;\vdots&\ddots\end{bmatrix}\begin{bmatrix}u_{0}\\\ u_{1}\\\
u_{2}\\\ \vdots\end{bmatrix}_{n}~{}.$ (82)

where $G_{i}(k)$ denotes non-zero segments of each propagator matrix that
corresponds to $k$ timesteps. Note that the dimension of the block matrix rows
indicates the physical dimension of the discretized earth models, that is
$i=N_{x}*N_{y}*N_{z}$ for a 3D earth model where $N_{x},N_{y},N_{z}$ are the
dimensions in the 3 earth model dimensions. The non-zero elements are
clustered around the main diagonal. For example, when $k=0$ (no superstep),
the matrix turns into a unit matrix and the output vector is the same as the
input vector. As superstepping size $k$ in $S(k)$ grows, the non-zero elements
in each row also grow.

The $G_{i}(k)$ non-zero segments in each column of the block matrices
constitute the propagator matrix corresponding the unit input for each element
in the input vector calculated by solving the PDE for $k$ time steps. This
requires knowing the earth model parameters only and therefore can be pre-
computed for each column and stored for subsequent use. The reciprocity rules
between input and output applies and the column based propagator matrices also
serve as propagator matrices for each row (Wapenaar [22]).

Figure 1: Sketch illustrating the propagating matrices (PM) in their physical
domain. 2 PMs $G(A)$ and $G(B)$ shown corresponding to 30 timesteps of
propagation from locations $A$ and $B$. The compact support for $G(A)$ and
$G(B)$ are shown as gray squares. The wavefield $u_{60}$ at time step 60 is
indicated in light gray shade and it is computed for the whole physical
domain.

Since the block matrix is sparse, we can represent the matrix-vector multiply
in a more compact form, keeping only the $G_{i}(k)$ non-zero segments of the
matrix that are pre-computed Green’s functions and the spatially corresponding
subset of $\mathbf{u}_{n}$ denoted as $\mathbf{\tilde{u}}_{n}$, written for
each row $i$ (that corresponds to a given spatial point $x,y,z$ in the model):

$\displaystyle\hat{u}_{i}(n+k)=G_{i}(k)\;\tilde{u}_{i}(n)~{}.$ (83)

Equation 83 corresponds to the dot product calculation for each vector element
$i$. A 2D illustration of the scheme is shown in Figure 1 where 2 PMs $G(A)$
and $G(B)$ shown corresponding to 30 timesteps of propagation from locations
$A$ and $B$. The compact support for $G(A)$ and $G(B)$ are shown as gray
squares. The wavefield $\mathbf{u}(t=60)$ at time step 60 is indicated in
light gray shade and it is computed for the whole physical domain.

Figure 2 shows the dot product calculation between the PM (top left) and the
corresponding wavefield segment (bottom left) resulting in the element-wise
products (bottom right). Then, these element-wise products are all summed
together to generate the resulting scalar of the dot product. Note that we
used color plots in the top row of Figure 2 to illustrate the extent of the
computations: the yellow regions are beyond where the clearly formed wavefront
is touched computationally, yet it is near-zero value and corresponds to
physically acausal area.

Figure 2: Illustration of the dot product calculation between a propagating
matrix and the corresponding wavefield. Top left: PM after 30 time steps. Top
right: Wavefield after 300 timesteps in the physical 2D domain. Bottom left:
Wavefield segment corresponding to the physical area of the PM. Bottom right:
Element-wise multiplication of the PM and the wavefield segment.

The PMs $G_{i}(k)$ for individual $i$ element are computed independently.
There is overlap in $\tilde{u}_{i}(n)$ for different $i$ values (locations).

From this we can describe the computation of of equations 72 and 79 as

Algorithm 1 Superstep Wavefield Propagation for Acoustic Wave Equation

SuperstepAcoustic
$k,n_{0},N_{t},N_{mod},\mathbf{u}(n_{0}),\mathbf{v}(n_{0}),G(k)$, $H(k)$,
$K(k)$, $L(k)\,$$n=n_{0}:k:N_{t}$$\triangleright$ iterate in superstep time
$i=1:N_{mod}$ $\triangleright$ iterate in model location

1:$G_{i}(k),H_{i}(k),K_{i}(k),L_{i}(k)=\textsc{GetPropMatrix}\left(i\right)$
$\triangleright$ Propagator matrices for the i-th location

2:$\tilde{u}_{i}(n)$,
$\hat{u}_{i}(n)=\textsc{CutWavefieldSegment}\left(\mathbf{u}(n)\right)$
$\triangleright$ Cut matching wavefield segments

3:$\tilde{v}_{i}(n)$,
$\hat{v}_{i}(n)=\textsc{CutWavefieldSegment}\left(\mathbf{v}(n)\right)$

4:$u_{i}(n+k)=G_{i}(k)\cdot\tilde{u}_{i}(n)+H_{i}(k)\cdot\tilde{v}_{i}(n)$
$\triangleright$ update $u_{i}$

5:$v_{i}(n+k)=K_{i}(k)\cdot\hat{u}_{i}(n)+L_{i}(k)\cdot\hat{v}_{i}(n)$
$\triangleright$ update $v_{i}$

6:return $\mathbf{u},\mathbf{v}$ $\triangleright$ return computed wavefield
components

where $N_{t}=T/(\Delta t)$ is the number of time steps for $T$ recording time,
$\Delta t$ is the temporal sampling rate, $k$ is the superstep size, $n_{0}$
is the initial time sample level for superstepping,
$N_{mod}=N_{x}*N_{y}*N_{z}$ is the dimension for 3D-sampled models,
$G_{i}(k),H_{i}(k),K_{i}(k),L_{i}(k)$ are propagator matrices at location $i$
corresponding to input-output pairings of $u_{i}\rightarrow
u_{i},v_{i}\rightarrow u_{i},u_{i}\rightarrow v_{i},v_{i}\rightarrow v_{i}$,
respectively,
$\tilde{u}_{i}(n),\tilde{v}_{i}(n),\hat{u}_{i}(n),\hat{v}_{i}(n)$ are the part
of the $\mathbf{u}(n)$ and $\mathbf{v}(n)$ state variables (wavefields) that
correspond to the compact size of the PMs above. Function GetPropMatrix
returns all the PMs for a given model location, and function
CutWavefieldSegment isolates the wavefield segments from the overall wavefield
that matches the spatial extent of their corresponding PM. The updated
variables $u_{i}(n+k)$ and $v_{i}(n+k)$ at time step $n+k$ are the sum of the
dot products of PMs and the corresponding state variables at the $n$-th time
step.

We can generalize the number of state variables for more than 2. While
$\mathbf{u}$ and $\mathbf{v}$ as described in equation 5 is valid, these
vectors can be multi-component themselves. For example, $\mathbf{u}$ can be
single-component if it represents a scalar variable such as pressure, or
3-component vector $\mathbf{u}=(\mathbf{u}_{x},\mathbf{u}_{y},\mathbf{u}_{z})$
if it represents a gradient, or a 6-component vector
$\mathbf{u}=(\mathbf{u}_{xx},\mathbf{u}_{yy},\mathbf{u}_{zz},\mathbf{u}_{xy},\mathbf{u}_{xz},\mathbf{u}_{yz})$
if it represents the components of more complex differential operators. The
number of components can be different for $\mathbf{u}$ and $\mathbf{v}$ as
seen in the examples shown in the Appendices. The full wavefield $\mathbf{w}$
is a combination of wavefields $\mathbf{u}$ and $\mathbf{v}$ and has $N_{c}$
components where $N_{c}$ is the sum of the components in $\mathbf{u}$ and
$\mathbf{v}$. At any location $i$ at any time step $n$ there are $N_{c}$ state
variable values and we can represent them by denoting them as $w_{i\alpha}$.
For scalar-component $\mathbf{u}$ and $\mathbf{v}$ this corresponds to
$N_{c}=2$ with $\mathbf{w}_{0}=\mathbf{u}$ and $\mathbf{w}_{1}=\mathbf{v}$.
Since the state variables are $N_{c}$-component for both the input and the
output, the corresponding system matrix is the size of $N_{c}\times N_{c}$
block matrices. Therefore we will need to calculate and use $N_{c}\times
N_{c}$ propagator matrices for each location. These PMs at a given location
corresponds to the various combinations of input and output components. We
will denote these PMs as $G_{i\alpha\beta}(k)$ where $i$ is the index of the
physical location, $\alpha$ is the index of the output component, $\beta$ is
the index of the input component and $k$ is the superstep size. These
generalizations can be written as:

Algorithm 2 Superstep Wavefield Propagation for General Wave Equations

SuperstepGeneral
$k,n_{0},N_{t},N_{mod},N_{c},\mathbf{w}(n_{0}),G(k)$$n=n_{0}:k:N_{t}$$\triangleright$
iterate in superstep time $i=1:N_{mod}$ $\triangleright$ iterate in model
location $\alpha=1:N_{c}$ $\triangleright$ iterate for output wavefield
components $\beta=1:N_{c}$ $\triangleright$ iterate for input wavefield
components

1:$G_{i\alpha\beta}(k)=\textsc{GetPropMatrix}\left(i\right)$ $\triangleright$
Propagator matrices for the i-th location

2:$\tilde{w}_{i\alpha\beta}(n)=\textsc{CutWavefieldSegment}\left(\mathbf{w}(n)\right)$
$\triangleright$ Cut matching wavefield segments

3:$w_{i\alpha}(n+k)=\sum_{\beta=1}^{N_{c}}G_{i\alpha\beta}(k)\cdot\tilde{w}_{i\alpha\beta}(n)$
$\triangleright$ update for $w_{i\alpha}$

4:return $\mathbf{w}$ $\triangleright$ return computed wavefield components

We can also generalize Algorithm 2 to multiple survey records (whole survey)
by defining tensor $w_{i\alpha j}$ where each column $w_{i\alpha}$ is a single
record as above and where index $j$ indicates the index of the record:

Algorithm 3 Multi-survey Superstep Wavefield Propagation for General Wave
Equation

SuperstepGeneralMultiSurvey
$k,n_{0},N_{t},N_{mod},N_{c},N_{s},\mathbf{w}(n_{0}),G(k)$$n=n_{0}:k:N_{t}$$\triangleright$
iterate in superstep time $i=1:N_{mod}$ $\triangleright$ iterate in model
location $\alpha=1:N_{c}$ $\triangleright$ iterate for output wavefield
components $j=1:N_{s}$ $\triangleright$ iterate in surveys $\beta=1:N_{c}$
$\triangleright$ iterate for input wavefield components

1:$G_{i\alpha\beta}(k)=\textsc{GetPropMatrix}\left(i\right)$ $\triangleright$
Propagator matrices for the i-th location

2:$\tilde{w}_{i\alpha\beta
j}(n)=\textsc{CutWavefieldSegment}\left(\mathbf{w}(n)\right)$ $\triangleright$
Cut matching wavefield segments

3:$w_{i\alpha
j}(n+k)=\sum_{\beta=1}^{N_{c}}G_{i\alpha\beta}(k)\cdot\tilde{w}_{i\alpha\beta
j}(n)$ $\triangleright$ update for $w_{i\alpha j}$

4:return $\mathbf{w}$ $\triangleright$ return computed wavefield components

where $N_{s}$ indicates the total number of surveys. Note that there is no
dependence on PM $G_{i\alpha\beta}(k)$ from a particular record $j$. As
readily seen above, the fundamental task, the unit of computation is the dot
product between a PM and a state variable wavefield component of the same
size. The other tasks are to get the precomputed PMs and to select the
appropriate wavefield components that correspond to those PMs. Thus we can
build a computational graph based on these tasks to account for the total
computation used during the simulation of the survey. The optimized
computational graph determines the execution order of the nested loop
elements.

##  4 Superstepping computational model

Superstepping creates two fundamentally different tasks to propagate
wavefields: (a) Precompute - creating PMs for a given model and (b) Compute -
using these PMs to advance the wavefields. In this section we examine the
computational consequences of these very different tasks.

###  4.1 Precompute: Propagator matrix generation

In equation 82, the system matrix consists of rows of PMs (flattened for
multidimensional PM) and the input and output vectors represent the current
and updated wavefields. Since the material (earth model) parameters are all
contained in the system matrix, therefore the PMs for each spatial location
represent all the physics applied in wavefield propagation, including the
physical model encoded in the partial differential equations, all the
numerical issues of physics such as dispersion, boundary conditions. In
addition, the distribution of the PMs in space and the number of PMs at each
spatial location can also be considered information about the physics.

Since the system matrices in equations 72 and 79 are fixed, the system is
linear and time invariant. As a consequence, the PMs can be precomputed and
applied to all measurement time ranges and all shot records. This presents an
opportunity to generate high quality PMs where all the desired physics is
incorporated. The PMs can also be efficiently generated in a non-overlapping
parallel way for typical model sizes.

One way to accelerate the computation of PMs is to reduce their size to their
physical domain of influence. As seen on Figure 3a, the computational domain
of influence is shown in the yellow diamond with the physical wavefront
already formed in the middle of the diamond. The corner areas in blue are
untouched by computation after these time steps. We can reduce the size of the
PMs by trimming around the wavefront while also preserving numerical
precision. This can be thought of bringing two different parts of physics into
the PMs before discarding physics: (a) computing the PM coefficients based on
physics and (b) determining the causal regions based on both math and physics,
further tightening the compact support of PMs.

Figure 3: Comparison between the numerical and physical domains of influence.
The same 30-step propagator matrix is plotted with ”Accent” color scale in (a)
and with gray in (b). The yellow diamond in (a) is the computational domain
while the formed ring of wave in (b) is the physical domain.

Another way to improve the computed PMs is to incorporate the physics of
handling boundary reflections into them. Figure 4 shows a PM example with
boundary reflections still in it (4a), boundary reflections removed from it
(4b) and their difference (4c). In the Numerical examples section below we
demonstrate the effect with and without boundary reflections.

Figure 4: A propagator matrix near a reflecting boundary (a) without
compensation for the boundary layer; (b) compensated for the boundary layer
and (c) their difference.

We can reduce dispersion in the PMs by employing high fidelity numerical
schemes during their computation. Since PMs are shared in different time
scales and among all survey records, the optimization efforts provide
significant value for production runs. Although in our examples we ran 2nd
order spatial numerical schemes, in production we would replace them with
higher spatial order schemes.

PMs represent the filter that moves the incoming wavefields $k$ timesteps
forward essential for computing gradients for full waveform inversion or
migration images in reverse time migration. These same PMs can also be used to
move the incoming wavefield $k$ timesteps backward. In practice this
corresponds to taking the inverse of matrix in equation 80 and we illustrate
this in Appendix E.

PMs can be further improved by generating representations to reduce grid
effects or improving the wavefield continuity. These representations can be
accomplished by either traditional signal processing or machine learning.
While we do not show examples of these representations in this paper, it is
well-understood that higher quality PMs generate better wavefields.

PMs can also be used for approximations of the full wave equations, such as
rays, dynamic rays or other simpler geometrical optics approximations as
smooth gradients with second order derivatives.

###  4.2 Compute: wavefield propagation

Wavefield propagation in superstepping is agnostic to physics as physics has
already been parameterized and incorporated into PMs. Therefore one can
consider wavefield propagation a computer science task that uses parameterized
PMs. From this perspective, the implementations need to be general to
efficiently map to a broad range of compute resources and at the same time
allow for efficient representations of PMs and their interactions with
wavefields.

Once a simulation task has been initiated on a set of compute resources, a
computational graph can be created that connects PMs, source locations and
receiver locations and the simulation time ranges with wavefields. This graph
represents the overall computational effort to connect these components and it
can be mapped to (distributed) compute resources programmatically using the
most efficient computer science methods.

The generality of mapping is facilitated by the simplification of tasks during
superstepping: (a) keeping PMs in memory awaiting wavefields; (b) select and
crop the wavefields to the size of PMs and bring them to the corresponding
memory and (c) compute the dot product between the PM and the cropped
wavefield. These operations capture the available utilities on a compute node,
namely, memory, interconnect and compute resources. Maximizing all these
resources is thus a computational optimization goal.

Since all of the PMs are available at the wavefield propagation stage, the
goal is to keep all of them in memory for fast availability. For large
wavefield propagation tasks, this can be on the order of thousands to tens of
thousands of compute nodes for a non-optimized representation. For example,
consider a typical CPU node with over a 100 GB or $10^{11}$ bytes of memory.
If we perform wavefield propagation on a large $1000^{3}$ gridded model, we
would have $10^{9}$ grid locations. In each grid location there would be 3 to
45 Green’s PMs where, 3 is for acoustic and 45 for full elastic models of wave
propagation. This comes from the physical representation of system matrices as
2x2 block matrices for acoustic and 9x9 block matrices for elastic equations.
Assume the superstep to be $k=30$ and each PM is of the size of up to $30^{3}$
or 2.7$\times 10^{4}$ bytes. Therefore the total size of PMs in a location is
in the order of up to $10^{6}$ bytes. To store all these PMs in memory takes
$10^{15}$ bytes. For this example we would need tens of thousands ($\approx
10^{4}$) nodes to keep all the PMs in memory. While this is a very large
computational resource for memory, in practice we can optimize the PM
representation significantly to reduce this to hundreds to thousands of nodes
that is the typical size of the task using a traditional finite-difference
implementation. Note that keeping the Green’s functions in the memory of many
nodes sets up a large distributed computation. This is facilitated by that the
location (index $i$) is the second outermost loop in Algorithm 3.

The wavefields used for propagation (2 for acoustic and 9 for elastic) are
also distributed among the nodes. Similarly to the traditional
implementations, some halo exchanges of the wavefields happen. As each node
can host hundreds to hundreds of thousands of locations, these locations can
share the same wavefields.

The computational task of the dot product between the PM and the cropped
wavefield is less intensive compared to the task of managing PMs and
wavefields. There are several ways to optimize this dot product and to
facilitate the above goals. First, develop memory-efficient representations of
PMs that improves the dot product calculation. All PMs can be processed up-
front before the dot product computation. The coefficients can be flattened
and compressed using run-length encoding and kept in memory in a compressed
form. The dot product calculation can be executed while the PM coefficients
are in the run-length encoded form. This allows more PMs kept in the memory of
a single node and it also reduces the number of halo exchanges for the
wavefields.

The PMs can be also decomposed by bit order and the dot product executed at
bit order level taking into account the different statistics of the higher
bits and the lower bits. The multiplication at bit level is just an addition
if 1 and no action if 0. At the high-end level one can train a machine
learning model and use the trained model to infer the PM for the dot product
computation. Second, we can optimize the dot product computation by utilizing
that many PMs are kept in a single node’s memory. From the previous example,
45 PMs are kept in memory per location (elastic). These can be combined into a
single simplified representation where their dot products are computed
jointly. In the prior example, $10^{5}$ locations will fit in single node’s
memory that is, a size of $45^{3}$ sub-cubes can be jointly optimized by
forming a reduced representation by truncated SVD or other mathematical
techniques.

It is worth mentioning that Algorithm 3 prescribes a distributed computing
model that is unusual from traditional finite-difference modeling perspective.
It favors parallelization by earth model locations and not by individual shot
record locations. Note that the earth model locations (loops indexed by $i$)
are outer to the shot records (loops indexed by $j$) in Algorithm 3.

##  5 Synthetic examples

In this section we illustrate some aspects of superstepping on a subset of the
well-known 2D Marmousi model [28]. First, we demonstrate that Green’s
functions can incorporate the physics of handling boundary reflections without
any specific boundary reflection treatment during superstepping. Later we show
an example of superstepping on Marmousi model compared to traditional
wavefield propagation.

Figure 5 shows a wavefield propagating in a homogeneous medium with 3500 m/s
velocity. The model size is 301x301 with 10m spatial sampling. The source is
at location (1000,1000), denoted by a white dot in the bottom left figure.
Several PMs are shown over the model indicating the relative size of PMs after
30 regular time steps. The top left figure shows the wavefield after 420
traditional single-step time stepping with no boundary condition applied.
Strong boundary reflections are observed in the top and left of the model. The
top right figure shows the wavefield after superstepping using the PMs in
Figure 3 where the boundary conditions were addressed during the PM
calculation. The PMs were calculated for $k=30$ superstep size. A total of 13
supersteps were applied to the initial wavefield at 30 traditional steps. No
boundary conditions and boundary zones were defined and used during
superstepping. The bottom right figure shows the difference between the
wavefield in the top left and right figures. As seen, the difference
corresponds to the boundary reflections.

The superstepping procedure can also be interpreted in terms of Huygens’
principle. All points of the wavefield $\mathbf{u}(n)$ after $n$ time steps
act as point sources that scatter waves for $k$ time steps. The sum of the
scattered waves result in the full wavefield at time step $n+k$. Including
both $\mathbf{u}(n)$ and $\mathbf{u}(n-1)$ in the superstepping procedure
provides the temporal direction so that the full wavefield propagates forward
in time. Equivalently, the superstepping procedure can be interpreted as a
gathering operation. In this case, for a homogeneous medium with pressure
field, the PMs $\mathbf{G}(k)$ after $k$ time steps are hollow spheres. The
zero lag correlation between $\mathbf{G}(k)$ and $\mathbf{u}(n)$ (vector dot
product) is collected at the PM source locations to yield wavefield
$\mathbf{u}(n+k)$, the only nonzero contribution from the superstepping
procedure corresponds to PMs tangent to $\mathbf{u}(n)$. Thus, the wavefield
is propagated by supersteps of size $k$. For a virtual boundary separating the
domain of interest from an infinite domain, PMs can be truncated by the
boundary either partially or totally if their center falls outside the domain
of interest. As a result, the wavefield portion that falls outside the domain
of interest will also be truncated without having any effect on the one
remaining inside. Considering the wavefield inside the domain of interest,
truncating the PMs by the virtual boundary is equivalent to propagating the
wavefield in the unbounded medium (without the boundary). Therefore, PM
truncation effectively simulates Sommerfeld radiating condition, e.g., the
wavefield only exists the domain (outward propagation) without reentry. For
more complex media with the possibility of back-scattering from the external
region, the effective internal regions needs to be enlarged until the
backscattering can be neglected.

Figure 5: Wavefield after 420 traditional time steps. Top left: Wavefield
after 420 steps calculated by traditional method; Top right: Wavefield after
420 steps calculated by superstep method. Initial wavefield at 30 steps was
advanced by 30 steps at once using propagator matrices; Bottom left:
Homogeneous velocity model overlaid by propagator matrices at proportional
size and the source location (white dot); Bottom right: Difference between the
1-step and superstep methods at the time scale of the wave field.

The superstepping scheme was applied in a similar manner to a 300x300 grid
section of the Marmousi model shown in the bottom left figure in Figure 6. The
source location is at the middle of the model at (1500m, 1500m). The top left
figure shows the wavefield after 300 traditional time steps and the top right
shows the corresponding superstep-propagated wavefield after 9 30-traditional-
step supersteps were applied to the initial wavefield. The bottom left figure
shows the difference between the wavefields in the top left and right figures.
As expected, the superstep results closely matched the results achieved using
traditional time stepping. Figure 7 shows the wavefields after 450 traditional
time steps or 14 superstep after the initial 30 traditional time steps. No
boundary layers were used and no boundary reflection compensation included in
the PMs either. As a result, both the traditional and the superstep-generated
wavefields contains significant boundary reflections and the boundary
reflections are identical. This result further confirms that superstepping
does not contain additional physics that is not already included in the PMs.

Figure 6: Wavefield after 300 traditional time steps. Top left: Wave field
after 300 steps calculated by traditional method; Top right: Wavefield after
300 steps calculated by superstep method. Initial wavefield at 30 steps was
advanced by 30 steps at once using propagator matrices; Bottom left: Velocity
model (301x301 subset of Marmousi); Bottom right: Difference between the
1-step and superstep methods at the time scale of the wave field.  Figure 7:
Wavefield after 450 traditional time steps. Sub-figures are same as in Figure
5. Please note that the boundary edge reflections are accurately reproduced by
the superstep method.

##  6 Discussion

We have developed a superstepping scheme for linear or linearized partial
differential equations where the effect of the repeated application of the
system matrix can be precomputed in an efficient manner. This result has
several important consequences.

The superstep formulation is a general scheme that applies to all all
discretized PDEs that are linear w.r.t. the solution. For example, Algorithm 3
is applicable to acoustic or elastic wave equation, with or without density
terms. It is also applicable to propagate wavefields for other, more complex
discretized equations. The difference between the various equations manifests
in the content of the PMs and the number of those PMs for a given location. In
the Appendices we gave examples for first and second-order acoustic schemes
and VTI and elastic first-order formulations.

The superstep formulation clearly separates tasks that are associated with the
physics and task associated with wavefield propagation. This was our original
goal when we initiated this effort. The first task can be precomputed and the
resulting PMs used in a later wavefield propagation. We have shown that all
physics can be accounted for during the generation of the PMs and these PMs
can be further enhanced. We have shown that the second, the wavefield
propagation task is well-suited for computer science optimization and compute
resource mapping. Specifically, the three steps in wavefield propagation: (a)
keeping PMs in memory awaiting wavefields; (b) selecting and cropping the
wavefields to the size of the PMs and bring them to the corresponding memory
and (c) compute the dot product between the PM and the cropped wavefield, are
fully utilizing the resources on compute nodes and allow computational
optimization without the need to know the underlying physics.

Propagator matrices are computed for each grid location in the earth models
and the locations are manifested in the outer loops in Algorithm 3. This
results in a distributed computing scheme where the main distribution factor
is the placement of sufficient number of PMs in a node’s memory. Efficient
memory utilization is the driving force for computational optimization and it
underscores that most large-scale wavefield propagation job are bandwidth
(memory or IO) limited.

It may be futile to advance wavefields via superstepping at many traditional
timesteps at once if in many cases we would need to interpolate back to all
time step for imaging (or adjoint-state-based gradient computation) purposes.
We have developed imaging schemes using superstepping that avoid these
interpolations and the wavefield propagation and imaging can proceed with
supersteps.

Since the wavefield propagation is heavily influenced by the computer science
optimization, the future success of this method fully depends on strong
domain-specific language implementations.

##  7 Acknowledgements

We are thankful for the Chevron Technical Center for the permission to publish
this paper. Ziyi Yin contributed to this work during an internship at Chevron
Technical Center.

foo.bib

## References

  * [1] E. Hairer, G. Wanner, Ch. Lubich, Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations, Springer, 2006.
  * [2] C.H. Chapman, Fundamentals of seismic wave propagation: Cambridge University Press, 2004.
  * [3] D.R. Durran, Numerical methods for wave equations in geophysical fluid dynamics: Springer, 1999.
  * [4] K.P. Bube, T. Nemeth, J.P. Stefani, R. Ergas, W. Liu, K.T. Nihei, and L. Zhang, On the instability in second-order systems for acoustic VTI and TTI media, Geophysics 77, 2012, T171-T186, https://doi.org/10.1190/geo2011-0250.1
  * [5] K.P. Bube, T. Nemeth, J.P. Stefani, W. Liu, K.T. Nihei, R. Ergas and L. Zhang, First-order systems for elastic and acoustic variable-tilt TI media, Geophysics 77, 2012, T157-T170, https://doi.org/10.1190/geo2011-0249.1
  * [6] K.P. Bube, J. Washbourne, R. Ergas, and T. Nemeth, Self-adjoint, energy-conserving second-order pseudoacoustic systems for VTI and TTI media for reverse time migration and full-waveform inversion, SEG International Exposition and 86th Annual Meeting, 2016, 1110-1114, https://doi.org/10.1190/segam2016-13878451.1
  * [7] J. Washbourne, K.P. Bube, R. Ergas, and T. Nemeth, System and method of estimating anisotropy properties of geological formations using a self-adjoint pseudoartistic wave propagator, 2019, U.S. Patent 10,215,869
  * [8] W.W. Symes, Reverse time migration with optimal checkpointing, Geophysics 72, 2007, SM213-SM221, https://doi.org/10.1190/1.2742686
  * [9] J.S. Kole, Solving seismic wave propagation in elastic media using the matrix exponential approach, Wave Motion 38, 2003, 38, 279-293, https://doi.org/10.1016/S0165-2125(03)00065-9
  * [10] E.S. Araujo and R.C. Pestana, Time evolution of the first-order linear acoustic/elastic wave equation using Lie product formula and Taylor expansion, Geophysical Prospecting 69, 2021, 70-84
  * [11] Z. Ahmad, R. Chowdhury, R. Das, P. Ganapathi, A. Gregory and Y. Zhu, Brief Announcement: Faster Stencil Computations using Gaussian Approximations, Proceedings of the 34th ACM Symposium on Parallelism in Algorithms and Architectures, 2022, 291-293, https://doi.org/10.1145/3490148.3538558
  * [12] J. Ji, An exact adjoint operation pair in time extrapolation and its application in least-squares reverse-time migration, Geophysics 74, 2009, 1SO-Z92, https://doi.org/10.1190/1.3173894
  * [13] J. Vamaraju and M.K. Sen, Unsupervised physics-based neural networks for seismic migration, Interpretation 7, 2019, SE189-SE200, https://doi.org/10.1190/INT-2018-0230.1
  * [14] T. Vdovina, S. Minkoff and O. Korostyshevskaya, Operator upscaling for the acoustic wave equation, Multiscale Modeling and Simulation, 4, 2005, 1305-1338, https://doi.org/10.1137/050622146.
  * [15] O. Korostyshevskaya, and S. Minkoff, A matrix analysis of operator based upscaling for the wave equation, SIAM Journal on Numerical Analysis, 44, 2006, 586-612, https://doi.org/10.1137/050625369.
  * [16] R.L. Gibson, K. Gao, E. Chung, and Y. Efendiev, Multiscale modeling of acoustic wave propagation in 2D media, Geophysics 79, 2014, T61-T75. https://doi.org/10.1190/geo2012-0208.1
  * [17] K. Wapenaar, and J. Fokkema, Greens function representations for seismic interferometry: Geophysics 71, 2006, SI33-SI46, https://doi.org/10.1190/1.2213955.
  * [18] G.T. Schuster, Seismic Interferometry, Cambridge Univ. Press, 2009.
  * [19] K. Wapenaar, J. Thorbecke, and J. van der Neut, A single-sided homogeneous Green’s function representation for holographic imaging, inverse scattering, time-reversal acoustics and interferometric Green’s function retrieval, Geophysical Supplements to the Monthly Notices of the Royal Astronomical Society, 205, 2016, 531-535, https://doi.org/10.1093/gji/ggw023.
  * [20] E. Slob, K.Wapenaar, F. Broggini, and R. Snieder, Seismic reflector imaging using internal multiples with Marchenko-type equations, Geophysics 79, 2014, S63-S76, https://doi.org/10.1190/geo2013-0095.1
  * [21] A. Lomas and A. Curtis, (2019), An introduction to Marchenko methods for imaging, Geophysics 84, 2019, F35-F45. https://doi.org/10.1190/geo2018-0068.1
  * [22] K. Wapenaar, Wave-field representations with Green’s functions, propagator matrices, and Marchenko-type focusing functions. The Journal of the Acoustical Society of America, 151(1), 2022, 587-608. https://doi.org/10.1121/10.0009236
  * [23] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, A. Anandkumar, Fourier neural operator for parametric partial differential equations. arXiv preprint arXiv:2010.08895, 2020.
  * [24] T.J. Grady, R. Khan, M. Louboutin, Z. Yin, P.A. Witte, R. Chandra, R.J. Hewett, F.J. Herrmann, Model-parallel Fourier neural operators as learned surrogates for large-scale parametric PDEs, Computers $\&$ Geosciences, Volume 178, 2023, https://doi.org/10.1016/j.cageo.2023.105402
  * [25] H. Kaur, S. Fomel and N. Pham, Automated hyperparameter optimization for simulating boundary conditions for acoustic and elastic wave propagation using deep learning, Geophysics 88, 2023, WA309-WA318, https://doi.org/10.1190/geo2022-0231.1
  * [26] O. Charlesworth, T. Nemeth, O. Pell, T. Johnsen, P. Akerberg, Unified application framework for finite-difference modeling, 2017, U.S. Patent 9576089 B2
  * [27] F. Luporini, M. Louboutin, M. Lange, N. Kukreja, P. Witte, J. Hückelheim, Ch. Yount, P.H.J. Kelly, F.J. Herrmann, and G.J. Gorman, Architecture and Performance of Devito, a System for Automated Stencil Computation, ACM Transactions on Mathematical Software, Volume 46, Issue 1, 2020, 1-28, https://doi.org/10.1145/3374916
  * [28] R. Versteeg, The Marmousi experience: Velocity model determination on a synthetic complex data set, The Leading Edge 13, 1994, 927-936. https://doi.org/10.1190/1.1437051

##  Appendix A Matrix exponential derivation of the first order wave equations

In this section we we show an alternative derivation of the superstepping
scheme by using a more general matrix exponential formulation of wave
propagation. The velocity-stress equation

$\displaystyle\left[\begin{array}[]{c}\mathbf{\dot{u}}\\\
\mathbf{\dot{v}}\end{array}\right]=\begin{bmatrix}0&A\\\
B&0\end{bmatrix}\left[\begin{array}[]{c}\mathbf{u}\\\
\mathbf{v}\end{array}\right]~{},$ (A.5)

can also be written as

$\displaystyle\mathbf{\dot{w}}=H\mathbf{w}~{},$ (A.6)

where $\mathbf{w}=(\mathbf{u},\mathbf{v})^{T}$ and the system matrix denoted
by $H$. The formal solution of this matrix equation is

$\displaystyle\mathbf{w}(t)=e^{tH}\mathbf{w}(0)~{}.$ (A.7)

The matrix exponential can be written in power series:

$\displaystyle
e^{tH}=\lim\limits_{k\to\infty}\left(I+\frac{tH}{k}\right)^{k}=\sum_{m=o}^{\infty}\frac{(tH)^{m}}{m!}=I+\frac{(tH)}{1!}+\frac{(tH)^{2}}{2!}+\frac{(tH)^{3}}{3!}+\dots$
(A.8)

where $I$ is a 2x2 block identity matrix. The power series describes numerical
integration from $0$ to $t$ for an arbitrary $t$ time value. Substituting the
specific matrix for $H$ and discretizing time as $t=n\Delta t$, we get

$\displaystyle e^{tH}=\begin{bmatrix}I&0\\\ 0&I\end{bmatrix}+n\Delta
t\begin{bmatrix}0&A\\\ B&0\end{bmatrix}+\frac{(n\Delta
t)^{2}}{2!}\begin{bmatrix}AB&0\\\ 0&BA\end{bmatrix}+\frac{(n\Delta
t)^{3}}{3!}\begin{bmatrix}0&ABA\\\ BAB&0\end{bmatrix}+\frac{(n\Delta
t)^{4}}{4!}\begin{bmatrix}ABAB&0\\\ 0&BABA\end{bmatrix}+\dots$ (A.9)

where $n$ is a natural number and $\Delta t$ is the time step that is chosen
based on physical considerations. Next, we limit the power series for a single
time step ($n=1$). The zeroth order approximation ($m=0$) is limiting as in
this case no update happens ($e^{tH}\approx I$). The first order approximation
($m=1$) is also limiting as neither state variable completes a round trip from
$\mathbf{u}_{0}$ to $\mathbf{u}_{1}$, or $\mathbf{v}_{0}$ to $\mathbf{v}_{1}$,
e.g., $\mathbf{u}_{0}$ does not contribute to $\mathbf{u}_{1}$,

$\displaystyle e^{\Delta tH}\approx\begin{bmatrix}1&0\\\
0&1\end{bmatrix}+\Delta t\begin{bmatrix}0&A\\\ B&0\end{bmatrix}$ (A.10)

resulting in the following partitioned explicit Euler finite-difference scheme

$\displaystyle\mathbf{u}_{n+1}$ $\displaystyle=$
$\displaystyle\mathbf{u}_{n}+A\,\mathbf{v}_{n}~{},$
$\displaystyle\mathbf{v}_{n+1}$ $\displaystyle=$
$\displaystyle\mathbf{v}_{n}+B\,\mathbf{u}_{n}~{}.$ (A.11)

It is the second order approximation ($m=2$) where such a round trip happens
first

$\displaystyle e^{\Delta tH}\approx\begin{bmatrix}1&0\\\
0&1\end{bmatrix}+\Delta t\begin{bmatrix}0&A\\\ B&0\end{bmatrix}+\frac{(\Delta
t)^{2}}{2!}\begin{bmatrix}AB&0\\\ 0&BA\end{bmatrix}$ (A.12)

but this partitioned implicit Euler method is not symplectic and not used in
practice. Instead, the most common practice is to use symplectic operator
splitting methods, such as Lie-Trotter, Strang and higher order composition
methods (Hairer [1]). If we remove one of the block matrices from the second
order term in equation A.12 and compensate for it by multiplying the second
oder term by 2, we create a partitioned explicit-implicit Euler approximation
that is symplectic and can be split by using the Lie-Trotter splitting scheme

$\displaystyle e^{\Delta tH}\approx\begin{bmatrix}I&\Delta tA\\\ \Delta
tB&I+(\Delta t)^{2}BA\end{bmatrix}=\begin{bmatrix}I&0\\\ \Delta
tB&I\end{bmatrix}\begin{bmatrix}I&\Delta tA\\\ 0&I\end{bmatrix}=e^{\Delta
tH_{B}}\circ e^{\Delta tH_{A}}~{}.$ (A.13)

This is the approximation we used in equation 2.1.

With $n$ time steps we need $n$ round trips for state variables so for an
appropriately chosen $\Delta t$ this corresponds to applying the Lie-Trotter
split-step operators chained $n$ times

$\displaystyle e^{n\Delta tH}\approx(e^{\Delta tH_{B}}\circ e^{\Delta
tH_{A}})^{n}~{}.$ (A.14)

##  Appendix B First-order system for the acoustic wave equation

In this section we derive the first-order system matrix for the acoustic wave
equation. This wave equation is expressed as based on equation 5:

$\displaystyle\left[\begin{array}[]{c}\dot{\mathbf{v}}\\\
\dot{p}\end{array}\right]=\begin{bmatrix}0&-1/\rho\;\nabla\\\
-\kappa\nabla\cdot&0\end{bmatrix}\left[\begin{array}[]{c}\mathbf{v}\\\
p\end{array}\right]~{},$ (B.5)

where $\mathbf{v}$ is the velocity vector, $p$ is pressure, $\rho$ is density
and $\kappa$ is the adiabatic compression modulus. Equation B.5 is typically
discretized the following way:

$\displaystyle\mathbf{v}_{n+1/2}$ $\displaystyle=$
$\displaystyle\mathbf{v}_{n-1/2}-\Delta t/\rho\;\nabla\,p_{n}~{},$
$\displaystyle p_{n+1}$ $\displaystyle=$ $\displaystyle p_{n}-\Delta
t\;\kappa\nabla\cdot\,\mathbf{v}_{n+1/2}~{},$ (B.6)

Note that $\mathbf{v}$ and $p$ are on staggered grid relative to each other
and the material properties are sampled on their respective grid. Expressing
the $k$ time step shifted components $\mathbf{v}_{n-1/2+k}$ and $p_{n+k}$ from
$\mathbf{v}_{n-1/2}$ and $p_{n}$ yields

$\displaystyle\left[\begin{array}[]{c}\mathbf{v}_{n-1/2+k}\\\
{p}_{n+k}\end{array}\right]=\begin{bmatrix}I&-\Delta t/\rho\;\nabla\\\ -\Delta
t\;\kappa\nabla\cdot&I+(\Delta
t)^{2}\;\kappa\nabla\cdot(1/\rho\;\nabla)\end{bmatrix}^{k}\left[\begin{array}[]{c}\mathbf{v}_{n-1/2}\\\
p_{n}\end{array}\right]~{},$ (B.11)

where the system matrix is now the superstep propagator matrix in equation 80:

$\displaystyle\left[\begin{array}[]{c}\mathbf{v}_{n-1/2+k}\\\
{p}_{n+k}\end{array}\right]=\begin{bmatrix}S_{11}(k)&S_{12}(k)\\\
S_{21}(k)&S_{22}(k)\end{bmatrix}\left[\begin{array}[]{c}\mathbf{v}_{n-1/2}\\\
p_{n}\end{array}\right]~{}.$ (B.16)

The last equation is formally equivalent to equation 72.

##  Appendix C First-order system for the elastic wave equation

In this section we derive the first-order system for the elastic wave
equation. The elastic wave equation is expressed as based on 5:

$\displaystyle\left[\begin{array}[]{c}\dot{v}_{x}\\\ \dot{v}_{y}\\\
\dot{v}_{z}\\\ \dot{\sigma}_{xx}\\\ \dot{\sigma}_{yy}\\\ \dot{\sigma}_{zz}\\\
\dot{\sigma}_{yz}\\\ \dot{\sigma}_{xz}\\\
\dot{\sigma}_{xy}\end{array}\right]=\begin{bmatrix}0&1/\rho\;\nabla_{3\times
6}\\\ C\;\nabla_{6\times 3}&0\end{bmatrix}\left[\begin{array}[]{c}{v}_{x}\\\
{v}_{y}\\\ {v}_{z}\\\ {\sigma}_{xx}\\\ {\sigma}_{yy}\\\ {\sigma}_{zz}\\\
{\sigma}_{yz}\\\ {\sigma}_{xz}\\\ {\sigma}_{xy}\end{array}\right]~{},$ (C.19)

where $\nabla_{3\times 6}$ and $\nabla_{6\times 3}$ are the matrix
differential operators

$\displaystyle\nabla_{3\times
6}=\begin{bmatrix}\partial_{x}&0&0&0&\partial_{z}&\partial_{y}\\\
0&\partial_{y}&0&\partial_{z}&0&\partial_{x}\\\
0&0&\partial_{z}&\partial_{y}&\partial_{x}&0\end{bmatrix}~{};\;\;\;\nabla_{6\times
3}=\begin{bmatrix}\partial_{x}&0&0\\\ 0&\partial_{y}&0\\\ 0&0&\partial_{z}\\\
0&\partial_{z}&\partial_{y}\\\ \partial_{z}&0&\partial_{x}\\\
\partial_{y}&\partial_{x}&0\end{bmatrix}~{};$ (C.20)

and $\mathbf{v}=[v_{x},v_{y},v_{z}]^{T}$ is the particle velocity vector,
$[\sigma_{xx},\sigma_{yy},\sigma_{zz},\sigma_{yz},\sigma_{xz},\sigma_{xy}]^{T}$
are the six independent components of the stress tensor $\bm{\sigma}$ and bold
upper case $\bf{C}$ will denote the $6\times 6$ matrix of stiffnesses in the
Voigt notation [2]. After discretization of equation C.19 and taking $k$
traditional time steps (that is, raising the system matrix to the $k$-th
power), the state vectors at the $n+k$ step are related to the state vectors
at the $n$-th steps as

$\displaystyle\left[\begin{array}[]{c}{v}_{x}\\\ {v}_{y}\\\ {v}_{z}\\\
{\sigma}_{xx}\\\ {\sigma}_{yy}\\\ {\sigma}_{zz}\\\ {\sigma}_{yz}\\\
{\sigma}_{xz}\\\ {\sigma}_{xy}\end{array}\right]_{n+k}=\begin{bmatrix}I&\Delta
t/\rho\;\nabla_{3\times 6}\\\ \Delta t\;C\;\nabla_{6\times 3}&I+(\Delta
t)^{2}C\;\nabla_{6\times 3}(1/\rho\;\nabla_{3\times
6})\end{bmatrix}^{k}\left[\begin{array}[]{c}{v}_{x}\\\ {v}_{y}\\\ {v}_{z}\\\
{\sigma}_{xx}\\\ {\sigma}_{yy}\\\ {\sigma}_{zz}\\\ {\sigma}_{yz}\\\
{\sigma}_{xz}\\\ {\sigma}_{xy}\end{array}\right]_{n}~{},$ (C.39)

where the time step value of the state vectors is denoted at the bottom tight
corner of their square bracket. Here the system matrix is now the superstep
propagator matrix in equation 80

$\displaystyle\left[\begin{array}[]{c}{v}_{x}\\\ {v}_{y}\\\ {v}_{z}\\\
{\sigma}_{xx}\\\ {\sigma}_{yy}\\\ {\sigma}_{zz}\\\ {\sigma}_{yz}\\\
{\sigma}_{xz}\\\
{\sigma}_{xy}\end{array}\right]_{n+k}=\begin{bmatrix}S_{11}(k)&S_{12}(k)\\\
S_{21}(k)&S_{22}(k)\end{bmatrix}\left[\begin{array}[]{c}{v}_{x}\\\ {v}_{y}\\\
{v}_{z}\\\ {\sigma}_{xx}\\\ {\sigma}_{yy}\\\ {\sigma}_{zz}\\\ {\sigma}_{yz}\\\
{\sigma}_{xz}\\\ {\sigma}_{xy}\end{array}\right]_{n}~{},$ (C.58)

and last equation is formally equivalent to equation 72.

##  Appendix D A 5 $\times 5$ first-order system for the VTI acoustic wave
equation

In this section we derive the 5 $\times 5$ first-order system for the VTI
acoustic wave equation. The 5 $\times 5$ first-order system is based on 5:

$\displaystyle\left[\begin{array}[]{c}\dot{v}_{x}\\\ \dot{v}_{y}\\\
\dot{v}_{z}\\\ \dot{\sigma}_{H}\\\
\dot{\sigma}_{M}\end{array}\right]=\begin{bmatrix}0&A\\\
B&0\end{bmatrix}\left[\begin{array}[]{c}{v}_{x}\\\ {v}_{y}\\\ {v}_{z}\\\
{\sigma}_{H}\\\ {\sigma}_{M}\end{array}\right]~{},$ (D.11)

where

$\displaystyle A=\begin{bmatrix}(1/\rho)\;\partial_{x}&0\\\
(1/\rho)\;\partial_{y}&0\\\
\frac{\partial_{z}}{\rho(1+2\eta)\sqrt{1+2\delta}}&\frac{\sqrt{2\eta}\partial_{z}}{\rho(1+2\eta)\sqrt{1+2\delta}}\end{bmatrix}~{},\;\;\;$
(D.12) $\displaystyle B=\begin{bmatrix}\rho
V^{2}_{Pz}(1+2\varepsilon)\partial_{x}&\rho
V^{2}_{Pz}(1+2\varepsilon)\partial_{y}&\rho
V^{2}_{Pz}\sqrt{1+2\delta}\;\partial_{z}\\\ 0&0&\rho
V^{2}_{Pz}\sqrt{2\eta}\sqrt{1+2\delta}\;\partial_{z}\end{bmatrix}~{},$

$\mathbf{v}=[v_{x},v_{y},v_{z}]^{T}$ is the particle velocity vector,
$\bm{\sigma}=[\sigma_{H}$, $\sigma_{M}]^{T}$ are the horizontal and mixed
horizontal and vertical components of the stress tensor $\bm{\sigma}$, $\rho$
is density, $\varepsilon$, $\delta$ and $\eta$ are the anisotropy parameters.
For details of definitions and derivations of equation D see Bube et al. [5].
Equation D.11 is typically discretized the following way:

$\displaystyle\mathbf{v}_{n+1/2}$ $\displaystyle=$
$\displaystyle\mathbf{v}_{n-1/2}+\Delta t\;A\;\bm{\sigma}_{n}~{},$
$\displaystyle\bm{\sigma}_{n+1}$ $\displaystyle=$
$\displaystyle\bm{\sigma}_{n}+\Delta t\;B\;\mathbf{v}_{n+1/2}~{},$ (D.13)

Note that $\mathbf{v}$ and $\bm{\sigma}$ are on staggered grid relative to
each other and the material properties are sampled on their respective grid.
Expressing the $k$ time step shifted components $\mathbf{v}_{n-1/2+k}$ and
$\bm{\sigma}_{n+k}$ from $\mathbf{v}_{n-1/2}$ and $\bm{\sigma}_{n}$ yields

$\displaystyle\left[\begin{array}[]{c}{v}_{x}\\\ {v}_{y}\\\ {v}_{z}\\\
{\sigma}_{H}\\\ {\sigma}_{M}\end{array}\right]_{n+k}=\begin{bmatrix}I&\Delta
t\;A\\\ \Delta t\;B&I+(\Delta
t)^{2}BA\end{bmatrix}^{k}\left[\begin{array}[]{c}{v}_{x}\\\ {v}_{y}\\\
{v}_{z}\\\ {\sigma}_{H}\\\ {\sigma}_{M}\end{array}\right]_{n}~{},$ (D.24)

where for simplicity we ignored the staggered grid notation when denoting the
state vector time steps. Here the system matrix is also the superstep
propagator matrix in equation 80

$\displaystyle\left[\begin{array}[]{c}{v}_{x}\\\ {v}_{y}\\\ {v}_{z}\\\
{\sigma}_{H}\\\
{\sigma}_{M}\end{array}\right]_{n+k}=\begin{bmatrix}S_{11}(k)&S_{12}(k)\\\
S_{21}(k)&S_{22}(k)\end{bmatrix}\left[\begin{array}[]{c}{v}_{x}\\\ {v}_{y}\\\
{v}_{z}\\\ {\sigma}_{H}\\\ {\sigma}_{M}\end{array}\right]_{n}~{},$ (D.35)

and last equation is formally equivalent to equation 72.

##  Appendix E Second-order system for the acoustic wave equation

The second-order system for the acoustic wave equation is expressed as based
on 45:

$\ddot{p}=V^{2}\Delta p$ (E.1)

where $p$ is pressure and $V$ is P-wave velocity. Equation E.1 is typically
discretized the following way:

$p_{n+1}=2\,p_{n}-p_{n-1}+\Delta t\;V^{2}\Delta p_{n}=Gp_{n}-p_{n-1}$ (E.2)

where $G=G_{1}=2+\Delta t\;V^{2}\Delta$ is a single step operator. In general,
let $G_{k}$ to be the propagator matrix for advancing $k$ time levels, taking
into account that we need to store two adjacent time levels to be able to
advance in time for second-order systems. $G_{k}$ is defined by

$G_{k}(p_{0})=p_{k}~{},$ (E.3)

where $p_{k}$ is the solution of equation E.2 at $n=k$ with initial conditions
$p_{-1}=0$ and $p_{0}$.

The $k^{th}$ time step based on 65:

$\displaystyle\left[\begin{array}[]{c}p_{n+k-1}\\\
p_{n+k}\end{array}\right]=\begin{bmatrix}0&I\\\
-I&G\end{bmatrix}^{k}\left[\begin{array}[]{c}p_{n-1}\\\
p_{n}\end{array}\right]~{},$ (E.8)

or

$\displaystyle\left[\begin{array}[]{c}p_{n+k-1}\\\
p_{n+k}\end{array}\right]=\begin{bmatrix}S_{11}(k)&S_{12}(k)\\\
S_{21}(k)&S_{22}(k)\end{bmatrix}\left[\begin{array}[]{c}p_{n-1}\\\
p_{n}\end{array}\right]~{},$ (E.13)

where the last equation is formally equivalent to equation 79.

Next, we show derive properties of $\bm{S}(k)$ and $G_{k}$ that will indicate
reverse time propagation using these rotated block matrices. Block matrix
multiplication from equation (E.8) gives

$\displaystyle S(2)=S(1)^{2}=\begin{bmatrix}-I&G\\\
-G&G^{2}-I\end{bmatrix}~{},$ (E.14) $\displaystyle
S(3)=S(1)^{3}=\begin{bmatrix}-G&G^{2}-I\\\
-(G^{2}-I)&G^{3}-2G\end{bmatrix}~{}.$ (E.15)

For $k\geq 1$, since $G$ is independent of $k$, we see from the definition of
$G_{k}$ that the lower right block of $S(k)$ is $G_{k}$: $p_{n+k}$ would be
$G_{k}(p_{n})$ if $p_{n-1}$ were zero. Similarly, the upper right block of
$S(k)$ is $G_{k-1}$: $p_{n+k-1}$ would be $G_{k-1}(p_{n})$ if $p_{n-1}$ were
zero. For $k\geq 2$, consider the solution of equation (E.2) with a given
$p_{n-1}$ and zero $p_{n}$. By equation (E.2), this solution has zero $p_{n}$
and $p_{n+1}=-p_{n-1}$; advancing this solution $k-1$ time steps, $p_{n+k}$
would be $G_{k-1}(p_{n+1})=-G_{k-1}(p_{n-1})$ and $p_{n+k-1}$ would be
$G_{k-2}(p_{n+1})=-G_{k-2}(p_{n-1})$. Thus the lower left block of $S(k)$ is
$-G_{k-1}$, the upper left block of $S(k)$ is $-G_{k-2}$, and for an arbitrary
$k$ time step we get

$\displaystyle S(k)=\begin{bmatrix}-G_{k-2}&G_{k-1}\\\
-G_{k-1}&G_{k}\end{bmatrix}~{}.$ (E.16)

Now we derive the recurrence relation for $G_{k}$. Consider time step $S(k+1)$

$\displaystyle S(k+1)=\begin{bmatrix}-G_{k-1}&G_{k}\\\
-G_{k}&G_{k+1}\end{bmatrix}=S(1)S(k)=\begin{bmatrix}0&I\\\
-I&G_{1}\end{bmatrix}\begin{bmatrix}-G_{k-2}&G_{k-1}\\\
-G_{k-1}&G_{k}\end{bmatrix}=\begin{bmatrix}-G_{k-1}&G_{k}\\\
G_{k-2}-G_{1}G_{k-1}&G_{1}G_{k}-G_{k-1}\end{bmatrix}$ (E.17)

we conclude that for $k\geq 1$

$G_{k+1}=G_{1}G_{k}-G_{k-1}~{}.$ (E.18)

Another form of equation E.18 is

$G^{2}_{k+1}-G_{k-2}G_{k}=I~{},$ (E.19)

that can be shown directly by induction using the recurrence relation E.18
twice.

Starting with $G_{0}=I$ and $G_{1}=2+V^{2}\Delta$, we get for $G_{k}$

$\displaystyle\begin{aligned} G_{2}=G_{1}^{2}-&I~{},\\\
G_{3}=G_{1}^{3}-&2G_{1}~{},\\\ G_{4}=G_{1}^{4}-&3G_{1}^{2}+I~{},\end{aligned}$
(E.20)

etc. We see that each $G_{k}$ is a polynomial $q_{k}(G_{1})$ in $G_{1}$, where
$q_{k}$ is a monic polynomial of degree $k$: $q_{0}(x)=1$, $q_{1}(x)=x$,
$q_{2}(x)=x^{2}-1$, $q_{3}(x)=x^{3}-2x$, and in general,
$q_{k+1}(x)=xq_{k}(x)-q_{k-1}(x)$ for $k\geq 1$. In particular, all of the
$G_{k}$’s commute with each other: for $k\geq 0$ and $n\geq 0$,

$G_{k}G_{n}=G_{n}G_{k}~{}.$ (E.21)

This monic polynomial can be illustrated in a triangle similar to Pascal’s
Triangle (Figure 8).

Figure 8: Wave equation triangle

Now we show that the inverse of the superstep propagator matrix is the
reverse-time superstep propagator. Consider equation

$\displaystyle S(1)^{-1}=\begin{bmatrix}G_{1}&-I\\\ I&0\end{bmatrix}~{}.$
(E.22)

Notice that $S(1)^{-1}$ is the same as $S(1)$ with the block columns switched
and the block rows switched or 180 degrees rotated. This observation can be
written in matrix form. Let $R$ be the 2 x 2 block matrix

$\displaystyle R=\begin{bmatrix}0&I\\\ I&0\end{bmatrix}~{};$ (E.23)

$R$ is a symmetric permutation matrix, and $R^{-1}=R^{T}=R$. Either from the
observation above about switching block columns and block rows, or by direct
computation, we have

$S(1)^{-1}=RS(1)R~{}.$ (E.24)

From this, we can compute $S(k)^{-1}$ for $k\geq 2$:

$\displaystyle
S(k)^{-1}=(S(1)^{k})^{-1}=(S(1)^{-1})^{k}=(RS(1)R)^{k}=(RS(1)R)(RS(1)R)\cdots(RS(1)R)=(RS(1)^{k}R)$
$\displaystyle=(RS(k)R)=\begin{bmatrix}0&I\\\
I&0\end{bmatrix}\begin{bmatrix}-G_{k-2}&G_{k-1}\\\
-G_{k-1}&G_{k}\end{bmatrix}\begin{bmatrix}0&I\\\
I&0\end{bmatrix}=\begin{bmatrix}G_{k}&-G_{k-1}\\\
G_{k-1}&-G_{k-2}\end{bmatrix}$ (E.25)

indicating reverse timestep propagation.