Massive Gauge Bosons in Yang-Mills Theory without Higgs Mechanism

Xin-Bing<EMAIL_ADDRESS>
Shanghai United Center for Astrophysics (SUCA),

Shanghai Normal University, No.100 Guilin Road, Shanghai 200234, China

Abstract

Two kinds of Yang-Mills fields are found upon the concepts of mass eigenstate
and nonmass eigenstate. The Yang-Mills fields of the first kind were proposed
by Yang and Mills, which couple to the mass eigenstates with the same rest
mass, whose gauge bosons are massless. I find that there are second kind of
Yang-Mills fields, which are constructed on a five-dimensional manifold. Only
the nonmass eigenstates couple to the Yang-Mills fields of the second kind,
which are the nonmass eigenstates as well and composed of mass eigenstates of
gauge bosons. The mass eigenstates of the Yang-Mills fields of the second kind
live in the four-dimensional spacetime, the corresponding gauge bosons of
which may be massive. The $SU(2)\times U(1)$ gauge fields of the second kind
are studied carefully, whose gauge bosons, which are the mass eigenstates, are
the $W^{\pm}$, $Z^{0}$ and photon fields. The rest masses of $W^{\pm}$ and
$Z^{0}$ obtained are the same as that given by the Glashow-Salam-Weinberg
model of electroweak interactions. It is discussed that this model should be
renormalizable.

PACS numbers: 11.15.-q, 11.10.Kk, 12.60.-i

55 years ago, Yang and Mills constructed the gauge field theory of non-Abelian
group, which has become the most fundamental content in quantum field theory.
Upon the principle that physical laws should be covariant under the local
isospin rotation they proposed the $SU(2)$ Yang-Mills theory [1]. But they
could not obtain the massive gauge bosons then. About 10 years later, an
ingenious trick called the Higgs mechanism was independently invented by Higgs
and Englert and Brout [2], who introduced a scalar field and the spontaneous
symmetry broken mechanism of vacuum by fixing a vacuum expectation value of
the scalar field and make the intermediate vector bosons obtain masses.

Based on the Yang-Mills fields and the Higgs mechanism, Glashow, Salam and
Weinberg etc. proposed a renormalizable theory unifying the weak and
electromagnetic interactions, namely $SU_{L}(2)\times U_{Y}(1)$ gauge theory
[3]. Although this electroweak theory had predicted the masses of intermediate
vector bosons, which were confirmed by experiments, there are still several
unconfirmed predictions or conflicting phenomena in it. e.g. Firstly,
experimenters have not found any hints of the Higgs boson till now; Secondly,
a lot of recent experiments imply that the neutrinos should be massive and be
mixed [4]. Here I discuss a model to give the massive gauge bosons in Yang-
Mills theory without Higgs mechanism.

In this letter, the signature of spacetime metric
$\eta_{\mu\nu}(\mu,\nu=0,1,2,3)$ is $(+,-,-,-)$, and the spacetime coordinates
are described by the contravariant four-vector $x^{\mu}$ ($\hbar=c=1$ is
adopted). In Ref.[5], the rest mass operator†††I use
$\partial_{z}\equiv\frac{\partial}{\partial
z}~{},~{}\partial_{\mu}\equiv\frac{\partial}{\partial x^{\mu}}$ and
$\partial_{\alpha}\equiv\frac{\partial}{\partial x^{\alpha}}$.

$\hat{m}=-i\partial_{z}$ (1)

is defined by introducing an extra parameter $z$ besides of the spacetime
coordinates $x^{\mu}$. From the mathematical point of view, $z$ and $x^{\mu}$
establish a five-dimensional manifold. The definition of the rest mass
operator leads to a theorem that a field ${\cal F}(x,z)$ is massless if and
only if ${\cal F}(x,z)$ is $z$-independent [5]. Hence the massless
gravitational field, the electromagnetic field and $SU(3)$ gauge fields in
$QCD$ are all $z$-independent, who live in the $z=0$ brane of five-dimensional
manifold.

The Lagrangian of a nonmass eigenstate ${\Phi}(x,z)$ of free
spin-$\frac{1}{2}$ fields is of the form‡‡‡${\cal L}_{1n}$, ${\cal L}_{1m}$
denote the Lagrangian of one nonmass eigenstate or one mass eigenstate
respectively. ${\cal L}_{2n}$, ${\cal L}_{2m}$ have the similar meanings.

$\displaystyle{\cal
L}_{1n}={\bar{\Phi}}(x,z)\left(i\gamma^{\mu}{\partial}_{\mu}+i{\partial}_{z}\right)\Phi(x,z)~{},$
(2)

here $\bar{\Phi}\equiv{\Phi}^{{\dagger}}\gamma^{0}$ is called the spinor
adjoint to $\Phi$. I indicated that the mass eigenstate of a
spin-$\frac{1}{2}$ field satisfies $\Phi(x,z)=e^{imz}\phi(x)$ in Ref.[5],
where $m$ is the rest mass. Therefore one can obtain the Lagrangian of the
mass eigenstate of a free spin-$\frac{1}{2}$ field from (2), that is

$\displaystyle{\cal
L}_{1m}={\bar{\phi}}(x)\left(i\gamma^{\mu}{\partial}_{\mu}-m\right)\phi(x)~{},$
(3)

where $\bar{\phi}\equiv{\phi}^{{\dagger}}\gamma^{0}$ is the spinor adjoint to
$\phi$.

Let’s consider a quantum field system in which two different nonmass
eigenstates $\Psi_{1}(x,z)$ and $\Psi_{2}(x,z)$ of free spin-$\frac{1}{2}$
fields form an isospin doublet as follows

$\displaystyle\Psi(x,z)=\left(\begin{array}[]{c}\Psi_{1}(x,z)\\\
\Psi_{2}(x,z)\end{array}\right)~{}.$ (6)

So the Lagrangian of two nonmass eigenstates of free spin-$\frac{1}{2}$ fields
is

$\displaystyle{\cal
L}_{2n}={\bar{\Psi}}(x,z)\left(i\gamma^{\mu}{\partial}_{\mu}+i{\partial}_{z}\right)\Psi(x,z)~{}.$
(7)

Here Let’s first consider a special case: if $\Psi_{1}(x,z)$ and
$\Psi_{2}(x,z)$ are mass eigenstates with the same rest mass, then
$\Psi_{1}(x,z)=e^{imz}\psi_{1}(x)$, and $\Psi_{2}(x,z)=e^{imz}\psi_{2}(x)$,
where $m$ is the rest mass. I can therefore obtain $\Psi(x,z)=e^{imz}\psi(x)$
by defining

$\displaystyle\psi(x)=\left(\begin{array}[]{c}\psi_{1}(x)\\\
\psi_{2}(x)\end{array}\right)~{}.$ (10)

Hence the Lagrangian of a quantum field system where two mass eigenstates
$e^{imz}\psi_{1}(x)$ and $e^{imz}\psi_{2}(x)$ form an isospin doublet is
acquired from (7), namely

$\displaystyle{\cal
L}_{2m}={\bar{\psi}}(x)\left(i\gamma^{\mu}{\partial}_{\mu}-m\right)\psi(x)~{}.$
(11)

The largest inner gauge symmetry group in this system is obviously
$SU(2)\times U(1)$. The total Lagrangian of this system reads

$\displaystyle{\cal
L}_{2mt}=i{\bar{\psi}}\gamma^{\mu}(\partial_{\mu}-ig^{\prime}{\bf T}\cdot{\bf
B}_{\mu})\psi-e{\bar{\psi}}\gamma^{\mu}A_{\mu}\psi$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}-m{\bar{\psi}}\psi-\frac{1}{4}~{}{\tilde{\bf
F}}_{\mu\nu}\cdot{\tilde{\bf
F}}^{\mu\nu}-~{}\frac{1}{4}~{}{\tilde{E}}_{\mu\nu}{\tilde{E}}^{\mu\nu}~{},$
(12)

where $e$, $g^{\prime}$ are the coupling constants of $U(1)$ and $SU(2)$ gauge
fields respectively, and the dot “$\cdot$” denotes a scalar product in the
isospace. In this case, ${\bf T}\cdot{\bf B}_{\mu}$ means

$\displaystyle{\bf T}\cdot{\bf
B}_{\mu}=T^{1}B^{1}_{\mu}+T^{2}B^{2}_{\mu}+T^{3}B^{3}_{\mu}~{},$ (13)

where $T^{a},~{}a=1,2,3$ are the generators of $SU(2)$ group, which are
written as

$T^{a}=\frac{1}{2}\tau^{a}~{},$ (14)

where $\tau^{a}$ are the traceless matrices

$\displaystyle\tau^{1}=\left(\begin{array}[]{cc}0&1\\\
1&0\end{array}\right)~{},~{}~{}\tau^{2}=\left(\begin{array}[]{cc}0&-i\\\
i&0\end{array}\right)~{},~{}~{}\tau^{3}=\left(\begin{array}[]{cc}1&0\\\
0&-1\end{array}\right)~{},$ (21)

known as the Pauli matrices. They obey the commutation relations

$[\tau^{a},\tau^{b}]=2i\sum^{3}_{c=1}\varepsilon_{abc}\tau^{c}~{}.$ (22)

Here $\varepsilon_{abc}$ is the totally antisymmetry tensor in 3-dimensions.

In Yang-Mills theory [1], ${\bf T}\cdot{\bf B}_{\mu}$ is called the $SU(2)$
gauge field, and its field strength tensor is of the form

$\displaystyle{\tilde{\bf F}}_{\mu\nu}\cdot{\bf T}=\partial_{\mu}({\bf
B}_{\nu}\cdot{\bf T})-\partial_{\nu}({\bf B}_{\mu}\cdot{\bf
T})-{i}g^{\prime}[{\bf B}_{\mu}\cdot{\bf T},{\bf B}_{\nu}\cdot{\bf T}]~{}.$
(23)

Hence the field strength ${\tilde{\bf F}}_{\mu\nu}$ satisfies

${\tilde{\bf F}}_{\mu\nu}=\partial_{\mu}{\bf B}_{\nu}-\partial_{\nu}{\bf
B}_{\mu}+g^{\prime}{\bf B}_{\mu}\times{\bf B}_{\nu}~{}.$ (24)

I use $A_{\mu}$ to denote the $U(1)$ gauge field. The field strength of the
$U(1)$ gauge field is defined by

${\tilde{E}}_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}~{}.$ (25)

We are very familiar with above $SU(2)$ gauge fields and $U(1)$ gauge field
which have been the fundamental content of quantum field theories. The gauge
bosons are massless. From a viewpoint of the rest mass operator, above $SU(2)$
gauge fields couple to the mass eigenstates with the same rest mass. I call
this Yang-Mills fields the first kind. The $SU(3)$ gauge fields in $QCD$
obviously belong to this kind.

From now on I will study another kind of Yang-Mills fields carefully. To make
the gauge invariance explicit, let’s formally introduce the extra dimension
$x^{4}$ as follows

$x^{4}=-x_{4}=z~{},~{}~{}~{}~{}\gamma^{4}=-\gamma_{4}={\bf 1}~{}.$ (26)

Hence the Lagrangian ${\cal L}_{2n}$ is rewritten as

$\displaystyle{\cal
L}_{2n}=i{\bar{\Psi}}(x,z)\gamma^{\alpha}{\partial}_{\alpha}\Psi(x,z)~{},~{}~{}~{}~{}\alpha=0,1,2,3,4~{}.$
(27)

Since two different nonmass eigenstates $\Psi_{1}(x,z)$ and $\Psi_{2}(x,z)$
form an isospin doublet, I can consider a local isospin rotation logically
similar to what Yang and Mills did in their original paper. That is

$\displaystyle{\Psi}^{\prime}(x,z)=S(x,z)\Psi(x,z)~{},$ (28)

where $S(x,z)$ is a $2\times 2$ matrix. To make sure that the probability
density ${\bar{\Psi}}(x,z)\Psi(x,z)$ is invariant under above rotation (28),
the matrix $S(x,z)$ must be unitary with unit determinant

$\displaystyle S^{{\dagger}}(x,z)S(x,z)=1~{}.$ (29)

All the matrices satisfy this condition generate the group $SU(2)$, which is a
non-Abelian Lie group. The transformation (28) directly means that

$\displaystyle{\bar{\Psi}}^{\prime}(x,z)={\bar{\Psi}}(x,z)S^{{\dagger}}(x,z)~{}.$
(30)

The matrix $S(x,z)$ can be written in the form

$\displaystyle
S(x,z)=\exp\left(i\sum_{a=1}^{3}\frac{{\tau}^{a}}{2}{\Theta}^{a}(x,z)\right)~{}.$
(31)

To discuss the gauge invariance, here I introduce the gauge-invariant
derivative

$\displaystyle\hat{D}_{\alpha}=\partial_{\alpha}-ig_{1}{\bf T}\cdot{\bf
W}_{\alpha}(x,z)~{},$ (32)

where $g_{1}$ is the coupling constant of $SU(2)$ gauge fields, and

$\displaystyle{\bf T}\cdot{\bf W}_{\alpha}(x,z)=\sum_{a=1}^{3}T^{a}{\bf
W}_{\alpha}^{a}(x,z)~{}.$ (33)

Invariance requires that

$\displaystyle(\partial_{\alpha}-ig_{1}{\bf T}\cdot{\bf
W}_{\alpha}^{\prime}){\Psi}^{\prime}=S(\partial_{\alpha}-ig_{1}{\bf
T}\cdot{\bf W}_{\alpha})\Psi~{}.$ (34)

Combining (28) and (34), I obtain the gauge transformation on ${\bf
W}_{\alpha}$:

$\displaystyle{\bf T}\cdot{\bf W}_{\alpha}^{\prime}=S{\bf T}\cdot{\bf
W}_{\alpha}S^{-1}+\frac{i}{g_{1}}S\left(\partial_{\alpha}S^{-1}\right)~{}.$
(35)

In analogy to the procedure of obtaining gauge invariant field strengths in
electromagnetic case, I define now

$\displaystyle{{\bf F}}_{\alpha\beta}\cdot{\bf T}=\sum_{a=1}^{3}{\bf
F}^{a}_{\alpha\beta}{T}^{a}=\hat{D}_{\alpha}({\bf T}\cdot{\bf
W}_{\beta})-\hat{D}_{\beta}({\bf T}\cdot{\bf W}_{\alpha})$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}=\partial_{\alpha}({\bf
W}_{\beta}\cdot{\bf T})-\partial_{\beta}({\bf W}_{\alpha}\cdot{\bf
T})-{i}g_{1}\left[{\bf W}_{\alpha}\cdot{\bf T},{\bf W}_{\beta}\cdot{\bf
T}\right]$ $\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}=(\partial_{\alpha}{\bf
W}_{\beta})\cdot{\bf T}-(\partial_{\beta}{\bf W}_{\alpha})\cdot{\bf
T}+g_{1}\sum_{abc}{\bf W}^{a}_{\alpha}{\bf
W}^{b}_{\beta}\varepsilon_{abc}T^{c}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}=(\partial_{\alpha}{\bf
W}_{\beta}-\partial_{\beta}{\bf W}_{\alpha}+g_{1}{\bf W}_{\alpha}\times{\bf
W}_{\beta})\cdot{\bf T}~{}.$ (36)

Therefore the isovector of field strengths is

$\displaystyle{{\bf F}}_{\alpha\beta}=\partial_{\alpha}{\bf
W}_{\beta}-\partial_{\beta}{\bf W}_{\alpha}+g_{1}{\bf W}_{\alpha}\times{\bf
W}_{\beta}~{}.$ (37)

One easily shows from the equation (35) that

$\displaystyle{\bf F}^{\prime}_{\alpha\beta}\cdot{\bf T}=S{{\bf
F}}_{\alpha\beta}\cdot{\bf T}S^{-1}~{}.$ (38)

I obtain a gauge invariant Lagrangian by performing the trace over the isospin
indices:

$\displaystyle{\cal L}_{SU(2)}=-\frac{1}{2}{\rm Tr}\\{({{\bf
F}}_{\alpha\beta}\cdot{\bf T})({{\bf F}}^{\alpha\beta}\cdot{\bf T})\\}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}=-\frac{1}{4}{{\bf
F}}_{\alpha\beta}\cdot{{\bf F}}^{\alpha\beta}=-\frac{1}{4}\sum_{a=1}^{3}{{\bf
F}}^{a}_{\alpha\beta}{{\bf F}}^{a\alpha\beta}~{}.$ (39)

Considering the couplings between fermions and gauge bosons and the self-
couplings of gauge bosons, one can get the complete Lagrangian as follows

$\displaystyle{\cal L}_{2nt}={\cal L}_{2n}+{\cal L}_{int}+{\cal L}_{SU(2)}$
$\displaystyle~{}~{}~{}~{}~{}~{}=i{\bar{\Psi}}\gamma^{\alpha}(\partial_{\alpha}-ig_{1}{\bf
T}\cdot{\bf W}_{\alpha})\Psi-\frac{1}{4}{{\bf F}}_{\alpha\beta}\cdot{{\bf
F}}^{\alpha\beta}~{}.$ (40)

In order to build a foundation for setting up an electroweak model without
Higgs mechanism, I discuss the $SU(2)\times U(1)$ gauge fields in this letter.
The $U(1)$ gauge field that couples to a nonmass eigenstate has been studied
in my preceding paper [5]. The Lagrangian (27) shows me that the maximal gauge
groups for this quantum field system are $SU(2)\times U(1)$. I have introduced
the $SU(2)$ gauge fields in this system, now I put in the $U(1)$ gauge field.
Let us multiply the nonmass eigenstates $\Psi(x,z)$ by a local phase
$e^{i\Theta(x,z)}$, namely

$\Psi^{\prime\prime}=~{}e^{i\Theta(x,z)}\Psi~{},~{}{\bar{\Psi}}^{\prime\prime}=~{}e^{-i\Theta(x,z)}{\bar{\Psi}}~{}.$
(41)

According to the discussion in Ref.[5], I introduce the $U(1)$ gauge field of
the second kind, that is ${\bf X}_{\alpha}(x,z)$. Under the transformation of
(41), ${\bf X}_{\alpha}(x,z)$ transforms as

${\bf X}^{\prime\prime}_{\alpha}(x,z)={\bf
X}_{\alpha}(x,z)+\frac{1}{g_{2}}{\partial_{\alpha}}\Theta(x,z)~{},$ (42)

here $g_{2}$ is the coupling constant of $U(1)$. The strength tensor of $U(1)$
gauge field is of the form

${\bf E}_{\alpha\beta}(x,z)={\partial}_{{\alpha}}{\bf
X}_{\beta}(x,z)-{\partial}_{\beta}{\bf X}_{\alpha}(x,z)~{},$ (43)

which is invariant under the transformations of (41) and (42). Therefore the
total Lagrangian including the $U(1)$ gauge field of the second kind is
written as

$\displaystyle{\cal L}_{total}={\cal L}_{2n}+{\cal L}_{int}+{\cal
L}_{U(1)}+{\cal L}_{SU(2)}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}=i{\bar{\Psi}}\gamma^{\alpha}(\partial_{\alpha}-ig_{2}{\bf
X}_{\alpha}-ig_{1}{\bf T}\cdot{\bf W}_{\alpha})\Psi$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-\frac{1}{4}{{\bf
E}}_{\alpha\beta}{\bf E}^{\alpha\beta}-\frac{1}{4}{{\bf
F}}_{\alpha\beta}\cdot{{\bf F}}^{\alpha\beta}~{}.$ (44)

Obviously the total Lagrangian is also invariant under the transformations of
(41) and (42). The gauge covariance requires that ${\bf X}_{\alpha}$ and ${\bf
T}\cdot{\bf W}_{\alpha}$ are all nonmass eigenstates.

Till now I merely constructed the $SU(2)\times U(1)$ gauge fields on a five-
dimensional manifold, which is quite the same as the gauge fields proposed by
Yang and Mills. Yes, from the five-dimensional point of view, the gauge bosons
are massless in above discussed $SU(2)\times U(1)$ gauge fields since there
are no mass term in the total Lagrangian (44). But, things will be quite
different when I discuss them from the viewpoint of $z=0$ brane.

I have pointed out that the gravitational field, the electromagnetic field and
$SU(3)$ gauge fields in $QCD$ are living in the $z=0$ brane of five-
dimensional manifold. Also I have proved that the $z$-independent
electromagnetic field, gravitational field and $SU(3)$ gauge fields only
couple to the mass eigenstates. Therefore I can find that the mass eigenstates
coupled by the gravitation, the electromagnetic field and the gluon fields are
also living in the $z=0$ 4-dimensional brane.

It is indicated that the nonmass eigenstate is composed of mass eigenstates
[5]. To discuss the physical properties of the mass eigenstates who compose
the gauge fields ${\bf X}_{\alpha}$ and ${\bf W}_{\alpha}$, I write out the
spacetime component and $z$-related component of gauge fields separately.
Therefore

$\displaystyle{\bf X}_{\alpha}(x,z)\equiv({\bf X}_{\mu}(x,z),{\bf
X}_{z}(x,z))~{},$ (45) $\displaystyle{\bf W}_{\alpha}(x,z)\equiv({\bf
W}_{\mu}(x,z),{\bf W}_{z}(x,z))~{}.$ (46)

To list the components of ${\bf W}_{\alpha}(x,z)$ manifestly, I rewrite (46)
as

${\bf W}^{a}_{\alpha}(x,z)\equiv({\bf W}^{a}_{\mu}(x,z),{\bf
W}^{a}_{z}(x,z))~{},~{}a=1,2,3~{}.$ (47)

After that, the strength tensor ${\bf E}_{\alpha\beta}(x,z)$ is
correspondingly divided into three parts

$\displaystyle{\bf E}_{\mu\nu}(x,z)=\partial_{\mu}{\bf
X}_{\nu}(x,z)-\partial_{\nu}{\bf X}_{\mu}(x,z)~{},$ (48) $\displaystyle{\bf
E}_{\mu z}(x,z)=-{\bf E}_{z\mu}(x,z)=\partial_{\mu}{\bf
X}_{z}(x,z)-\partial_{z}{\bf X}_{\mu}(x,z)~{},$ (49) $\displaystyle{\bf
E}_{zz}(x,z)=\partial_{z}{\bf X}_{z}(x,z)-\partial_{z}{\bf X}_{z}(x,z)\equiv
0~{}.$ (50)

Surely one can also get the decomposition of the strength tensor ${\bf
F}_{\alpha\beta}(x,z)$ as follows

$\displaystyle{{\bf F}}_{\mu\nu}(x,z)=\partial_{\mu}{\bf
W}_{\nu}-\partial_{\nu}{\bf W}_{\mu}+g_{1}{\bf W}_{\mu}\times{\bf
W}_{\nu}~{},$ (51) $\displaystyle{\bf F}_{\mu z}(x,z)=-{\bf
F}_{z\mu}(x,z)=\partial_{\mu}{\bf W}_{z}-\partial_{z}{\bf W}_{\mu}+g_{1}{\bf
W}_{\mu}\times{\bf W}_{z}~{},$ (52) $\displaystyle{\bf F}_{zz}(x,z)\equiv
0~{}.$ (53)

Now let’s consider the movement of gauge bosons. Firstly, the interaction term
${\cal L}_{int}$ in total Lagrangian (44) shows that the movement of gauge
bosons is decided by the momentum of $\Psi$ and $\bar{\Psi}$. The $\Psi$ and
$\bar{\Psi}$ are nonmass eigenstates, who are composed of mass eigenstates
that are living in the $z=0$ brane and moving along the $z=0$ brane, hence
gauge bosons must move along the $z=0$ brane. Secondly, once the gauge bosons
are produced, they are constrained by gravitation, which is living in the
$z=0$ brane. Consequently

${\bf W}_{z}(x,z)=0~{},~{}{\bf X}_{z}(x,z)=0~{}.$ (54)

Hence the equations (49) and (52) reduce to

$\displaystyle{\bf E}_{\mu z}(x,z)=-{\bf E}_{z\mu}(x,z)=-\partial_{z}{\bf
X}_{\mu}(x,z)~{},$ (55) $\displaystyle{\bf F}_{\mu z}(x,z)=-{\bf
F}_{z\mu}(x,z)=-\partial_{z}{\bf W}_{\mu}(x,z)~{}.$ (56)

The spacetime components ${\bf X}_{\mu}(x,z)$ and ${\bf W}_{\mu}(x,z)$ are
nonmass eigenstates, which are linear combinations of mass eigenstates. I
define the mass eigenstates of bosons $W^{\pm}$ by

$\displaystyle{\bf W}^{+}_{\mu}(x,z)=\frac{1}{\sqrt{2}}\left({\bf
W}^{1}_{\mu}(x,z)-i{\bf W}^{2}_{\mu}(x,z)\right)~{},$ (57) $\displaystyle{\bf
W}^{-}_{\mu}(x,z)=\frac{1}{\sqrt{2}}\left({\bf W}^{1}_{\mu}(x,z)+i{\bf
W}^{2}_{\mu}(x,z)\right)~{}.$ (58)

When the mass eigenstates of $W^{\pm}$ are expressed by

$\displaystyle{\bf W}^{+}_{\mu}(x,z)=e^{im_{W}z}{W}^{+}_{\mu}(x)~{},~{}{\bf
W}^{-}_{\mu}(x,z)=e^{im_{W}z}{W}^{-}_{\mu}(x)~{},$ (59)

$m_{W}$ being the rest mass of $W^{\pm}$, the nonmass eigenstates ${\bf
W}^{1}_{\mu}(x,z)$ and ${\bf W}^{2}_{\mu}(x,z)$ are manifestly given by

$\displaystyle{\bf W}^{1}_{\mu}(x,z)=\frac{1}{\sqrt{2}}\left({\bf
W}^{+}_{\mu}(x,z)+{\bf W}^{-}_{\mu}(x,z)\right)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}=\frac{1}{\sqrt{2}}e^{im_{W}z}\left({W}^{+}_{\mu}(x)+{W}^{-}_{\mu}(x)\right)~{},$
(60)

and

$\displaystyle{\bf W}^{2}_{\mu}(x,z)=\frac{i}{\sqrt{2}}\left({\bf
W}^{+}_{\mu}(x,z)-{\bf W}^{-}_{\mu}(x,z)\right)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}=\frac{i}{\sqrt{2}}e^{im_{W}z}\left({W}^{+}_{\mu}(x)-{W}^{-}_{\mu}(x)\right)~{}.$
(61)

The boson fields ${\bf W}^{+}_{\mu}(x,z)$, ${\bf W}^{-}_{\mu}(x,z)$, ${\bf
Z}_{\mu}(x,z)$ and photon field $A_{\mu}(x)$, which are mass eigenstates,
constitute a complete Hilbert space. From Ref.[5], I know that

$\displaystyle{\bf Z}_{\mu}(x,z)=e^{im_{Z}z}Z_{\mu}(x)~{},$ (62)

here $m_{Z}$ is the rest mass of boson $Z^{0}$. The nonmass eigenstates ${\bf
W}^{3}_{\mu}(x,z)$ and ${\bf X}_{\mu}(x,z)$ are the linear combinations of
${\bf Z}_{\mu}(x,z)$ and $A_{\mu}(x)$, namely

$\displaystyle{\bf W}^{3}_{\mu}(x,z)=\sin\theta_{W}A_{\mu}+\cos\theta_{W}{\bf
Z}_{\mu}(x,z)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}=\sin\theta_{W}A_{\mu}+\cos\theta_{W}e^{im_{Z}z}Z_{\mu}~{},$
(63) $\displaystyle{\bf X}_{\mu}(x,z)=\cos\theta_{W}A_{\mu}-\sin\theta_{W}{\bf
Z}_{\mu}(x,z)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}=\cos\theta_{W}A_{\mu}-\sin\theta_{W}e^{im_{Z}z}Z_{\mu}~{},$
(64)

where $\theta_{W}$ is the Weinberg angle.

The nonmass eigenstates ${\bf W}_{\alpha}^{a}(a=1,2,3)$ in ${\bf T}\cdot{\bf
W}_{\alpha}$ must have the same rest mass because of two reasons: Each ${\bf
W}_{\alpha}^{a}$ plays the quite equal role in gauge field ${\bf T}\cdot{\bf
W}_{\alpha}$; The model must be $SU(2)$ gauge invariant. Consequently

$\displaystyle m_{{\bf W}^{1}_{\mu}}=m_{{\bf W}^{2}_{\mu}}=m_{{\bf
W}^{3}_{\mu}}~{}.$ (65)

In Ref.[5], it is indicated that the rest mass squared of nonmass eigenstate
of vector fields can be calculated, that is

$\displaystyle m_{{\bf V}_{\mu}}^{2}=\sum_{j=1}^{n}a_{j}a^{*}_{j}m_{j}^{2}$
(66)

is right if and only if

$\displaystyle{\bf V}_{\mu}=\sum_{j=1}^{n}a_{j}[{\bf
V}_{\mu}]_{j}=\sum_{j=1}^{n}a_{j}e^{im_{j}z}[V_{\mu}]_{j}~{}.$ (67)

Then one can easily obtain the rest masses of ${\bf W}^{1}_{\mu}$ and ${\bf
W}^{2}_{\mu}$ from (60) and (61) respectively

$\displaystyle m_{{\bf W}^{1}_{\mu}}^{2}=m_{{\bf
W}^{2}_{\mu}}^{2}=m_{W}^{2}~{},$ (68)

also get the rest mass of ${\bf W}^{3}_{\mu}$ from (63)

$\displaystyle m_{{\bf W}^{3}_{\mu}}^{2}=m_{Z}^{2}(\cos\theta_{W})^{2}~{}.$
(69)

Therefore combining (65), (68) and (69), I obtain the following relation

$\displaystyle m_{W}=m_{Z}\cos\theta_{W}~{}.$ (70)

In the total Lagrangian (44), the kinetic terms of Fermions and the
interaction terms will be discussed carefully in my forthcoming paper [6], in
this letter I only discuss the self-couplings of $SU(2)\times U(1)$ gauge
fields, namely the terms ${\cal L}_{U(1)}+{\cal L}_{SU(2)}$ in (44). It has
been pointed out that the gauge bosons in my model merely propagate along the
$z=0$ brane, therefore ${\bf W}_{z}(x,z)=0,~{}{\bf X}_{z}(x,z)=0$. In this
case, substituting (64) into (55) yields

$\displaystyle{{\bf E}}_{z\mu}=-im_{Z}\sin\theta_{W}e^{im_{Z}z}Z_{\mu}~{}.$
(71)

Substituting (60), (61) and (63) into (56) respectively, I obtain

$\displaystyle{{\bf
F}}_{z\mu}^{1}=\frac{i}{\sqrt{2}}m_{W}e^{im_{W}z}\left({W}^{+}_{\mu}+{W}^{-}_{\mu}\right)~{},$
(72) $\displaystyle{{\bf
F}}_{z\mu}^{2}=-\frac{1}{\sqrt{2}}m_{W}e^{im_{W}z}\left({W}^{+}_{\mu}-{W}^{-}_{\mu}\right)~{},$
(73) $\displaystyle{{\bf
F}}_{z\mu}^{3}=im_{Z}\cos\theta_{W}e^{im_{Z}z}Z_{\mu}~{}.$ (74)

Substituting (71), (72), (73) and (74) into ${\cal L}_{U(1)}+{\cal
L}_{SU(2)}$, I find that the self-coupling terms of $SU(2)\times U(1)$ gauge
fields become

$\displaystyle~{}~{}~{}~{}{\cal L}_{U(1)}+{\cal L}_{SU(2)}$
$\displaystyle=-\frac{1}{4}{{\bf E}}_{\alpha\beta}{\bf
E}^{\alpha\beta}-\frac{1}{4}{{\bf F}}_{\alpha\beta}\cdot{{\bf
F}}^{\alpha\beta}$ $\displaystyle=-\frac{1}{4}{{\bf E}}_{\mu\nu}{\bf
E}^{\mu\nu}-\frac{1}{4}{{\bf F}}_{\mu\nu}\cdot{{\bf
F}}^{\mu\nu}-\frac{1}{2}\left({{\bf E}}_{z\mu}{\bf E}^{z\mu}+{{\bf
F}}_{z\mu}\cdot{{\bf F}}^{z\mu}\right)$ $\displaystyle=-\frac{1}{4}{{\bf
E}}_{\mu\nu}{\bf E}^{\mu\nu}-\frac{1}{4}{{\bf F}}_{\mu\nu}\cdot{{\bf
F}}^{\mu\nu}-\frac{1}{2}\left({{\bf E}}_{z\mu}{\bf E}^{z\mu}+{\bf
F}^{1}_{z\mu}{{\bf F}}^{1z\mu}+{\bf F}^{2}_{z\mu}{{\bf F}}^{2z\mu}+{\bf
F}^{3}_{z\mu}{{\bf F}}^{3z\mu}\right)$ $\displaystyle=-\frac{1}{4}{{\bf
E}}_{\mu\nu}{\bf E}^{\mu\nu}-\frac{1}{4}{{\bf F}}_{\mu\nu}\cdot{{\bf
F}}^{\mu\nu}+\frac{1}{2}m_{Z}^{2}e^{2im_{Z}z}Z_{\mu}Z^{\mu}+m_{W}^{2}e^{2im_{W}z}W_{\mu}^{+}W^{\mu-}~{}.$
(75)

It is indicated that all the mass eigenstates coupled by the gravitation, the
electromagnetic field and the gluon fields are living in the $z=0$ brane.
Hence, expressed by the mass eigenstates in the $z=0$ brane, the self-coupling
terms of $SU(2)\times U(1)$ gauge fields reduce to

$\displaystyle~{}~{}~{}~{}{\cal L}_{U(1),z=0}+{\cal L}_{SU(2),z=0}$
$\displaystyle=-\frac{1}{4}{E}_{\mu\nu}{E}^{\mu\nu}-\frac{1}{4}{F}_{\mu\nu}\cdot{F}^{\mu\nu}+\frac{1}{2}m_{Z}^{2}Z_{\mu}Z^{\mu}+m_{W}^{2}W_{\mu}^{+}W^{\mu-}~{},$
(76)

where
${F}_{\mu\nu}\equiv\\{{F}^{1}_{\mu\nu},{F}^{2}_{\mu\nu},{F}^{3}_{\mu\nu}\\}$,
and ${E}_{\mu\nu}$ is formulated by

$\displaystyle{E}_{\mu\nu}(x)=\partial_{\mu}{X}_{\nu}(x)-\partial_{\nu}{X}_{\mu}(x)~{},$
(77)

and ${F}_{\mu\nu}$ is given by

$\displaystyle{F}_{\mu\nu}(x)=\partial_{\mu}{W}_{\nu}(x)-\partial_{\nu}{W}_{\mu}(x)+g_{1}{W}_{\mu}(x)\times{W}_{\nu}(x)~{},$
(78)

in which
${W}_{\mu}(x)\equiv\\{{W}^{1}_{\mu}(x),{W}^{2}_{\mu}(x),{W}^{3}_{\mu}(x)\\}$.
The four-dimensional fields ${X}_{\mu}(x)$ and ${W}_{\mu}(x)$ are composed of
mass eigenstates which are constrained in the $z=0$ brane. From (60), (61),
(63) and (64), one can easily obtain the expressions of them in the following

$\displaystyle{W}^{1}_{\mu}(x)=\frac{1}{\sqrt{2}}\left({W}^{+}_{\mu}(x)+{W}^{-}_{\mu}(x)\right)~{},$
(79)
$\displaystyle{W}^{2}_{\mu}(x)=\frac{i}{\sqrt{2}}\left({W}^{+}_{\mu}(x)-{W}^{-}_{\mu}(x)\right)~{},$
(80)
$\displaystyle{W}^{3}_{\mu}(x)=\sin\theta_{W}A_{\mu}(x)+\cos\theta_{W}Z_{\mu}(x)~{},$
(81)
$\displaystyle{X}_{\mu}(x)=\cos\theta_{W}A_{\mu}(x)-\sin\theta_{W}Z_{\mu}(x)~{}.$
(82)

Obviously they have the same forms as the definitions of gauge bosons of
$SU(2)\times U(1)$ gauge fields in the Glashow-Salam-Weinberg model [7]. The
fields ${W}^{+}_{\mu}(x)$, ${W}^{-}_{\mu}(x)$, $Z_{\mu}(x)$ and $A_{\mu}(x)$
in above expressions are mass eigenstates that are constrained in the $z=0$
brane.

The $SU(2)\times U(1)$ gauge fields of the second kind merely couple to the
nonmass eigenstates, which are the nonmass eigenstates as well, hence cannot
be observed directly. The nonmass eigenstates of gauge fields are composed of
the mass eigenstates that are constrained in the $z=0$ brane. When I reexpress
the Lagrangian ${\cal L}_{U(1)}+{\cal L}_{SU(2)}$ by the mass eigenstates of
gauge bosons who live in four-dimensional spacetime, I find that the fields
${W}^{+}_{\mu}(x)$, ${W}^{-}_{\mu}(x)$ and $Z_{\mu}(x)$ can be treated as
massive gauge bosons from the four-dimensional point of view since their mass
terms automatically appear in the four-dimensional Lagrangian (76).

The mass terms of gauge bosons who are living in four-dimensional spacetime
aren’t inserted by hands, which is produced automatically. From five-
dimensional point of view, the $SU(2)\times U(1)$ gauge fields of the second
kind are massless, which are the usual gauge fields that we are very familiar
with, since there are no mass term in the total Lagrangian (44). Therefore,
the gauge fields in this model should be renormalizable.

let me explicitly explain this model again: The general equation of a nonmass
eigenstate of spin-$\frac{1}{2}$ fields is built on a five-dimensional
manifold, therefore the gauge fields of the second kinds, who couple to the
nonmass eigenstates of spin-$\frac{1}{2}$ fields, are constructed on a five-
dimensional manifold as well. But the nonmass eigenstates are composed of mass
eigenstates, which are physically observable. The mass eigenstates are merely
coupled by the electromagnetic field, the gravitation and the gluon fields,
who are living in the $z=0$ brane. Hence the initial momentum of the nonmass
eigenstates of spin-$\frac{1}{2}$ fields are along the $z=0$ brane, which
decides that the gauge fields of the second kind must propagate along the
$z=0$ spacetime as well. When the gauge fields of the second kind are
reexpressed by their mass eigenstates that living in the four-dimensional
spacetime, it is found that the gauge bosons who are mass eigenstates can be
treated as massive vector fields.

## References

  * [1] C. N. Yang and R. L. Mills, Phys. Rev. 96, 191 (1954).
  * [2] P. W. Higgs, Phys. Rev. Lett. 13, 508 (1964); P. W. Higgs, Phys. Rev. 145, 1156 (1966); F. Englert and R. Brout, Phys. Rev. Lett. 13, 321 (1964).
  * [3] S. L. Glashow, Nucl. Phys. 22, 579 (1960); J. Goldstone, A. Salam and S. Weinberg, Phys. Rev. 127, 965 (1962); S. Weinberg, Phys. Rev. Lett. 19, 1264 (1967); S. L. Glashow, J. Iliopoulos and L. Maiani, Phys. Rev. D 2, 1285 (1970).
  * [4] For recent reviews on neutrino masses and mixing angles, see Z. Z. Xing, Int. J. Mod. Phys. A 19, 1 (2004); R. D. McKeown and P. Vogel, Phys. Rept. 394, 315 (2004); M. C. Gonzalez-Garcia and Y. Nir, Rev. Mod. Phys. 75, 345 (2003); M.-C. Chen and K. T. Mahanthappa, Int. J. Mod. Phys. A 18, 5819 (2003).
  * [5] X.-B. Huang, “Nonmass Eigenstates of Boson and Fermion Fields”, [arXiv:hep-th/0906.2441].
  * [6] X.-B. Huang, “An Electroweak Model without Higgs Mechanism”, in preparation.
  * [7] W. Greiner and B. Müller, Gauge Theory of Weak Interactions, (3rd. edition), (Springer-Verlag, 2000).


# Neural network augmented wave-equation simulation

Ali Siahkoohi, Mathias Louboutin, and Felix J. Herrmann
School of Computational Science and Engineering
Georgia Institute of Technology
{alisk, mlouboutin3<EMAIL_ADDRESS>
###### Abstract

Accurate forward modeling is important for solving inverse problems. An
inaccurate wave-equation simulation, as a forward operator, will offset the
results obtained via inversion. In this work, we consider the case where we
deal with incomplete physics. One proxy of incomplete physics is an inaccurate
discretization of Laplacian in simulation of wave equation via finite-
difference method. We exploit intrinsic one-to-one similarities between
timestepping algorithm with Convolutional Neural Networks (CNNs), and propose
to intersperse CNNs between low-fidelity timesteps. Augmenting neural networks
with low-fidelity timestepping algorithms may allow us to take large timesteps
while limiting the numerical dispersion artifacts. While simulating the wave-
equation with low-fidelity timestepping algorithm, by correcting the wavefield
several time during propagation, we hope to limit the numerical dispersion
artifact introduced by a poor discretization of the Laplacian. As a proof of
concept, we demonstrate this principle by correcting for numerical dispersion
by keeping the velocity model fixed, and varying the source locations to
generate training and testing pairs for our supervised learning algorithm.

##  1 Introduction

In inverse problem, we heavily rely on having an accurate forward modeling
operators. Often, we can not afford being physically or numerically accurate.
In other words, being numerically inaccurate can be due to computationally
complexity of accurate methods, or incomplete knowledge of the underlying data
generation process. In either case, motivated by Ruthotto and Haber [1], we
propose to intersperse CNNs between timestepping algorithm for simulating
acoustic wave-equation. We mimic incomplete/inaccurate physics by simulating
wave equation with finite-difference method, while utilizing a poor (second-
order) discretization of Laplacian.

Conventional method for solving partial differential equations (PDEs), e.g.,
finite-difference and finite-element method, given enough computational
recourses, are able to simulate high-fidelity solutions to PDEs. On one hand,
as long as the Courant–Friedrichs–Lewy conditions for stability are satisfied,
finite-difference methods are able to compute solutions to PDE, regardless of
medium parameters, with arbitrary precision. On the other hand, finite-element
method requires careful meshing of the medium in order to carry out the
simulation. Another motivation behind this work is to exploit the fact that in
seismic, wave simulations are usually carried out for specific families of
velocity models and source/receiver distributions. We hope our proposed method
meets halfway between two mentioned extremes—i.e., being too generic (finite-
difference method) and being too problem specific (finite-element method).

There are several attempts to exploit learning methods in wave-equation
simulation. Raissi [2] approximates the solution to a nonlinear PDE with a
neural network. The neural network, given points on the computational grid as
input, computes the solution of PDE. Training data is obtained by computing
the solution of PDE on several points. Moseley et al. [3] completely ignore
the Laplacian and they solely rely on predicting the next timestep from the
previous two timesteps by learning the action of the spatially varying
velocity and Laplacian. While possible in principle, their approach needs to
train for long times to provide reasonable simulations on relatively simple
models. Siahkoohi et al. [4] instead of ignoring the physics, relies on low-
fidelity wave-equation simulation, and by exploiting transfer learning, they
utilize a single CNN to correct wavefield snapshots simulated on a “nearby”
velocity mode for numerical dispersion at any given timestep. In this work, we
extend ideas in Siahkoohi et al. [4] and propose using multiple CNNs,
interspersed between low-fidelity timesteps. Finally, Rizzuti et al. [5]
propose interspersing Krylov-subspace iterations and neural nets while
inverting the Helmholtz equation. They show improvement in convergence by
“propagating” an approximated wavefield, obtained from a limited number of
iterations, with the aid of a trained convolutional neural net. This technique
can be seen as the frequency-domain counterpart of our proposed method.

Our paper is organized as follows. First, we describe our approach in detail
by first, describing our formulation for learned wave simulation. Next, we
introduce our training objective function. Due to dependencies of CNN
parameters, we devised an training heuristic that we describe. Before
explaining our numerical experiments, we state used CNN architecture and
training details. Next, we describe our three numerical experiments we conduct
and discuss effectiveness of the proposed method.

##  2 Theory

We describe how we augment low-fidelity physics with learning techniques to
handle incomplete and/or inaccurate physics, where the low-fidelity physics is
modeled via finite-difference method with a poor discretization of the
Laplacian. To ensure accuracy, the temporal and spatial discretization in
high-fidelity wave-equation simulations have to be chosen very fine, typically
one to two orders of magnitude smaller than Nyquist sampling rate. As
mentioned earlier, we will utilize a poor discretization (only second order)
for the Laplacian to carry out low-fidelity wave-equation simulations, but the
scheme can be extended to other proxies of incomplete or inaccurate physics.

###  2.1 Simulations by timestepping

After discretization of the acoustic wave equation, a single timestep of of
scalar wavefields, simulated on $0\leq t\leq T$, can be written as below:

$\mathbf{u}_{j+1}=2\mathbf{u}_{j}-\mathbf{u}_{j-1}+\delta
t^{2}\mathbf{c}^{2}\mathbf{\Delta}\mathbf{u}_{j},\quad j=0,1,\ldots,N-1,$ (1)

where $\mathbf{u}_{j}$ is the high-fidelity scalar wavefield at
$j^{\text{th}}$ timestep, $\delta t$ is the temporal discretization
(timestep), $\mathbf{c}$ is the spatially varying velocity in the medium, and
$\mathbf{\Delta}$ is the high-order discretization of Laplacian. Similar to
Equation 1, the low-fidelity timestepping equation can be formulated as

$\mathbf{\bar{u}}_{j+1}=2\mathbf{\bar{u}}_{j}-\mathbf{\bar{u}}_{j-1}+\delta
T^{2}\mathbf{\bar{c}}^{2}\mathbf{\bar{\Delta}}\mathbf{\bar{u}}_{j},\quad
j=0,1,\ldots,M-1,$ (2)

where $\mathbf{\bar{u}}_{j}$ is the low-fidelity scalar wavefield, $\delta T$
is the coarse timestep, $\mathbf{\bar{c}}$ is the coarse spatially varying
velocity, and $\mathbf{\bar{\Delta}}$ is the coarse (only second second order)
discretized Laplacian.

Motivated by Ruthotto and Haber [1], we consider every timestep as a single
layer in a neural network, where the discretized Laplacian is a linear
operator, followed by the (nonlinear) action of the spatial varying velocity.
Moreover, noticing the additional terms in the Equation 1, each timestep is
similar to a residual block introduced by Szegedy et al. [6]. Figures 1a and
1b schematically indicate each timestep as a block, corresponding to high- and
low-fidelity discretization of wave equation. respectively. The similarity of
high- and low -fidelity timestepping method and CNNs can be perceived from
Figures 2a and low-fidelity-step, respectively, where red and yellow blocks
correspond to high- and low-fidelity timestepping equations, respectively.

(a)

(b)

Figure 1: Comparing a single low and high-fidelity timestep. a) High-fidelity
timestep. b) Low-fidelity timestep.

(a)

(b)

Figure 2: Comparing low and high-fidelity discretized wave equation
simulations. a) High-fidelity simulation. b) Low-fidelity simulation.

As it can be seen, high-fidelity simulation of the wave equation up to time
$t=T$ requires a lot of high-fidelity timesteps. On the other hand, Figure 2b
shows that the low-fidelity simulations can be done with much less low-
fidelity timesteps, due to course time sampling, which each timestep is
cheaper than the high-fidelity timesteps due to the coarse discretization of
Laplacian. Although computationally cheap, the low-fidelity wave-equation
simulations suffer from numerical dispersion artifacts.

###  2.2 Learned wave simulations

Depending on the domain of application, we can assume wave simulations are
typically carried out for specific families of velocity models and
source/receiver distributions. This motivates us to deploy a data-driven wave
simulation algorithm which is coupled with low-fidelity and cheap physics and
hope to recover high-fidelity wave simulations on a family of velocity models.
In our method, we propagate the coarse-grained wavefields according to
Equation 2 with a coarsened Laplacian. After $k$ timesteps, where $k$ is a
hyperparameter, we apply a correction with a CNN, $\mathcal{G}_{\theta_{i}}$,
parameterized by $\theta_{i}$, to the obtained wavefield at $j^{\text{th}}$
timestep and proceed with the timestepping. The proposed data-driven
timestepping wave simulation method is formalized in Equation 3.

$\mathbf{\bar{u}}_{j+1}=\begin{cases}\mathcal{G}_{\theta_{i}}\left(2\mathbf{\bar{u}}_{j}-\mathbf{\bar{u}}_{j-1}+\delta
T^{2}\mathbf{\bar{c}}^{2}\mathbf{\bar{\Delta}}\mathbf{\bar{u}}_{j}\right),\
i=\lfloor\frac{j}{k}\rfloor&\text{ if }j\equiv k-1\ (\bmod\ k),\\\
2\mathbf{\bar{u}}_{j}-\mathbf{\bar{u}}_{j-1}+\delta
T^{2}\mathbf{\bar{c}}^{2}\mathbf{\bar{\Delta}}\mathbf{\bar{u}}_{j}&\text{ else
}\end{cases}$ (3)

where $j=0,1,\ldots,M-1$. The schematic representation of Equation 3 is
illustrated in Figure 3. Yellow blocks represent low-fidelity timesteps (see
Equation 2 and Figure 1b) and blue blocks correspond to CNNs,
$\mathcal{G}_{\theta_{i}},\ i=0,1,\ldots,\lfloor\frac{M-1}{k}\rfloor$.

Figure 3: A schematic representation of the proposed method.

The CNNs correct for the effects of inaccurate physics—i.e., numerical
dispersion in our experiments, at every $k^{\text{th}}$ low-fidelity timestep.
In this work, the parameters, $\theta_{i}$, are not shared among the CNNs. We
have not explored the possibility of shared weights among the CNNs in
different stages of wave propagation.

Note that although parameters of the CNNs are not shared, they are not
independent—i.e., after $j^{\text{th}}$ timestep, CNN
$\mathcal{G}_{\theta_{i}},i=\lfloor\frac{j}{k}\rfloor$, corrects for errors in
the wavefield introduced by low-fidelity timestepping and imperfections
present in the output of $(i-1)^{\text{th}}$ CNN, which have been propagated
through timestepping. Therefore, a small perturbation in the parameters of a
CNN in the initial stages of neural network augmented timestepping causes
noticeable differences in the input of the CNNs in later stages of the wave
propagation. The described dependencies among the CNN parameters introduces
difficulties in optimizing the parameters of the CNNs. Below we describe our
heuristic for training the CNNs.

###  2.3 Training objective

During training, we train all the CNNs toward the high-fidelity solution of
wave-equation, at the corresponding timestep, obtained by solving Equation 1.
As it can be seen from Equation 3, after $j^{\text{th}}$ timestep, CNN
$\mathcal{G}_{\theta_{i}},i=\lfloor\frac{j}{k}\rfloor$, is tasked to correct
the effects of low-fidelity timestepping. During training, $i^{\text{th}}$ CNN
maps its input, $2\mathbf{\bar{u}}_{j}-\mathbf{\bar{u}}_{j-1}+\delta
T^{2}\mathbf{\bar{c}}^{2}\mathbf{\bar{\Delta}}\mathbf{\bar{u}}_{j}$, to
$\mathbf{u}_{j+1}$, result of $j^{\text{th}}$ timestep using high-fidelity
timestepping, obtained by Equation 1.

Define function $\bar{F}_{k}(.)$ as the action of $k$ low-fidelity
timesteps—i.e., $\bar{F}_{k}(.)$ represents $k$ consecutive low-fidelity time
stepping blocks, depicted in Figure 1b. Clearly, $\bar{F}_{k}$ is only a
function of $k$, $\delta T$, $\mathbf{\bar{c}}$, and $\mathbf{\bar{\Delta}}$.
Using the defined notation, we can write the input to $i^{\text{th}}$ CNN,
$\mathbf{\hat{u}}_{i},\ i=0,1,\ldots,\lfloor\frac{M-1}{k}\rfloor$, as follows:

$\displaystyle\
\mathbf{\hat{u}}_{i}=\bar{F}_{k}(\mathcal{G}_{\theta_{i-1}}(\mathbf{\hat{u}}_{i-1})),\quad
i=1,2,\ldots,\lfloor\frac{M-1}{k}\rfloor,$ (4) $\displaystyle\
\mathbf{\hat{u}}_{0}=\bar{F}_{k}(\mathbf{q}),$

where $\mathbf{q}$ is the source. Also let $\mathbf{u}_{\tau_{i}}$ denote the
wavefield obtained at $j^{\text{th}}$ timestep of high-fidelity timestepping
(Equation 1), where $\tau_{i}=j=(k+1)i-1$.

The input-output pair of the $i^{\text{th}}$ CNN is
$(\mathbf{\hat{u}}_{i},\mathbf{u}_{\tau_{i}})$. We can generated multiple
training pairs for CNNs by simulating
$(\mathbf{\hat{u}}_{i},\mathbf{u}_{\tau_{i}})$ pairs, for various velocity
models and source locations. Assume we have $n$ pairs of of training data for
CNNs, namely, $(\mathbf{\hat{u}}_{i}^{(p)},\mathbf{u}_{\tau_{i}}^{(p)}),\
p=0,1,\ldots,n-1$. The objective function for optimizing $i^{\text{th}}$ CNN
can be written as follows:

$\begin{split}\mathcal{L}_{i}=\frac{1}{n}\sum_{p=0}^{n-1}\left\|\mathcal{G}_{\theta_{i}}(\mathbf{\hat{u}}_{i}^{(p)})-\mathbf{u}_{\tau_{i}}^{(p)}\right\|_{1},\quad
i=0,1,\ldots,\lfloor\frac{M-1}{k}\rfloor.\\\ \end{split}$ (5)

In the past, in a similar attempt, we used Generative Adversarial Networks
[GANs, 7] to train a CNN in order to remove numerical dispersion from
wavefield snapshots [8, 4]. In this work we choose to use $\ell_{1}$ norm as
the misfit function based on two reasons. First, training GANs is
computationally expensive since it requires training an additional neural
network that discerns between high-fidelity wavefield snapshots and corrected
ones. The computational complexity of the proposed method in this work is
significantly higher than our previous attempts [8, 4], because it involves
training multiple CNNs. Based on the mentioned facts, for limiting the
computation time we chose to use $\ell_{1}$ misfit function. Second, motivated
by a numerical experiment performed by Hu et al. [9], $\ell_{1}$ norm misfit
function yields the second best results after the misfit function utilizing a
combination of GANs and $\ell_{1}$ norm misfit. In the next section, we
describe our heuristic for training the CNNs.

###  2.4 Training heuristic

To overcome complexities caused by dependencies between parameters of CNNs, we
optimize the objective functions $\mathcal{L}_{i}$ with a heuristic described
below. We minimize $\mathcal{L}_{i},\
i=0,1,\ldots,\lfloor\frac{M-1}{k}\rfloor$, with respect to $\theta_{i}$,
respectively. In other words, we minimize $\mathcal{L}_{i}$ with respect to
$\theta_{i}$, by keeping the rest of the parameters fixed. We keep updating
all the set of parameters, in a cyclic fashion—i.e., once we updated all the
parameters, $\mathcal{L}_{i},\ i=0,1,\ldots,\lfloor\frac{M-1}{k}\rfloor$, we
start over and update them again, in order, until a stopping criteria is
achieved. We will describe the stopping criteria used in our experiments
later.

We minimize objective functions the above objectives 5 with a variant of
Stochastic Gradient Descent known as the Adam optimizer [10] with momentum
parameter $\beta=0.9$ and a linearly decaying stepsize with initial value
$\mu=2\times 10^{-4}$ for both the generator and discriminator networks.
During each iteration of Adam, the gradient $\mathcal{L}_{i}$ is approximated
by a single randomly selected training pair. These pairs are selected without
replacement. Once all the training pairs have been selected, we start over by
randomly picking training pairs, without replacement from the entire training
set.

The optimization carries out for a predetermined number of total _iterations_
, where each iterations consists of drawing a random training pair, without
replacement, and updating parameters of a CNN. Additionally, while optimizing
$\theta_{i}$ by keeping the rest of the parameters fixed, before proceeding to
the next set of parameters, we carry out the optimization to update
$\theta_{i}$ for number of iterations, which we refer to it as _mini-
iterations_. Algorithm 2.4 indicates the steps for optimizing objective
functions 5.

1\. INPUT:
$MaxItr\quad$ // total number of iterations to carry out the optimization
$MaxMiniItr\quad$ // mini-iterations before proceeding to the next CNN
$\bar{F}_{k}(.)\quad$ // k consecutive low-fidelity time stepping blocks
$\mathbf{q}^{(p)}\ p=0,1,\ldots,n-1\quad$ \\\ sources corresponding to
different training pairs
$\mathbf{u}_{\tau_{i}}^{(p)},\ p=0,1,\ldots,n-1,\
i=0,1,\ldots,\lfloor\frac{M-1}{k}\rfloor\quad$ \\\ high-fidelity snapshots
${\theta^{0}_{i}},\ i=0,1,\ldots,\lfloor\frac{M-1}{k}\rfloor\quad$ // randomly
initialized parameters
2\. $ItrNum\leftarrow 0$
3\. FOR $i=0:\lfloor\frac{M-1}{k}\rfloor$ DO
4\. $\theta_{i}=\theta^{0}_{i}$
5\. FOR $p=0:n-1$ DO
6\. $\mathbf{\hat{u}}_{0}^{(p)}=\bar{F}_{k}(\mathbf{q}^{(p)})$
7\. WHILE $itrNum<MaxItr$ DO
8\. FOR $i=0:\lfloor\frac{M-1}{k}\rfloor$ DO
9\. IF $i>0$ DO
10\. FOR $p=0:n-1$ DO
11\.
$\mathbf{\hat{u}}_{i}^{(p)}=\bar{F}_{k}(\mathcal{G}_{\theta_{i-1}}(\mathbf{\hat{u}}_{i-1}^{(p)}))$
12\. FOR $miniItrNum=1:MaxMiniItr$ DO
13\. $p\leftarrow SampleWithoutReplacement(\left\\{0,1,\ldots,n-1\right\\})$
14\.
$\theta_{i}\leftarrow\mathop{\arg\min}_{\theta_{i}}\left\|\mathcal{G}_{\theta_{i}}(\mathbf{\hat{u}}_{i}^{(p)})-\mathbf{u}_{\tau_{i}}^{(p)}\right\|_{1}$
15\. $ItrNum\leftarrow ItrNum+1$
16\. RETURN ${\theta_{i}},\ i=0,1,\ldots,\lfloor\frac{M-1}{k}\rfloor$

Algorithm 1 Heuristic for optimizing CNNs $\mathcal{G}_{\theta_{i}},\
i=0,1,\ldots,\lfloor\frac{M-1}{k}\rfloor$.

###  2.5 CNN architecture

Motivated by our previous attempts for numerical dispersion removal from
wavefield snapshots [8, 4], we use the exact architecture provided by Johnson
et al. [11], which includes Residual Blocks, the main building block of
ResNets, introduced by He et al. [12], for all the CNNs
$\mathcal{G}_{\theta_{i}},\ i=0,1,\ldots,\lfloor\frac{M-1}{k}\rfloor$.

###  2.6 Training details and implementation

While CNNs are known to generalize well—i.e., maintain the quality of
performance when applied to unseen data, they can only be successfully applied
to a data set drawn from the same distribution as the training data. Because
of the Earth’s heterogeneity and complex geological structures present in
realistic-looking models, training a neural network that can generalize well
when applied to another velocity model can become challenging. While we have
successfully demonstrated that transfer learning [13] can be used in
situations where the neural network is initially trained on data from a
proximal survey [4], we chose in this contribution, as a proof of concept, to
keep the velocity model fixed, and vary the source locations to generate
different training/testing pairs.

We use the Marmousi velocity model and out of $401$ available shot locations
with $7.5$ $\mathrm{m}$ spacing, we allocate half of the shot locations to
training and use the rest of the shot locations to generate testing pairs, for
evaluation purposes. The maximum simulation time in our experiments in $1.1$
$\mathrm{s}$.

We designed and implemented our deep architectures in
TensorFlow111https://www.tensorflow.org/. To carry out our wave-equation
simulations with finite differences, we used
Devito222https://www.devitoproject.org/ [14, 15]. We used the functionality of
Operator Discretization Library333https://odlgroup.github.io/odl/ to wrap
Devito operators into a TensorFlow layers. Our implementation can be found on
GitHub444https://github.com/alisiahkoohi/NN-augmented-wave-sim.

We ran our algorithm on Amazon Web Services’ g3.4xlarge instance, where we
optimize the CNN parameters on a NVIDIA Tesla M60 GPU and Devito utilizes $16$
CPU cores to perform finite-difference wave-equation simulations. Initially,
we simulate the high-fidelity training wavefield snapshots,
$\mathbf{u}_{\tau_{i}}^{(p)},\ p=0,1,\ldots,n-1,\
i=0,1,\ldots,\lfloor\frac{M-1}{k}\rfloor$, only once, in the beginning, and
store them. In order to limit CPU-GPU communication, before utilizing the GPU
to to update ${\theta_{i}},\ i=0,1,\ldots,\lfloor\frac{M-1}{k}\rfloor$, we
generate the input to $i^{\text{th}}$ CNN, $\mathbf{\hat{u}}_{i}^{(p)},\
p=0,1,\ldots,n-1$ all at once, and store them. Afterwards, $i^{\text{th}}$ CNN
can be (re)trained using the stored input/output wavefield snapshot pairs for
several mini-iterations.

##  3 Numerical experiments

We want to indicate that neural networks, when augmented with inaccurate
physics, e.g., a poor discretization of Laplacian, are able to approximate the
wavefields obtained by an accurate approximation to wave equation. To
demonstrate this, we conduct three numerical experiments in which we keep the
velocity model fixed, and vary the source locations to generate different
training/testing pairs. The experiments differ in the number of CCNs used
throughout learned wave propagation. We use three, five, and ten CNNs while
keeping the total number of iterations fixed. This implies that an experiment
with more CNNs, optimizes each CNN with a smaller number of iterations per
CNN, because, $\text{iterations per CNN}\times\text{number of
CNNs}=\text{total number of iterations}$.

A neural network augmented wave simulator with $n_{1}$ CNNs needs more
training iterations and possibly more training data to perform equally as well
as a neural network augmented wave simulator with $n_{2}$ CNNs, when
$n_{1}>n_{2}$. For a fixed number of total iterations, iterations per CNN is
inversely proportional to number of CNNs utilized. Therefore, the first
$n_{2}$ CNNs in the neural network augmented wave simulator with $n_{1}$ CNNs
will perform worse than the CNNs in the wave propagator with $n_{2}$ CNNs.
Consequently, the error accumulated by the poor performance of first $n_{2}$
CNNs, combined with artifacts introduced by low-fidelity timestepping makes
the matters worse for the later CNNs in the more complex learned wave
propagator. Therefore, in our experiments, since the total number of
iterations is fixed, we expect to see the quality of dispersion removal
degrade as the number of CNNs increase in a learned wave propagator. Table 1
summarizes the total number of iterations, iterations per CNN, training pairs,
training time, and number of tunable parameters for the three different
experiments.

CNNs | Iterations | Iterations per CNN | Pairs per CNN | Time | Param. count
---|---|---|---|---|---
$3$ | $100500$ | $33500$ | $201$ |  $17.99$ hours | $34150272$
$5$ | $100500$ | $20100$ | $201$ |  $19.79$ hours | $56917120$
$10$ | $100500$ | $10050$ | $201$ |  $49.24$ hours | $113834240$
Table 1: Summary of details in the three neural network augmented wave-
equation simulation experiments.

As described earlier and presented in Table 1, the MaxItr variable used in the
While condition in line $8$ of Algorithm 2.4 is set to $500$ for all our
experiments. Figures 4 $-$ 6 demonstrate the values of objective function
presented in Equation 5 in orange, and the wavefield correction signal-to-
noise ratio (SNR) in blue, evaluated on testing data pairs during training,
for experiments with three, five, and ten CNNs, respectively. Note that the
SNR curves have not been used to determine when to stop training and they are
only depicted for demonstration purposes.

Figures 4a, 4c, and 4e show the wavefield correction SNR for first, second,
and third CNN, respectively, in the neural network augmented wave simulator
that includes three CNNs. Similarly, Figures 4b, 4d, and 4f depict the
training objective values throughout training for first, second, and third
CNNs. As it can be seen from objective function curves, the raining heuristic
has been effective and the objective function values have a decreasing trend.
Note that CNNs has been trained for $33500$ iterations, on average, with a
total of $100500$ iterations. Several equispaced spikes can be noticed on the
objective function value curves. For instance, see the objective value
function curve of the third CNN, in Figure 4f, at $6030,8040,10050,\text{and
}12060$ iterations. The mini-batch we use in this experiment is $10$. Those
spikes occur in moments in training when we have started retraining the third
CNN, after updating the first and second CNNs. As discussed before, a change
in the parameters of the CNNs preceding a CNN causes changes in the input of
the later CNN, and consequently the objective function becomes large when
starting to retrain the CNNs in later stages again.

Similar objective function value and SNR curves for other two neural network
augmented wave propagators, utilizing five and ten CNNs, can be found in
Figures 5 and 6, respectively. First column in Figures 5 and 6, from top to
bottom, indicate SNR of wavefield correction obtained by the first to the last
CNN, evaluated on testing pairs while training, respectively. The second
column of Figures 5 and 6 indicate the objective function value curves
throughout optimization of Equation 5 for training neural network augmented
wave propagators, utilizing five and ten CNNs, respectively. In both columns,
from top to bottom, the objective function value curves correspond to CNNs
from beginning to the end of the learned wave propagators, in order.

We make two main observations from Figures 5 and 6. First, the objective
function values indicate overall decreases, validating effectiveness of the
introduced heuristic. Also, the spikes on the objective function value curve
can be seen, which are correlated with the stages in training when Algorithm
2.4 revisits a CNN after updating the rest of the CNN parameters. As explained
before, spikes are caused by change in parameters of preceding CNNs to a CNN,
which in turn alters the input training wavefields of the CNN. Second, due to
decrease in iterations per CNN as the number of CNNs increases, the SNR curves
converge to a lower value when the number of CNNs increase.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4: Neural network augmented wave simulation with three CNNs. First
column from top to bottom: SNR curves, evaluated on testing pairs during
training, for a) the first to e) the last CNN, in order. Second column from
top to bottom: training objective function value curves, evaluated on training
pairs, for b) the first to the f) last CNN, in order.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

Figure 5: Neural network augmented wave simulation with five CNNs. First
column from top to bottom: SNR curves, evaluated on testing pairs during
training, for a) the first to i) the last CNN, in order. Second column from
top to bottom: training objective function value curves, evaluated on training
pairs, for b) the first to the j) last CNN, in order.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

(q)

(r)

(s)

(t)

Figure 6: Neural network augmented wave simulation with ten CNNs. First column
from top to bottom: SNR curves, evaluated on testing pairs during training,
for a) the first to s) the last CNN, in order. Second column from top to
bottom: training objective function value curves, evaluated on training pairs,
for b) the first to the t) last CNN, in order.

Next, we will demonstrate the the corrected wavefields in three conducted
experiments evaluated over one testing shot location. For each experiment, we
show the high-fidelity wavefield snapshots, $\mathbf{u}_{\tau_{i}},\
i=0,1,\ldots,\lfloor\frac{M-1}{k}\rfloor$, where $i$ iterates over the CNNs,
numerically dispersed low-fidelity wavefields, and the corrected wavefield
snapshots by the CNNs. To evaluate the performance of each correction, we also
depict the correction error—i.e., difference between the high-fidelity and
corrected wavefield snapshots. Figure 7 shows the mentioned wavefield
snapshots for the neural network augmented wave simulator with three CNNs.
First column shows the high-fidelity wavefields by solving Equation 1, second
column depicts low-fidelity simulations by solving Equation 2, third column
indicates the result of neural network augmented wavefield simulations, and
the fourth column is the learned wave simulation error—i.e., difference
between the first and last column in Figure 7. Similarly, Figures 8 $-$ 10
show the high- and low- fidelity and learned wave simulation wavefield
snapshots in the first three columns, in order, for the neural network
augmented wave simulator with five and ten CNNs, respectively.

As expected because of the reasons stated before, we observe that the quality
of neural network augmented wave-equation simulation degrades as the number of
CNNs increases. On the other hands, the high quality of learned wave
simulation with few CNNs (see Figure 7) suggests the quality of the simulation
with more CNNs might be improved by increasing the number of iterations. As it
can be seen in the last column of Figures 9 and 10, the learned wave
simulation with ten CNNs has the lowest quality. It can be seen that the
learned wave simulation has the least accuracy in direct wave, which happen to
be the events with largest amplitudes. Also, it appears that most of the
numerical dispersion has been removed, the phase has been recovered, and
residual is mostly amplitude differences.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Figure 7: Neural network augmented wave simulation with three CNNs. First
column from top to bottom: a, e, i) high-fidelity wavefield snapshots, in
order. Second column from top to bottom: b, f, j) low-fidelity wavefield
snapshots simulated by solving Equation 2 with the same simulation time as
high-fidelity wavefields, in order. Third column from top to bottom: c, g, k)
result of neural network augmented wave-equation simulation. Output of the
first, second, and the last CNN, in order. Fourth column from top to bottom:
d, h, l) difference between first and third column, in order.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

(q)

(r)

(s)

(t)

Figure 8: Neural network augmented wave simulation with five CNNs. First
column from top to bottom: a) to q) high-fidelity wavefield snapshots, in
order. Second column from top to bottom: b) to r) low-fidelity wavefield
snapshots simulated by solving Equation 2 with the same simulation time as
high-fidelity wavefields, in order. Third column from top to bottom: c) to s)
result of neural network augmented wave-equation simulation. Output of the
first to the last CNN, in order. Fourth column from top to bottom: d) to t)
difference between first and third column, in order.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

(q)

(r)

(s)

(t)

Figure 9: Neural network augmented wave simulation with ten CNNs, first part.
First column from top to bottom: a) to q) high-fidelity wavefield snapshots,
in order. Second column from top to bottom: b) to r) low-fidelity wavefield
snapshots simulated by solving Equation 2 with the same simulation time as
high-fidelity wavefields, in order. Third column from top to bottom: c) to s)
result of neural network augmented wave-equation simulation. Output of the
first to the fifth CNN, in order. Fourth column from top to bottom: d) to t)
difference between first and third column, in order.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

(q)

(r)

(s)

(t)

Figure 10: Neural network augmented wave simulation with ten CNNs, second
part. First column from top to bottom: a) to q) high-fidelity wavefield
snapshots, in order. Second column from top to bottom: b) to r) low-fidelity
wavefield snapshots simulated by solving Equation 2 with the same simulation
time as high-fidelity wavefields, in order. Third column from top to bottom:
c) to s) result of neural network augmented wave-equation simulation. Output
of the sixth to the last CNN, in order. Fourth column from top to bottom: d)
to t) difference between first and third column, in order.

###  3.1 Performance comparison: Single CNN low-to-high-fidelity mapping

In order to evaluate the effectiveness of the proposed method we also train a
single CNN similar to our previous attempt to remove numerical dispersion from
wavefield snapshots [8, 4] and compare the result of numerical dispersion
removal with the proposed method. To be more precise, for each presented
neural network augmented wave-equation simulation experiment, where we use
three, five, and ten CNNs, we train a single CNN, $\mathcal{G}_{\theta}$, with
the same architecture as the architecture used in learned wave propagators, in
order to remove numerical dispersion from all the low-fidelity wavefield
snapshots simulated by solving Equation 2 for $j\equiv k-1\ (\bmod\ k)$, on
training shot locations. Likewise to previous examples, here we also use half
of the available shot locations to simulate training pairs, and the rest is
used to evaluate the performance of the trained CNN. The input to
$\mathcal{G}_{\theta}$ during training can be written as follows (compare with
Equation 4):

$\displaystyle\
\mathbf{\tilde{u}}_{i}=\bar{F}_{k}(\mathbf{\tilde{u}}_{i-1}),\quad
i=1,2,\ldots,\lfloor\frac{M-1}{k}\rfloor,$ (6) $\displaystyle\
\mathbf{\tilde{u}}_{0}=\bar{F}_{k}(\mathbf{q}),$

The desired output for the mentioned CNN is the high-fidelity wavefield
snapshots simulated on training shot locations, $\mathbf{u}_{\tau_{i}}^{(p)},\
p=0,1,\ldots,n-1,\ i=0,1,\ldots,\lfloor\frac{M-1}{k}\rfloor$. The objective
function for the mentioned CNN can be represented as follows:

$\begin{split}\mathcal{L}=\frac{1}{n(\lfloor\frac{M-1}{k}\rfloor+1)}\sum_{i=0}^{\lfloor\frac{M-1}{k}\rfloor}\sum_{p=0}^{n-1}\left\|\mathcal{G}_{\theta}(\mathbf{\tilde{u}}_{i}^{(p)})-\mathbf{u}_{\tau_{i}}^{(p)}\right\|_{1}.\\\
\end{split}$ (7)

We minimize objective function 7 over $\theta$ with Adam optimizes, using the
same maximum number of iterations as before, this time by combining all the
training pairs associated with different CNNs in the learned wavefield
simulation example. As mentioned before, in order to compare with the proposed
method, we minimize objective function 7 over three different set of input-
output pairs, each corresponding to our presented experiments with varying
number of CNNs. Table 2 summarizes the total number of iterations, training
pairs, training time, and number of tunable parameters for three different
cases, which differ in number of timesteps which we choose to correct the
numerical dispersion. This selected timesteps are associated with the
timesteps that the CNNs operated on, in our three previous examples.

# Timesteps to correct | Iterations | Pairs per CNN | Time | Param. count
---|---|---|---|---
$3$ | $100500$ | $603$ |  $13.85$ hours | $11383424$
$5$ | $100500$ | $1005$ |  $13.96$ hours | $11383424$
$10$ | $100500$ | $2010$ |  $14.56$ hours | $11383424$
Table 2: Summary of details in three conducted experiments

The slight difference in runtime among three different cases provided in Table
2 is partly due to different number of training pairs needed to be generated.
Figure 11 depicts the wavefield snapshot correction SNR curves, evaluated on
testing pairs while training, and the value of objective function 7, in single
CNN low-to-high-fidelity mapping experiment, as a function of number of
iterations. Figures 11a, 11c, and 11e show the SNR curves, when we trained the
CNN on wavefield snapshots correspond to three, five, and ten, timesteps,
respectively. Similarly, Figures 11b, 11d, and 11f depict the training
objective function value (Equation 7), when the CNN is trained on wavefield
snapshots correspond to three, five, and ten, timesteps, respectively. SNR
curves depicted in first column of Figure 11 show the evolution of wavefield
correction SNR evaluated on randomly selected testing wavefield wavefields
form the wavefield snapshots combined from different timesteps. Therefore,
Figures 11a, 11c, and 11e indicate that the three different CNNs converge to a
wavefield correction SNR around $20$ dB, regardless of number of timesteps
they are correcting for. Although this does not suggest that the performance
will stay the same as we increase the number of timesteps needed to be
corrected. By comparing Figure 11a with first column of Figure 4 (SNR curves
for neural network augmented wave-equation simulation with three CNNs), we
observe that, on average the two methods are performing equally well.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 11: Single CNN low-to-high-fidelity mapping. a) wavefield snapshot
correction SNR curve and b) objective function value curve when CNN is trained
on wavefield snapshot pairs corresponding to learned wavefield simulation with
three CNNs. c) wavefield snapshot correction SNR curve and d) objective
function value curve when CNN is trained on wavefield snapshot pairs
corresponding to learned wavefield simulation with five CNNs. e) wavefield
snapshot correction SNR curve and f) objective function value curve when CNN
is trained on wavefield snapshot pairs corresponding to learned wavefield
simulation with ten CNNs.

Finally, we will show the wavefield corrected by the single CNN low-to-high-
fidelity mapping method, for comparison with our proposed method. Figures 12
and 13 indicate the corrected wavefields, for cases where three and five
timesteps need to be corrected. Figures 14 and 15 demonstrate the corrected
wavefields, for the case where ten timesteps need to be corrected, in two
parts. In Figures 12 $-$ 15, the first, second, third, and fourth columns
depict the high-fidelity and low fidelity wavefield snapshots, corrected low-
fidelity wavefield snapshots, and the error in numerical dispersion removal,
respectively. In each columns, from top to bottom, the simulation time
increases.

For comparison between our proposed method and single CNN low-to-high-fidelity
mapping, compare Figures 7 with 12, 8 with 13, 9 with 14, 10 with 15. As it
can be seen, the single CNN low-to-high-fidelity mapping method maintains the
quality of its performance when the number of timesteps that need to be
corrected increases. On the other hand, as the number of CNNs in neural
network augmented wave-equation simulation increases, the performance drops,
by keeping the maximum number of iterations fixed. Also, by comparing Tables 1
and 2, we observe that the training time needed for single CNN low-to-high-
fidelity mapping, when number of timesteps needed to be corrected increases,
for fixed number of maximum iterations, grows very slowly compared to the
training time required for neural network augmented wave-equation simulation.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Figure 12: Single CNN low-to-high-fidelity mapping with three timesteps to be
corrected. First column from top to bottom: a, e, i) high-fidelity wavefield
snapshots, in order. Second column from top to bottom: b, f, j) low-fidelity
wavefield snapshots simulated by solving Equation 2 with the same simulation
time as high-fidelity wavefields, in order. Third column from top to bottom:
c, g, k) result of single CNN low-to-high-fidelity mapping. Fourth column from
top to bottom: d, h, l) difference between first and third column, in order.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

(q)

(r)

(s)

(t)

Figure 13: Single CNN low-to-high-fidelity mapping with five timesteps to be
corrected. First column from top to bottom: a) to q) high-fidelity wavefield
snapshots, in order. Second column from top to bottom: b) to r) low-fidelity
wavefield snapshots simulated by solving Equation 2 with the same simulation
time as high-fidelity wavefields, in order. Third column from top to bottom:
c) to s) result of single CNN low-to-high-fidelity mapping. Fourth column from
top to bottom: d) to t) difference between first and third column, in order.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

(q)

(r)

(s)

(t)

Figure 14: Single CNN low-to-high-fidelity mapping with ten timesteps to be
corrected, first part. First column from top to bottom: a) to q) high-fidelity
wavefield snapshots, in order. Second column from top to bottom: b) to r) low-
fidelity wavefield snapshots simulated by solving Equation 2 with the same
simulation time as high-fidelity wavefields, in order. Third column from top
to bottom: c) to s) result of single CNN low-to-high-fidelity mapping. Fourth
column from top to bottom: d) to t) difference between first and third column,
in order.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

(q)

(r)

(s)

(t)

Figure 15: Single CNN low-to-high-fidelity mapping with ten timesteps to be
corrected, second part. First column from top to bottom: a) to q) high-
fidelity wavefield snapshots, in order. Second column from top to bottom: b)
to r) low-fidelity wavefield snapshots simulated by solving Equation 2 with
the same simulation time as high-fidelity wavefields, in order. Third column
from top to bottom: c) to s) result of single CNN low-to-high-fidelity
mapping. Fourth column from top to bottom: d) to t) difference between first
and third column, in order.

##  4 Conclusions

Our numerical experiments demonstrate that, given suitable training data, the
well-trained neural network augmented wave-equation simulator is capable of
approximating wavefield snapshots simulated by high-fidelity simulation. In
this work, as a proxy of inaccurate physics, we simulate wave-equation with
finite-difference method, using a poor discretization of Laplacian. Although
not computationally favorable to high-fidelity wave simulation, we showed that
the learned wave simulator deals with inaccurate physics. An important
observation we made is that training time of the proposed method gets quickly
very long and to achieve high accuracy, it may not be possible to utilize too
many CNNs. On the other hand, training time required for the single CNN low-
to-high-fidelity mapping experiments, conducted for the sake of comparison,
grows very slowly as the number of timesteps to be corrected increases. In
future, we intend to initialize the CNN parameters in the proposed method with
the parameters of a CNN trained by the single CNN low-to-high-fidelity mapping
algorithm. The initialization may significantly reduce the training time
needed for the neural network augmented wave-equation simulation method, and
may give the chance to fine-tune the CNNs to the specific timestep that each
CNN is assigned to correct.

##  5 Acknowledgments

The authors thank Xiaowei Hu for his open-access
repository555https://github.com/xhujoy/CycleGAN-tensorflow on GitHub. Our
software implementation built on this work.

## References

  * Ruthotto and Haber [2018] Lars Ruthotto and Eldad Haber.  Deep neural networks motivated by partial differential equations.  _CoRR_ , abs/1804.04272, 2018.  URL http://arxiv.org/abs/1804.04272.
  * Raissi [2018] Maziar Raissi.  Deep hidden physics models: Deep learning of nonlinear partial differential equations.  _The Journal of Machine Learning Research_ , 19(1):932–955, 2018.
  * Moseley et al. [2018] Benjamin Moseley, Andrew Markham, and Tarje Nissen-Meyer.  Fast approximate simulation of seismic waves with deep learning.  _arXiv preprint arXiv:1807.06873_ , 2018.
  * Siahkoohi et al. [2019] Ali Siahkoohi, Mathias Louboutin, and Felix J. Herrmann.  The importance of transfer learning in seismic modeling and imaging.  2019\.  Submitted to GEOPHYSICS in February 2019.
  * Rizzuti et al. [2019] Gabrio Rizzuti, Ali Siahkoohi, and Felix J. Herrmann.  Learned iterative solvers for the helmholtz equation.  In _EAGE Annual Conference Proceedings_ , 06 2019.  doi:

10.3997/2214-4609.201901542. URL
https://slim.gatech.edu/Publications/Public/Conferences/EAGE/2019/rizzuti2019EAGElis/rizzuti2019EAGElis.pdf.
(EAGE, Copenhagen).

* Szegedy et al. [2017] Christian Szegedy, Sergey Ioffe, Vincent Vanhoucke, and Alexander A Alemi.  Inception-v4, Inception-ResNet and the Impact of Residual Connections on Learning.  In _Proceedings of the Thirty-First Association for the Advancement of Artificial Intelligence Conference on Artificial Intelligence (AAAI-17)_ , volume 4, pages 4278–4284, 2017.  URL http://aaai.org/ocs/index.php/AAAI/AAAI17/paper/view/14806.
* Goodfellow et al. [2014] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative Adversarial Nets.  _Advances in neural information processing systems_ , pages 2672–2680, 2014.
* Siahkoohi et al. [2018] Ali Siahkoohi, Mathias Louboutin, Rajiv Kumar, and Felix J. Herrmann.  Deep-convolutional neural networks in prestack seismic: Two exploratory examples.  _SEG Technical Program Expanded Abstracts 2018_ , pages 2196–2200, 2018.  doi: 10.1190/segam2018-2998599.1.  URL https://library.seg.org/doi/abs/10.1190/segam2018-2998599.1.
* Hu et al. [2019] Tao Hu, Zhizhong Han, Abhinav Shrivastava, and Matthias Zwicker.  Render4completion: Synthesizing multi-view depth maps for 3d shape completion.  _arXiv preprint arXiv:1904.08366_ , 2019.
* Kingma and Ba [2015] Diederik Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  In _International Conference on Learning Representations_ , 2015.
* Johnson et al. [2016] Justin Johnson, Alexandre Alahi, and Li Fei-Fei.  Perceptual Losses for Real-Time Style Transfer and Super-Resolution.  In _Computer Vision – European Conference on Computer Vision (ECCV) 2016_ , pages 694–711. Springer International Publishing, 2016.  doi: 10.1007/978-3-319-46475-6_43.  URL https://link.springer.com/chapter/10.1007%2F978-3-319-46475-6_43.
* He et al. [2016] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep Residual Learning for Image Recognition.  In _The IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 770–778, June 2016.  doi: 10.1109/CVPR.2016.90.  URL https://ieeexplore.ieee.org/document/7780459.
* Yosinski et al. [2014] Jason Yosinski, Jeff Clune, Yoshua Bengio, and Hod Lipson.  How transferable are features in deep neural networks?  In _Advances in neural information processing systems_ , pages 3320–3328, 2014.
* Louboutin et al. [2018] M. Louboutin, M. Lange, F. Luporini, N. Kukreja, P. A. Witte, F. J. Herrmann, P. Velesko, and G. J. Gorman.  Devito: an embedded domain-specific language for finite differences and geophysical exploration.  _CoRR_ , abs/1808.01995, Aug 2018.  URL https://arxiv.org/abs/1808.01995.
* Luporini et al. [2018] F. Luporini, M. Lange, M. Louboutin, N. Kukreja, J. Hückelheim, C. Yount, P. Witte, P. H. J. Kelly, G. J. Gorman, and F. J. Herrmann.  Architecture and performance of devito, a system for automated stencil computation.  _CoRR_ , abs/1807.03032, jul 2018.  URL http://arxiv.org/abs/1807.03032.