#  Physics-guided Full Waveform Inversion using Encoder-Solver Convolutional
Neural Networks

††The Israel Science Foundation supported this research grant no. 656/23 and
by the Lynn and William Frankel Center for Computer Science at BGU. MG is
supported by the Kreitman High-Tech scholarship. Matan M. Goren, Eran Treister
††Computer Science Department, Ben-Gurion University of the Negev, Beer Sheva,
Israel.
<EMAIL_ADDRESS>erant@cs.bgu.ac.il)2

###### Abstract

Full Waveform Inversion (FWI) is an inverse problem for estimating the wave
velocity distribution in a given domain, based on observed data on the
boundaries. The inversion is computationally demanding because we are required
to solve multiple forward problems, either in time or frequency domains, to
simulate data that are then iteratively fitted to the observed data. We
consider FWI in the frequency domain, where the Helmholtz equation is used as
a forward model, and its repeated solution is the main computational
bottleneck of the inversion process. To ease this cost, we integrate a
learning process of an encoder-solver preconditioner that is based on
convolutional neural networks (CNNs). The encoder-solver is trained to
effectively precondition the discretized Helmholtz operator given velocity
medium parameters. Then, by re-training the CNN between the iterations of the
optimization process, the encoder-solver is adapted to the iteratively
evolving velocity medium as part of the inversion. Without retraining, the
performance of the solver deteriorates as the medium changes. Using our light
retraining procedures, we obtain the forward simulations effectively
throughout the process. We demonstrate our approach to solving FWI problems
using 2D geophysical models with high-frequency data.

Keywords: Full Waveform Inversion, Helmholtz equation, Multigrid,
Preconditioner, Convolutional Neural Networks, U-Net, Encoder-Solver

##  1 Introduction

Full waveform inversion (FWI) is an inverse problem where we estimate the wave
propagation velocity based on recorded seismic data on the boundaries. The
solution of FWI was originally developed for geophysics purposes such as sub-
surface mapping and seismic exploration of oil and gas reservoirs, but was
shown to be applicable for other imaging tasks including brain imaging [29],
ultrasound imaging [47] and optical diffraction tomography [14]. The problem
is formulated as an optimization problem and is solved by a descent algorithm
which iteratively updates the underlying velocity model to reduce the misfit
function, measuring the difference between the recorded seismic data and the
simulated waveforms. At each iteration, the waveform simulation is obtained
for multiple sources for a given velocity model. This is done by either time
integrating through the wave equation (in the time domain) or solving the
Helmholtz equation (in the frequency domain) for multiple frequencies. That
is, in the frequency domain, which is the focus of this work, one needs to
solve the Helmholtz equation for multiple right-hand sides (sources) and
multiple frequencies at each iteration.

The efficient numerical solution of the Helmholtz equation with high and
spatially-dependent wavenumber is a difficult computational task. The linear
system obtained from the discretization of the Helmholtz equation involves a
very large, complex-valued, sparse, and indefinite matrix. The equation is
most difficult in three dimensions, where a direct solver is highly memory-
consuming and impractical in most scenarios. On the other hand, there are
multiple available iterative methods to solve the equation, but they typically
require many iterations at high wavenumbers. Either way, the data simulations
are highly cumbersome in FWI and are the main computational bottleneck of the
inversion.

Beyond the high computational costs of the forward simulations, FWI is a
highly non-linear and ill-posed problem, especially in the absence of low-
frequency data, which are hard to acquire because of physical limitations
[33]. Therefore, sophisticated regularization techniques [11], loss functions
[25], and careful minimization strategies [56, 7] are often considered. Still,
even after many years of research, the solution process of FWI is sensitive to
the inversion hyper-parameters, and cumbersome computationally.

In recent years, data-driven deep learning (DL) methods [28, 30, 54] have
revolutionized many scientific fields, most notably computer vision and
natural language processing. Inverse problems were also treated by DL, and FWI
in particular, because of their difficulty and non-linearity. In the context
of FWI, two main approaches are considered: (1) data-guided approaches (end-
to-end) where deep networks are trained to consume data on the one hand and
output a velocity model on the other [60, 18, 32], and (2) physics-guided
approaches where the data-fitting structure of the problem is kept and we make
sure that the resulting model fits the observed data in simulations [48, 13].
Here, we focus on the physics-guided approach, since the first approach is
highly dependent on the setup of the problem used during training, which may
be incompatible with real-life scenarios. This includes locations of sources
and receivers, types of media in the training set, missing data due to
technical problems, etc. In other words, end-to-end training may not be
flexible enough to handle different experiments. In the physics-guided
approach, we can be flexible with respect to the setup of the problem, and we
are faithful to the classical approach of fitting numerical simulations. The
latter is important for applications where the confidence and correctness of
the velocity estimation are more important than the computational cost of the
inversion, like in oil and gas exploration. In such applications, in
particular, the estimation is not required in near-real-time, like in other
vision applications that exploit the efficiency of the end-to-end approach.
Choosing the physics-guided approach, we integrate deep learning within the
optimization process of solving FWI.

In this work, we tackle the computational cost of the physical simulations in
FWI, by using the multigrid-augmented encoder-solver convolutional neural
network (CNN) approach to solve the Helmholtz equation, where two
convolutional neural networks are used to solve the equation [6, 31]. The
encoder network is used once-per-medium to encode context vectors, and the
solver network is used as a preconditioner to solve the system for the many
right-hand-sides, given the context vectors. The main observation in [6] is
that for a single given medium, the solver CNN can solve the linear system
efficiently to any accuracy, but it struggles when it is expected to
generalize over the velocity medium, which is generally unknown during
training. Hence, the advantage of the approach in our case stems from two
properties of the inverse solution:

  1. 1.

At each iteration of the inverse solution, we solve multiple systems with the
same medium but with different right-hand-sides. That is true even if we use
acceleration techniques such as simultaneous sources [23], and even more so if
we use high order methods like Gauss-Newton [42].

  2. 2.

During the inversion, the velocity model changes gradually, hence the
network’s weights do not need to change dramatically to cope with that and
maintain the solver’s efficiency. Hence, small training sessions between the
inversion iterations can dramatically speed up the solution process, relieving
the need for the CNN to generalize over velocity mediums.

To summarize, in this work we use the encoder-solver deep learning approach to
accelerate physics-guided FWI. We accelerate the simulations to ease the
overall inversion process significantly (e.g., inversion parameter tuning). We
do note that DL approaches can be highly beneficial in other aspects of FWI as
well, like regularization and data pre-processing [33], in addition to our
approach.

The rest of the work is organized as follows. In the following sections, we
review related work and discuss preliminaries and background that are needed
to introduce our approach. Next, we present our main contribution: a light-
weight data efficient retraining procedure of the encoder-solver model.
Lastly, we demonstrate the improvements and effectiveness of our approach when
solving the FWI problem. We include the results of experiments performed on
several 2D models with high-frequency data.

##  2 Related Work

#### End-to-end neural network FWI solvers.

With the advances of DL methods for the solution of inverse problems [36, 34],
and the concurrent rise in recorded seismic data, efficient DL models have
been proposed to solve the FWI problem. The DL approach involves training a
neural network and the solution is obtained by the network’s forward
application. That is, the DL model is viewed as an end-to-end solver and
entirely replaces the classical inversion process. Various models have been
proposed, including Multi-Layer Perceptron (MLP) [53, 26], Convolutional
Neural Network (CNN) [4, 40, 60, 24], Recurrent Neural Network (RNN) [45, 1]
and encoder-decoder models [20, 21, 18, 32]. However, end-to-end data driven
methods are highly dependent on the data and the training process of the
network. As a result, this approach suffers from a generalization gap when
faced with out-of-distribution data. Moreover, acquiring rich and expressive
training data is a difficult task that impacts the generalization capabilities
of the proposed models.

#### Physics-informed neural network.

To overcome the challenges of training the neural network posed by the data
acquisition process, [43] proposed a physics-informed neural network (PINN)
instead of using a pure data-mapping loss objective. In PINN, the underlying
physical laws is incorporated in the loss function to reconstruct NN-based
functional solutions to PDEs. That is, given spatial or temporal input data,
the network estimates the solution to the PDE at the coordinates of interest.
PINNs use automatic differentiation to represent all the differential
operators [9] and hence there is no explicit need for a mesh or grid data
representation [8, 55]. In geophysical applications, among others, the PINN
approach was shown to be effective for solving the forward problem, given by
the Eikonal equation [58] or by the Helmholtz equation [2, 49, 50, 51], and
for solving the inversion problem [53, 38, 44].

#### Classical inversion enhanced by deep learning.

A different approach aims to improve the performance of classical methods by
integrating DL networks in the inversion process. One of the proposed
improvements is to mitigate the issue of cycle-skipping in FWI [57]. [52, 39]
proposed a CNN-based model for low-frequency data extrapolation. [7] proposed
to perform wavelet source transformation [59] using a CNN model, which results
in a more robust inversion process. A different approach was proposed by [45],
where an RNN is utilized to model the forward wave equation and the
optimization is solved using stochastic approaches, such as ADAM [27] or SGD.
However, the inversion process is still computationally demanding and many
iterations are required to obtain a good velocity model estimation. Here, we
propose to improve the forward simulation, obtained by the solution of the
Helmholtz equation, using the encoder-solver architecture proposed by [6].
Furthermore, we introduce a data-efficient light-weight retraining process to
maintain the network performance throughout the inversion process and reduce
the overall computational complexity and solve-time of the inversion process.

##  3 Preliminaries and Background

###  3.1 The Full Waveform Inversion problem

Full waveform inversion (FWI) is a challenging inverse problem in which we
estimate the wave propagation velocity of a given medium using observed
wavefield data on the boundaries. We consider the frequency domain version of
FWI.

#### The forward problem: the Helmholtz equation.

To model the waveforms, we consider the discrete acoustic Helmholtz equation
as a forward problem, assuming a constant density media

$-\Delta u(x,\omega)-\omega^{2}m(x)(1-\gamma i)u(x,\omega)=g(x,\omega),\quad
x\in\Omega.$ (3.1)

The unknown $u(x,\omega)$ represents the pressure field in the frequency
domain, $\Delta$ is the Laplacian operator, $m(x)$ is the heterogeneous
squared slowness of the medium which is the physical parameter we wish to
estimate in this work. That is $m(x)=1/c(x)^{2}$ where $c(x)$ is the acoustic
wave velocity, which varies with position, and $\omega=2\pi f$ denotes angular
frequency, a scalar measure of rotation rate ($f$ is the frequency, measured
in Hertz in our context). The source term is denoted by $g$, and is assumed to
be a discretization of a delta function, $\delta(\hat{x}-x)$ located at $x$.
The parameter $\gamma$ indicates the fraction of attenuation (or damping) in
the medium, and can possibly be considered as varying attenuation, i.e.,
$\gamma(x)$. $i=\sqrt{-1}$ is the imaginary unit. Free space or absorbing
(non-reflecting) boundary conditions are applied at the boundaries
$\Gamma=\partial\Omega$.

Equation (3.1) is given in continuous space, and to solve it we discretize it
on a grid. We use the finite difference discretization, obtained by a standard
central difference scheme which leads to a system of linear equations

$\mathcal{H}({\bf m},\omega){\bf u}({\bf x},\omega)={\bf g}({\bf x},\omega).$
(3.2)

Here, boldface letters, e.g., ${\bf u}$, denote the functions sampled on the
grid. $\mathcal{H}$ is the linear operator defined by the finite difference
stencil of the Helmholtz equation. When discretizing the Helmholtz equation we
must use a small enough grid size such a wavelength is accomodated on several
grid nodes, depending on the discretization order (typically, 10 grid nodes
per wavelength are used). This typically requires a fine mesh for high
wavenumbers, making the grid and system $\mathcal{H}{\bf u}={\bf g}$ above
huge, involving possibly hundreds of millions of unknowns.

#### The inverse problem.

In FWI, sources are located at many locations on the top part of the domain,
and the waveform that is generated by each source is recorded at the locations
where receivers are placed. We denote a discrete solution of (3.2) as ${\bf
u}({\bf m},\omega,x_{s})$ for ${\bf g}$ that is a discretization of a source
at location $x_{s}$. In our formulation, each observed data sample
corresponding to a source $x_{s}$ and frequency $\omega_{j}$, is assumed to be
given by

${\bf d}^{\text{obs}}({\bf m}_{\text{true}},x_{s},\omega_{j})={\bf
P}_{s}^{\top}{\bf u}({\bf m}_{\text{true}},\omega_{j},x_{s})+\epsilon_{sj}$
(3.3)

where ${\bf P}_{s}^{\top}$ is a sampling matrix that measures the wave field
${\bf u}$ that is generated by the source at $x_{s}$ at the locations of the
receivers, and $\epsilon_{sj}$ is noise, which we assume to be i.i.d sampled
from a Gaussian distribution with zero mean and variance $\Sigma_{sj}$. Given
data for many sources and multiple frequencies, we aim to estimate the true
underlying velocity model, $\textbf{m}_{\text{true}}$, by minimizing the
difference between the measured and simulated data.

This may be done by solving the following PDE-constrained optimization problem

$\displaystyle\underset{{\bf m},\\{{\bf
u}_{sj}\\}_{s=1,j=1}^{n_{s},n_{f}}}{\arg\min}\Phi(\textbf{m},\\{\textbf{u}_{sj}\\})=\sum_{j=1}^{n_{f}}\sum_{s=1}^{n_{s}}\|{\bf
P}_{s}^{\top}{\bf u}_{sj}-{\bf
d}_{sj}^{\text{obs}}\|_{\Sigma_{sj}^{-1}}^{2}+\alpha\mathcal{R}({\bf m})$
(3.4) $\displaystyle\text{s.t. }\mathcal{H}({\bf m},\omega){\bf u}_{sj}={\bf
g}_{s}$

where ${\bf u}_{sj}$ is the simulated waveform for source $x_{s}$ and
frequency $\omega_{j}$, which is predicted for a given velocity model ${\bf
m}$, according to the forward problem defined in (3.1), the data terms
$\textbf{d}_{sj}^{\text{obs}}$ are the corresponding observed data as in
(3.3). Due to the fact that we assume ${\bf m}$ to be layered model we include
a regularization term, $\mathcal{R}({\bf m})$, which promotes piece-wise
smooth functions like total variation regularization [46]. Later, we detail
the regularization terms we use in our experiments).

There are several approaches for solving constrained optimization problem such
as (3.4). One of which is the Lagrange multipliers approach, where a Lagrange
variable will be introduced to each pair $s,j$ [22, 37]. Another common
approach, which we follow in this paper, is to eliminate the PDE constraint
and include it in the optimization formulation itself. Resulting in the
unconstrained formulation

$\underset{\textbf{m}}{\arg\min}\,\,\Phi({\bf
m})=\sum_{j=1}^{n_{f}}\sum_{s=1}^{n_{s}}\|{\bf P}_{s}^{\top}\mathcal{H}({\bf
m},\omega_{j})^{-1}{\bf g}_{s}-{\bf
d}_{sj}^{\text{obs}}\|_{\Sigma_{sj}^{-1}}^{2}+\alpha\mathcal{R}({\bf m}).$
(3.5)

###  3.2 The Simultaneous Sources Approach

A computationally efficient approach for solving (3.5) utilizes the concept of
simultaneous sources, which is obtained by trace estimation [23]. That is,
given a matrix ${\bf A}\in\mathbb{R}^{m\times n}$ we can approximate its
Forbenius norm by [5]

$\|{\bf A}\|_{F}^{2}=trace({\bf A}^{T}{\bf A})=\mathbb{E}_{x}\|{\bf A}{\bf
x}\|_{2}^{2}\approx\frac{1}{p}\sum_{i=1}^{p}\|{\bf A}{\bf
x}_{i}\|_{2}^{2}=\frac{1}{p}\|{\bf A}{\bf X}\|_{F}^{2}$ (3.6)

where each ${\bf x}_{i}\in\mathbb{R}^{n}$ is chosen from Rademacher
distribution and ${\bf X}$ is $n\times p$ matrix whose $i$-th column is ${\bf
x}_{i}$. This approximation allows us to approximate the norm of a large
$m\times n$ matrix using an $m\times p$ matrix, which can significantly reduce
computational cost when $p\ll n$.

Utilizing (3.6) for approximating the data term in (3.5) is only possible when
$\textbf{P}_{s}^{\top}$ and $\Sigma_{sj}$ are not source dependent, i.e., when
all the receivers record the waveform from all the sources. To mitigate this
issue, [35] introduces a new set of variables $\hat{\textbf{d}}_{sj}$ that are
defined over all the locations of the receivers for all the sources, and an
additional penalty term is added to ensure that $\hat{\textbf{d}}_{sj}$ is
similar to the original $\textbf{d}_{sj}^{\text{obs}}$ at the locations of the
receivers that actually record the waveform for each source $s$. Thus, (3.5)
can be formulated as

$\underset{{\bf m},\\{\hat{{\bf
d}}_{sj}\\}_{s=1,j=1}^{n_{s},n_{f}}}{\arg\min}\Phi({\bf
m})=\sum_{j=1}^{n_{f}}\sum_{s=1}^{n_{s}}\|{\bf P}^{\top}\mathcal{H}({\bf
m},\omega_{j})^{-1}{\bf g}_{s}-\hat{{\bf
d}}_{sj}\|_{\hat{\Sigma}_{j}^{-1}}^{2}+\eta\|\hat{{\bf P}}_{s}^{\top}\hat{{\bf
d}}_{sj}-{\bf
d}_{sj}^{\text{obs}}\|_{\Sigma_{sj}^{-1}}^{2}+\alpha\mathcal{R}({\bf m})$
(3.7)

where ${\bf P}$ is an operator that projects a vector onto the union of the
receivers’ locations, and $\hat{{\bf P}}_{s}$ chooses the subset of receivers
for source $s$ out of that union, i.e., ${\bf P}_{s}^{\top}=\hat{{\bf
P}}_{s}{\bf P}^{\top}$. To apply the simultaneous sources technique to (3.7)
we define a $n_{s}\times p$ matrix ${\bf X}$ from Radermacher distribution,
and by (3.6) we obtain

$\displaystyle\sum_{s=1}^{n_{s}}\|{\bf P}^{\top}\mathcal{H}({\bf
m},\omega_{j})^{-1}{\bf g}_{s}-\hat{{\bf d}}_{sj}\|_{\Sigma_{sj}^{-1}}^{2}$
$\displaystyle=\|{\bf P}^{\top}\mathcal{H}({\bf m},\omega_{j})^{-1}{\bf
G}-\hat{{\bf D}}_{j}\|_{\Sigma_{j}^{-1}}^{2}$ (3.8)
$\displaystyle\approx\frac{1}{p}\|{\bf P}^{\top}\mathcal{H}({\bf
m},\omega_{j})^{-1}{\bf G}{\bf X}-\hat{{\bf D}}_{j}{\bf
X}\|_{\Sigma_{j}^{-1}}^{2}$

This effectively reduces the number of sources to $p$ in every iteration.

###  3.3 Frequency continuation and Gauss-Newton

Since FWI is ill-posed and highly non-convex, local minima can pose major
obstacles when solving FWI by optimization. That is, converging to a local
minimum may lead to an non-plausible estimated slowness model. Therefore,
avoiding local minima can be done effectively by using the method of frequency
continuation (FC) [42, 41]. As shown in algorithm 1, the FC algorithm iterates
over ranges of frequencies, called cycles, each defined by an initial and
final frequency. In each such cycle, we solve the optimization problem for the
given frequencies in a sliding window manner using the model estimated by the
previous cycle as an initial guess for the current cycle. Solving the
optimization problem with respect to $\mathbf{m}$ is done iteratively using
Gauss-Newton (GN) [42]. Specifically, at each iteration $t$ of GN we obtain
the linear approximation

${\bf u}_{sj}({\bf m}^{(t)}+\delta{\bf m})={\bf u}_{sj}({\bf m}^{(t)})+{\bf
J}_{sj}({\bf m}^{(t)})\delta{\bf m}$ (3.9)

where ${\bf u}_{sj}({\bf m}^{(t)})=\mathcal{H}({\bf
m}^{(t)},\omega_{j})^{-1}({\bf g}_{s})$ and the Jacobian matrix ${\bf
J}_{sj}({\bf m}^{(t)})$ is defined as

${\bf J}_{sj}({\bf m}^{(t)})=\nabla_{{\bf m}}{\bf
u}_{sj}=-\omega_{j}^{2}\,\mathcal{H}({\bf
m},\omega_{j})^{-1}\text{diag}(\mathcal{H}({\bf m},\omega_{j})^{-1}{\bf
g}_{s})$ (3.10)

Input:

$\omega_{1}<\ldots<\omega_{n_{f}}$: frequencies

$ws$: window size

$\mathcal{C}=[(i_{start}^{(1)},i_{end}^{(1)}),\dots,(i_{start}^{(k)},i_{end}^{(k)})]$:
set of cycles

Algorithm:

for  _$(i_{start},i_{end})\in\mathcal{C}$_ do

Initialize m (or use the model from previous cycle)

for  _$i=i_{start}:i_{end}$_ do

$\textbf{m}\leftarrow$ Approximately solve $\Phi$ using data for the
frequencies.

$\omega_{\max\\{i-ws,1\\}}\ldots\omega_{i}$, starting from the previous model
m.

end for

end for

Algorithm 1 Frequency Continuation

At each iteration we solve the forward problem (3.1) for each pair of
frequency and source, whether encoded through the simultaneous sources method
or not. At each GN iteration one has to solve the linear system ${\bf
H}\delta{\bf m}=-\nabla_{{\bf m}}\Phi({\bf m}^{(t)})$, where ${\bf H}$ is the
Gauss-Newton Hessian and is defined by

${\bf H}=\sum_{s,j}{\bf J}_{s,j}({\bf m}^{(t)})^{\top}{\bf P}{\bf
P}^{\top}{\bf J}_{s,j}({\bf m}^{(t)})+\alpha\nabla^{2}\mathcal{R}({\bf
m}^{(t)}).$ (3.11)

The system here is solved by the Conjugate Gradient algorithms, involving many
applications of the Hessian to vectors, which requires solving the forward and
the adjoint systems involving $\mathcal{H}$ and $\mathcal{H}^{*}$,
respectively. We note that the right-hand-side vectors for the Jacobians are
arbitrary (that is, not the source functions) as they come from the
optimization process.

###  3.4 Shifted Laplace multigrid-augmented CNN preconditioner

The Shifted Laplace operator in [17] was proposed to allow the solution of the
Helmholtz equation using multigrid, in a Krylov subspace method (e.g. GMRES).
The preconditioner is given by:

$P{\bf u}=-\Delta_{h}{\bf u}-\omega^{2}{\bf m}({\bf x})(\alpha-\beta i){\bf
u},\quad\alpha,\beta\in\mathbb{R}$ (3.12)

where typically, $\alpha=1$ to maintain similarity to $\mathcal{H}$, and
$\beta=0.5$ which was shown in [17] to lead to an efficient and robust
preconditioning operator. However, the SL method (similar to relaxation
methods) efficiently handles error modes with rather large absolute
eigenvalues, as the shift does not perturb those eigenvalues by a significant
factor. Therefore, [6] introduced a CNN-based augmentation to the SL
preconditioner, such that the augmentation focuses on ”algebraically smooth”
errors. The augmentation is obtained by an Encoder-Solver model, depicted in
Figure 1, responsible for mapping residual vectors generated by a Krylov
method to the corresponding error vectors, which is then followed by a
multigrid cycle [12].

Figure 1: The Encoder-Solver architecture. The Encoder receives the slowness
(squared) and progressively coarsens the feature map using a convolutional
layer followed by 2 consecutive ResNet layers. On the way up we up-sample the
coarse vector using a convolutional layer, concatenate the feature vector of
the corresponding size, and apply a ResNet layer. The Solver network maps a
residual vector ${\bf r}$ to an error ${\bf e}$. The Solver integrates the
learned feature vector of the Encoder at each level.

The Encoder-Solver model consists of two separate U-Net models. The Encoder
U-Net receives only the squared slowness model ${\bf m}$ as input and
generates context vectors

$\\{{\bf z}_{i}\\}=\text{EncoderNet}({\bf m};\boldsymbol{\theta}_{e}),$ (3.13)

aiming to generalize solely over the slowness model. The Solver U-Net
integrates the learned feature maps $\\{{\bf z}_{i}\\}$ from the encoder when
generalizing over the residual vectors. This way the solver network has more
information on the slowness model, and it allows us to compute the feature
maps $\\{{\bf z}_{i}\\}$ of a given slowness model only once, and solve for
any right-hand-side residual vector. This property is highly desirable in the
case of FWI where the Helmholtz equation is solved for numerous right-hand-
side vectors for a given slowness model. Overall, the preconditioning process,
denoted as $M_{VU}$, is given by

$M_{VU}(\textbf{r})=\text{V-cycle}(\textbf{e}_{0}=\text{SolverNet}(\textbf{r},\gamma;\\{{\bf
z}_{i}\\},\boldsymbol{\theta}_{s}),\textbf{r},\text{levels}=3)$ (3.14)

acting on a residual vector ${\bf r}$ generated by a Krylov method.

##  4 Method

Our method incorporates the U-Net Encoder-Solver preconditioner of the
Helmholtz equation in the inversion process, both in the data simulations and
in the sensitivity (Jacobian) computations in equation (3.10). The CNN network
in [6] was proposed in order to solve general Helmholtz problems, and hence,
was trained over a general data set of slowness models for a particular
setting of the problem (e.g., grid-size, domain, and proportion between the
horizontal and vertical grid sizes $h_{x}$ and $h_{y}$). Any deviation from
these settings can hinder the effectiveness of the approach. Furthermore, it
is shown that the performance of the model deteriorates when solving an
equation corresponding to an out-of-distribution slowness model. Namely, the
generalization over the slowness model hurts the performance. It is also shown
that briefly re-training the network on the slowness model of interest
reconstitutes accuracy for this model.

Here, we wish to exploit the fact that we solve an inverse problem to estimate
a single (yet, unknown) slowness model ${\bf m}({\bf x})$ for a given problem
setting (e.g., domain, grid-size), and that we know the initial guess of the
inversion process. We build on the fact that the model is built gradually and
does not change dramatically between the frequency continuation iterations.
Hence, to cope with the gradual change in the underlying slowness model we
introduce and integrate a light-weight retraining procedure of the Encoder-
Solver model in between the inversion iterations so that its weights are
updated together with the inversion to solve the forward problems well.

In summary, we employ the Encoder-Solver model as a preconditioner to the
Helmholtz equation, where the model has been trained to consume general right-
hand-sides (residuals) generated by a Krylov method until convergence. The
advantages of this approach in our context are as follows:

  1. 1.

Since there are many sources in typical FWI scenarios, we wish to train one
CNN model for all sources. Also, the CNN model should allow us to solve the
discrete Helmholtz equation to any desired accuracy as the estimation of deep
reflectors requires quite accurate simulations. Accurate solutions are also
required for a consistent gradient computation in (3.10). Applying the CNN as
a preconditioner allows us to reach any desired accuracy.

  2. 2.

It is highly preferable that the model can solve problems with general right-
hand sides. This is required for the multiplication of the Jacobian in (3.10)
and its transpose with a vector when computing the gradient or when solving
the inner Newton problem with CG. Furthermore, the forward solutions for
random sources are also needed for the data simulation in the simultaneous
sources method detailed in section 3.2, where in each iteration we encode the
sources using a random matrix.

We further elaborate on the training process of the Encoder-Solver model and
the data generation scheme in Section 4.1. Next, in Section 4.2 we present our
retraining integrated frequency continuation procedure. Lastly, in Section 4.3
we show how the Encoder-Solver model is able to solve the adjoint system
efficiently.

###  4.1 Training and data generation

#### Training.

The training of the Encoder-Solver model is performed in a supervised manner.
Specifically, let

${\bf e}^{\text{net}}({\bf r},{\bf
m},\gamma;\boldsymbol{\theta})=\text{SolverNet}({\bf
r},\gamma;\text{EncoderNet}({\bf
m};\boldsymbol{\theta}_{e}),\boldsymbol{\theta}_{s})$ (4.1)

be a forward application of the network for a given residual vector ${\bf r}$,
slowness (squared) model ${\bf m}({\bf x})$ and attenuation model $\gamma({\bf
x})$ (which is the same for every sample in the dataset).
$\boldsymbol{\theta}$ denotes the set of trainable weights of both networks
$\\{\boldsymbol{\theta}_{e},\boldsymbol{\theta}_{s}\\}$. We seek to minimize
the mean squared error (MSE)

$\underset{\boldsymbol{\boldsymbol{\theta}}}{\min}\,\frac{1}{m}\sum_{i=1}^{m}\|{\bf
e}^{\text{net}}({\bf r}_{i},{\bf m}_{i},\gamma;\boldsymbol{\theta})-{\bf
e}_{i}^{\text{true}}\|_{2}^{2}$ (4.2)

for each batch of data samples $\\{({\bf e}_{i}^{\text{true}},{\bf r}_{i},{\bf
m}_{i},\gamma)\\}_{i=1}^{m}$, where $\mathcal{H}_{i}^{h}{\bf
e}_{i}^{\text{true}}={\bf r}_{i}$. $\mathcal{H}_{i}^{h}$ is the Helmholtz
equation descretized operator as defined in 3.2 for the given slowness model
${\bf m}_{i}(\textbf{x})$. The minimization is performed by a stochastic
gradient descent optimizer, wherein we sweep through all the training examples
in batches each epoch. We use the ADAM optimizer [27] with learning rate
decay.

#### Data Generation.

Creating the dataset for training our network requires creating data samples
of the form $({\bf e}_{i}^{\text{true}},{\bf r}_{i},{\bf m}_{i},\gamma)$.
Where $({\bf r}_{i},{\bf m}_{i},\gamma)$ is the input to the network and ${\bf
e}_{i}^{\text{true}}$ is the corresponding error vector used to compute the
misfit term in (4.2). $\gamma({\bf x})$ is an absorbing layer (ABL) [15, 16,
17] that we use to model open domains in (3.1). Specifically, in our case,
$\gamma({\bf x})$ is a two dimensional function that goes from $0$ to $1$
towards the boundaries, except for the top boundary, which indicates the sea
layer, which we allow reflections from. $\\{{\bf m}_{i}({\bf x})\\}$ are
spatially dependent random linear slowness (squared) models (e.g., see Fig.
2(b)) and are generated with correspondence to the initial guess used for the
FWI optimization process. That is, we randomly select physically applicable
values for the top and bottom layers, and the slowness model is obtained by
uniform intervals between these two values. Lastly, since the relation between
the error vector and the residual vector is described as a linear system,
$\mathcal{H}_{i}^{h}{\bf e}_{i}^{\text{true}}={\bf r}_{i}$, the creation of
${\bf r}_{i}$ and ${\bf e}_{i}^{\text{true}}$ pairs is straightforward, and
does not involve solving the PDE directly to obtain ground-truth solution.
However, because we wish to focus on algebraically smooth errors (errors with
low residuals), we adapt our data by applying a random number of GMRES
iterations with a V-cycle preconditioner. That is, we first generate a random
vector ${\bf x}_{i}$ from a complex normal distribution, which is in
correspondence to the random right-hand-sides vectors we solve for throughout
the optimization process. Then we compute the right-hand-side vector ${\bf
b}_{i}=\mathcal{H}_{i}^{h}{\bf x}_{i}$ and apply

$\tilde{{\bf x}}_{i}=\text{GMRES}(\mathcal{H}_{i}^{h},M=\text{V-cycle},{\bf
b}_{i},{\bf x}^{(0)}=\textbf{0},iter\in\\{1,\dots,10\\}).$ (4.3)

Following that, we compute ${\bf r}_{i}={\bf
b}_{i}-\mathcal{H}_{i}^{h}\tilde{{\bf x}}_{i}$ and ${\bf
e}_{i}^{\text{true}}={\bf x}_{i}-\tilde{{\bf x}}_{i}$.

###  4.2 Retraining integrated frequency continuation

Here we describe the key idea for utilizing the Encoder-Solver described above
for solving FWI most effectively. The problem with the CNN-based solver above
is that it cannot be highly effective on all slowness models, but given one
slowness model and proper training, it performs well. We wish to exploit this
property through mini-retraining sessions of the solver against the gradually
evolving inversion iterates as detailed below.

In our method, we solve FWI using the frequency continuation method in
Algorithm 1. As the algorithm progresses the slowness model is gradually
changed to fit the observed data. Thus, the slowness model gradually deviates
from the initial linear model distribution that our solver was trained on. As
a result, the Encoder-Solver model is required to generalize over unseen
slowness models, which leads to poor performance. However, when a light-weight
retraining procedure is applied, the performance of the model is kept and the
iteration count needed for the solution of the Helmholtz equation remains low.
To this end, we incorporate a retraining procedure when a new slowness model
is obtained, i.e., after each inner iteration within the frequency
continuation algorithm.

At this point, clearly, there is a trade-off between the cost of the
retraining procedure and the improvement of our solver network. The more
resources we invest in the training procedure the more our solver is able to
adapt to the new slowness model. However, our main goal is to control the
computational resources used throughout the inversion process, which is one of
the main obstacles of FWI. In addition to the cost of the training itself a
major part of the computational cost is due to the process of data creation,
which is described in Subsection 4.1. Therefore, our proposed retraining
procedure is data efficient and computationally light, and is described in
Algorithm 2.

Input:

epochs: total number of epochs.

data_epochs: number of initial epochs that also generate data.

initial_size: the size of the initial data set.

Algorithm:

Create an initial dataset $\mathcal{D}$ corresponding to ${\bf m}$ of size
initial_size

for  _$i=1:epochs$_ do

Train model on $\mathcal{D}$

if  _$i\leq data\\_epochs$_ then

foreach  _$({\bf r},{\bf e},\gamma,{\bf m})\in\mathcal{D}$_ do

${\bf e}_{\text{new}}\leftarrow\text{Encoder-Solver}({\bf r},{\bf m},\gamma)$

$\tilde{{\bf e}}\leftarrow\text{GMRES}(\mathcal{H}^{h},M=\text{V-cycle},{\bf
r},{\bf x}^{(0)}={\bf e}_{\text{new}},iter=1)$

$\tilde{{\bf r}}\leftarrow\mathcal{H}^{h}\tilde{{\bf e}}$

Append $(\tilde{{\bf r}},\tilde{{\bf e}},\gamma,{\bf m})$ to $\mathcal{D}$

end foreach

end if

end for

Algorithm 2 Light-weight Retraining Procedure

To start the algorithm we create a small dataset, and the rest of the data
samples are gradually created throughout the retraining process in a manner
that is compatible with our solver inaccuracies. That is, for several initial
epochs we re-use the network outputs obtained in the training as new initial
guesses in the process of creating a new $({\bf r},{\bf e})$ corrected pair.
The process resembles ideas presented in [19], where the network is used to
generate more data to obtain a fixed-point process of implicit (or, recurrent)
neural networks. Specifically, in our method we initially create dataset of
size $initial\\_size=128$ and gradually create the rest of the data samples
for the first $data\\_epochs=3$ epochs, resulting in a final dataset of $1024$
samples. Then, we continue to train our model for the rest of the epochs until
reaching our epoch limit of 30. Note that in the data creation epochs we
obtain a new error vector using GMRES, but with one iteration only, which help
reduce the computational cost and still improve our Encoder-Solver
preconditioner performance. Moreover, the final dataset size (1024) is
considered small compared to the initial dataset size that was used to train
the network, keeping the retraining process at each iteration computationally
light-weight.

###  4.3 The adjoint system

Throughout the inversion process, one also has to solve the adjoint system
involving $\mathcal{H}^{\ast}$, as shown in (3.11) and (3.10). However, the
Encoder-Solver networks are trained on data created by the system obtained by
$\mathcal{H}$, and we do not wish to train two sets of networks to support
both the standard and adjoint Helmholtz systems. Thus, in order to adapt our
solver for the adjoint system we alter the right-hand-side vector of the
system rather than training the Encoder-Solver on more data or have two sets
of networks. That is, suppose we have the adjoint system

$\mathcal{H}^{\ast}{\bf x}={\bf b}$ (4.4)

we can re-write it with the real and imaginary terms, resulting in the
equation

$\left(\text{re}(\mathcal{H}^{\top})-\text{im}(\mathcal{H}^{\top})i\right)(\text{re}({\bf
x})+\text{im}({\bf x})i)=\text{re}({\bf b})+\text{im}({\bf b})i$ (4.5)

where $\text{re}(\cdot),\text{im}(\cdot)$ represent the real and imaginary
parts, respectively, and $i$ is the imaginary unit. The Helmholtz operator
matrix $\mathcal{H}$ has symmetric real and imaginary parts, i.e.,
$\text{re}(\mathcal{H}^{\top})=\text{re}(\mathcal{H})$ and
$\text{im}(\mathcal{H}^{\top})=\text{im}(\mathcal{H})$. Thus, equation (4.5)
translates to

$\left(\text{re}(\mathcal{H})-\text{im}(\mathcal{H})i\right)(\text{re}({\bf
x})+\text{im}({\bf x})i)=\text{re}({\bf b})+\text{im}({\bf b})i$ (4.6)

and we obtain a system of equations, for the real and imaginary part,
separately

$\displaystyle(1)\quad\text{re}(\mathcal{H})\text{re}({\bf
x})+\text{im}(\mathcal{H})\text{im}({\bf x})=\text{re}({\bf b})$ (4.7)
$\displaystyle(2)\quad\text{re}(\mathcal{H})\text{im}({\bf
x})-\text{im}(\mathcal{H})\text{re}({\bf x})=\text{im}({\bf b}).$

Now, performing $(1)-i(2)$ yields the equation

$\mathcal{H}{\bf x}^{\ast}={\bf b}^{\ast}.$ (4.8)

which includes the standard $\mathcal{H}$ with conjugate solution and right-
hand-side. Therefore, the result for the original adjoint system presented in
(4.4) is the conjugate of the vector obtained by solving (4.8), i.e., $({\bf
x}^{\ast})^{\ast}={\bf x}$, as required. Since our Encoder-Solver model was
trained on random right-hand-side vectors obtained by the operator
$\mathcal{H}$, equation (4.8) is compatible with our solver.

##  5 Numerical Results

In this section, we report our results and demonstrate the efficiency of the
proposed technique, focusing mainly on the improved performance of our forward
solver network throughout the inversion process. In our experiments, we
perform FWI with simultaneous sources only, and we do not augment the
inversions with other modalities or techniques known in the literature. In
principle, such techniques can be used in addition to our approach in more
complicated real-life experiments. We conduct our experiments on the 2D
SEG/EAGE salt model [3] and on the Marmousi model [10]. We compare the
iteration count and solve time for solving the FWI problem with and without
our proposed retraining procedure. All of our experiments and time
measurements were conducted using NVIDIA RTX-3090 GPU operated by the Julia
Flux framework.

Initial training details: We train our network on a dataset consisting of
$20,000$ samples, and a validation set of $1,000$ samples, using grid size of
$608\times 304$. We train for $120$ epochs, the batch size is set to $16$ and
the initial learning rate is set to $10^{-4}$.

Retraining details: We perform the retraining procedure after each window in
the frequency continuation algorithm 1 (after each inner iteration).
Therefore, we train the Encoder-Solver model using data samples with grid size
that corresponds to the highest frequency of the window. Namely, grid size
that obtains 10 grid points per wavelength. This further reduce the
computational cost of the retraining procedure while assuring the accuracy of
the Encoder-Solver model.

###  5.1 Regularization terms used for the reconstruction

The objective functions in the FWI formulations in equations (3.4), (3.5), and
(3.7) contain a regularization term $\mathcal{R}({\bf m})$. Based on the work
[56] we apply two regularization functions in our experiments. The first of
which is a high order regularization called spline smoothing, given by

$\mathcal{R}_{high}({\bf m})=\|\Delta_{h}({\bf m}-{\bf m}_{ref})\|_{2}^{2}$
(5.1)

The goal of this regularization is to create a smooth model from high-
frequency data. The reference model ${\bf m}_{ref}$ is set to be the initial
guess for the inversion, and while using this regularization, we keep ${\bf
m}_{ref}$ fixed. We use (5.1) to obtain a good smooth model, so that the rest
of the process will result in a plausible reconstruction and avoid local
minima. The second regularization function is a standard diffusion
regularization

$\mathcal{R}_{diff}({\bf m})=\|\nabla_{h}({\bf m}-{\bf m}_{ref})\|_{2}^{2}$
(5.2)

where $\nabla_{h}$ represents a discrete gradient on a nodal grid. Using (5.2)
results in a rather sharp reconstruction and is suitable after a smooth model
is constructed.

###  5.2 SEG/EAGE salt model setup and results

Our first set of experiments is conducted using the SEG/EAGE model presented
in Fig. 2(a). The model is described by a $608\times 304$ grid, representing
an area of size $13.5km\times 4.2km$. The data is generated by first solving
(3.1) for frequencies $\omega_{j}=2\pi f_{j}$, where

$\\{f_{j}\\}=\\{\underbrace{2.6,2.9,3.2,3.5}_{low},\underbrace{3.9,4.2,4.5,4.8,5.2,5.5}_{middle},\underbrace{5.8,6.1,6.4,6.7}_{high}\\}$
(5.3)

We then add $1\%$ Gaussian noise to the data. The initial model for all the
inversions is given in Fig. 2(b), and is initially used as ${\bf m}_{ref}$ in
the smoothing regularization. We apply five frequency continuation sweeps \-
the first three using the low frequencies with the regularization (5.1),
followed by two sweeps using low and middle frequencies with the
regularization (5.2) to obtain a sharp model. For the first three sweeps we
used 10 GN iterations per outer iteration, and for the last two we used 15 GN
iterations. To finalize the inversion, we applied up to 15 additional GN
iterations involving the four highest frequencies, in order to further sharpen
the reconstructed model. For all the sweeps we used 5 Conjugate Gradient (CG)
iterations in each formulations.

(a) True velocity model

(b) Reference initial model

Figure 2: The 2D SEG/EAGE salt velocity model and initial guess

###  5.3 FWI with simultaneous sources for the SEG/EAGE model

For the SEG/EAGE salt model experiment, we use simultaneous sources. We set
the size of X in (3.6) to be $n_{s}\times p$, with $p=16$. We see that the
reconstructed model, shown in Fig. 3(a), estimates the ground truth salt
model, 2(a), with high accuracy. This similarity is also evident in the total
misfit history in Fig. 3(b).

(a) FWI estimated SEG/EAGE model

(b) Misfit history

Figure 3: SEG/EAGE FWI result and the corresponding total misfit history

###  5.4 SEG/EAGE computational cost and solve time

To compare the computational cost and the solve time of the FWI problem with
and without our proposed retraining procedure we compare the iteration count
and solve time of each optimization problem in the frequency continuation
algorithm 1. As described, we perform multiple frequency continuation cycles
where in each cycle, we iterate over the available frequencies in a sliding
window fashion and use those frequencies to solve an optimization problem. In
Fig. 4, we compare the average iteration count and average solve-time of each
window in the inversion process. We plot the average cost, with the
corresponding standard deviation of each window, for the different tolerance
thresholds, $10^{-6}$ and $10^{-4}$, used throughout the solution process for
the forward model computation and sensitivities, respectively. The drop
depicted in the graphs is due to the transition between the last window in the
fourth cycle, containing the highest middle part frequencies, and the first
window of the fifth cycle, containing the four lowest frequencies (low part).
We can clearly see that by applying our light-weight retraining procedure, we
obtain a computational and time efficient FWI solution process. Moreover, we
notice that by applying the retraining procedure, our network has less
variance, which indicates a unifying process of adapting the network to the
different grid resolutions.

In addition, we compare the total iteration count and solve-time needed for
the entire FWI optimization process in Fig. 5. The graphs plot the accumulated
iteration count and solve-time of the FWI process, including the time of the
retraining procedure. By applying the retraining procedure we obtained the
solution of the FWI process in $3.73\cdot 10^{5}$ iterations compared to
$1.13\cdot 10^{6}$ without retraining, and, in a total solve-time of
$6.708\cdot 10^{4}$ seconds compared to $1.846\cdot 10^{5}$. That is, we
obtained an improvement of $7.57\cdot 10^{5}$ iterations and $1.2419\cdot
10^{5}$ seconds, which amounts to a $67\%$ and $63\%$ speed-up for iteration
count and solve-time, respectively.

Figure 4: SEG/EAGE salt model FWI average comparison. In the left column, we
compare the average iteration count with and without training our model. The
right column shows a comparison of the average solve time. Both comparisons
were conducted for the different tolerance values. Figure 5: SEG/EAGE salt
model FWI total iteration count and solve-time comparison.

###  5.5 Marmousi model setup and results

Our next set of experiments is conducted using the Marmousi model presented in
Fig. 6(a). The model is described by a $544\times 272$ grid, representing an
area of size $9.334km\times 2.904km$. The data is generated by first solving
(3.1) for frequencies $\omega_{j}=2\pi f_{j}$, where

$\\{f_{j}\\}=\\{\underbrace{2.9,3.3,3.7,4.2}_{low},\underbrace{4.6,5.0,5.4,5.8,6.2,6.7,7.1}_{middle},\underbrace{7.5,7.9,8.3,8.7}_{high}\\}.$
(5.4)

We then add $1\%$ Gaussian noise to the data. The initial model for all the
inversions is given in Fig. 6(b), and is initially used as ${\bf m}_{ref}$ in
the smoothing regularization. We apply three frequency continuation sweeps \-
the first sweep using the first four frequencies, followed by two sweeps using
all frequencies. To finalize the inversion, we applied up to 30 additional GN
iterations involving the four highest frequencies, in order to further sharpen
the reconstructed model. For all the sweeps we used 10 GN iterations and 5
Conjugate Gradient (CG) iterations with the diffusion regularization term
(5.2).

(a) True velocity model

(b) Reference initial model

Figure 6: The 2D Marmousi velocity model and initial guess

###  5.6 FWI with simultaneous sources for the Marmousi model

For the Marmousi model experiment, we use simultaneous sources and we follow
the same parameters as in the SEG/EAGE model experiments. Specifically, we set
$p=16$ for the simultaneous sources formulation. As before, we see that the
reconstructed model, shown in Fig. 7(a), estimates the ground truth model,
6(a), with high accuracy. This similarity is also evident in the misfit
history in Fig. 7(b).

(a) FWI estimated Marmousi model

(b) Misfit history

Figure 7: Marmousi FWI result and the corresponding misfit history

###  5.7 Marmousi Computational Cost and Solve Time

Similarly to previous experiments, we show the comparison of solving the FWI
problem with and without our retraining procedure for the Marmousi model. The
Marmousi model contains high contrast areas and is highly heterogeneous. As
depicted in Fig. 8 and Fig. 9, without adjusting the weights of the Encoder-
Solver the solver reaches stagnation at an early stage of the inversion
process. In our experiments we have limited the iterative solver of the
Helmholtz equation to $300$ iterations, which is computationally expensive and
is considered unfeasible for large scale data. However, posing a lower
iteration limit to the solver will lead to a non-plausible reconstruction due
to the Encoder-Solver failing to converge to the required threshold accuracy.
On the other hand, when incorporating our light-weight retraining procedure,
not only the Encoder-Solver does not stagnate, it also obtains low iteration
counts and significantly low solve time. We see that even for highly
heterogeneous medium the solver is able to cope with the gradual changes and
obtain high accuracy with low iteration count. Moreover, similarly to the
SEG/EAGE experiments, we notice that a unifying process of the network’s
performance is obtained here as well. Namely, the variance of the iteration
count and solve-time remains significantly low, indicating the ability of the
retraining procedure to adapt the network to different grid resolutions for
the case of highly heterogeneous slowness models such as the Marmousi model.

Figure 8: Marmousi model FWI test. In the left column we compare the average
iteration count with and without training our model. The right column shows a
comparison of the average solve time. Both comparisons were conducted for the
different tolerance values. The experiments without retraining were stopped
due to the effective stagnation of the solver. Figure 9: Marmousi model FWI
total iteration count and solve-time comparison. The experiments without
retraining were stopped due to the effective stagnation of the solver.

##  6 Conclusion

In this work, we considered the integration of DL methods in the inversion
process of FWI under the physics-guided approach. Namely, we integrated the
Encoder-Solver model, presented in [6], in the data simulation throughout the
inversion process. Our solver is applied to both the forward problem and its
adjoint through manipulating the right-hand-sides. We also presented a
lightweight data-efficient retraining procedure to the Encoder-Solver model.
The retraining procedure allows the efficient integration of the Encoder-
Solver model in the inversion process of the FWI problem. Our proposed
retraining procedure demonstrated significant improvements to the Encoder-
Solver performance, especially when solving for highly contrasted mediums such
as the Marmousi model. Together with source-encoding, this yields the solution
of physics-guided FWI at a more reasonable cost than traditional methods.

## References

  * [1] Amir Adler, Mauricio Araya-Polo, and Tomaso Poggio.  Deep recurrent architectures for seismic tomography.  In 81st EAGE Conference and Exhibition 2019, pages 1–5. European Association of Geoscientists & Engineers, 2019.
  * [2] Tariq Alkhalifah, Chao Song, U Bin Waheed, and Qi Hao.  Wavefield solutions from machine learned functions that approximately satisfy the wave equation.  In EAGE 2020 annual conference & exhibition online, pages 1–5. European Association of Geoscientists & Engineers, 2020.
  * [3] F Aminzadeh, J Brac, and T Kunz.  3-d salt and overthrust models: Soc.  Expl. Geophy, 1997.
  * [4] Mauricio Araya-Polo, Joseph Jennings, Amir Adler, and Taylor Dahlke.  Deep-learning tomography.  The Leading Edge, 37(1):58–66, 2018.
  * [5] Haim Avron and Sivan Toledo.  Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix.  Journal of the ACM (JACM), 58(2):1–34, 2011.
  * [6] Yael Azulay and Eran Treister.  Multigrid-augmented deep learning preconditioners for the Helmholtz equation.  SIAM Journal on Scientific Computing, 0:S127–S151, 2022.
  * [7] Chenglong Bao, Lingyun Qiu, and Rongqian Wang.  Robust full waveform inversion: A source wavelet manipulation perspective.  SIAM Journal on Scientific Computing, 45(6):B753–B775, 2023.
  * [8] Leah Bar and Nir Sochen.  Strong solutions for pde-based tomography by unsupervised learning.  SIAM Journal on Imaging Sciences, 14(1):128–155, 2021.
  * [9] Atilim Gunes Baydin, Barak A Pearlmutter, Alexey Andreyevich Radul, and Jeffrey Mark Siskind.  Automatic differentiation in machine learning: a survey.  Journal of machine learning research, 18(153):1–43, 2018.
  * [10] A Brougois, Marielle Bourget, Patriek Lailly, Michel Poulet, Patrice Ricarte, and Roelof Versteeg.  Marmousi, model and data.  In EAEG workshop-practical aspects of seismic data inversion, pages cp–108. European Association of Geoscientists & Engineers, 1990.
  * [11] Sagi Buchatsky and Eran Treister.  Full waveform inversion using extended and simultaneous sources.  SIAM Journal on Scientific Computing, 2021.
  * [12] Henri Calandra, Serge Gratton, Xavier Pinel, and Xavier Vasseur.  An improved two-grid preconditioner for the solution of three-dimensional helmholtz problems in heterogeneous media.  Numerical Linear Algebra with Applications, 20(4):663–688, 2013.
  * [13] Daniele Colombo, Ersan Turkoglu, Ernesto Sandoval-Curiel, and Taqi Alyousuf.  Self-supervised, active learning seismic full waveform inversion.  Geophysics, 89(2):1–93, 2023.
  * [14] E. Soubies, T.-A. Pham, and M. Unser.  Efficient inversion of multiple-scattering model for optical diffraction tomography.  Optics express, 25:21786–21800, 2017.
  * [15] Björn Engquist and Andrew Majda.  Absorbing boundary conditions for numerical simulation of waves.  Proceedings of the National Academy of Sciences, 74(5):1765–1766, 1977.
  * [16] Bjorn Engquist and Andrew Majda.  Radiation boundary conditions for acoustic and elastic wave calculations.  Communications on pure and applied mathematics, 32(3):313–357, 1979.
  * [17] Yogi A Erlangga, Cornelis W Oosterlee, and Cornelis Vuik.  A novel multigrid based preconditioner for heterogeneous helmholtz problems.  SIAM Journal on Scientific Computing, 27(4):1471–1492, 2006.
  * [18] Yinan Feng, Yinpeng Chen, Peng Jin, Shihang Feng, Zicheng Liu, and Youzuo Lin.  Simplifying full waveform inversion via domain-independent self-supervised learning.  arXiv preprint arXiv:2305.13314, 2023.
  * [19] Samy Wu Fung, Howard Heaton, Qiuwei Li, Daniel McKenzie, Stanley Osher, and Wotao Yin.  JFB: Jacobian-free backpropagation for implicit networks.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 36, pages 6648–6656, 2022.
  * [20] Maayan Gelboim, Amir Adler, Yen Sun, and Mauricio Araya-Polo.  Encoder–decoder architecture for 3d seismic inversion.  Sensors, 23(1):61, 2022.
  * [21] Gelboim, Maayan and Adler, Amir and Sun, Yen and Araya-Polo, Mauricio.  Deep compressed learning for 3d seismic inversion.  In SEG International Exposition and Annual Meeting, pages SEG–2023. SEG, 2023.
  * [22] Eldad Haber and Uri M Ascher.  Preconditioned all-at-once methods for large, sparse parameter estimation problems.  Inverse Problems, 17(6):1847, 2001.
  * [23] Eldad Haber, Matthias Chung, and Felix Herrmann.  An effective method for parameter estimation with pde constraints with multiple right-hand sides.  SIAM Journal on Optimization, 22(3):739–757, 2012.
  * [24] Peng Jin, Xitong Zhang, Yinpeng Chen, Sharon Xiaolei Huang, Zicheng Liu, and Youzuo Lin.  Unsupervised learning of full-waveform inversion: Connecting cnn and partial differential equation in a loop.  arXiv preprint arXiv:2110.07584, 2021.
  * [25] Justin Johnson, Alexandre Alahi, and Li Fei-Fei.  Perceptual losses for real-time style transfer and super-resolution.  European Conference on Computer Vision, pages 694–711, 2016.
  * [26] Yuji Kim and Nori Nakata.  Geophysical inversion versus machine learning in inverse problems.  The Leading Edge, 37(12):894–901, 2018.
  * [27] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [28] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton.  Imagenet classification with deep convolutional neural networks.  Communications of the ACM, 60(6):84–90, 2017.
  * [29] L. Guasch, O. C. Agudo, M.-X. Tang, P. Nachev, and M. Warner.  Full-waveform inversion imaging of the human brain.  NPJ digital medicine, pages 1–12, 2020.
  * [30] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep learning.  nature, 521(7553):436–444, 2015.
  * [31] Bar Lerer, Ido Ben-Yair, and Eran Treister.  Multigrid-augmented deep learning preconditioners for the helmholtz equation using compact implicit layers, 2024.
  * [32] Shucai Li, Bin Liu, Yuxiao Ren, Yangkang Chen, Senlin Yang, Yunhai Wang, and Peng Jiang.  Deep-learning inversion of seismic data.  IEEE Transactions on Geoscience and Remote Sensing, 58(3):2135–2149, 2020.
  * [33] Yunyue Elita Li and Laurent Demanet.  Full-waveform inversion with extrapolated low-frequency data.  Geophysics, 81(6):R339–R348, 2016.
  * [34] Dong Liang, Jing Cheng, Ziwen Ke, and Leslie Ying.  Deep magnetic resonance image reconstruction: Inverse problems meet neural networks.  IEEE Signal Processing Magazine, 37(1):141–151, 2020.
  * [35] Michelle Liu, Rajiv Kumar, Eldad Haber, and Aleksandr Aravkin.  Simultaneous-shot inversion for pde-constrained optimization problems with missing data.  Inverse Problems, 35(2):025003, 2018.
  * [36] Alice Lucas, Michael Iliadis, Rafael Molina, and Aggelos K Katsaggelos.  Using deep neural networks for inverse problems in imaging: beyond analytical methods.  IEEE Signal Processing Magazine, 35(1):20–36, 2018.
  * [37] Ludovic Métivier, Romain Brossier, Stéphane Operto, and Jean Virieux.  Full waveform inversion and the truncated newton method.  SIAM review, 59(1):153–195, 2017.
  * [38] Ana PO Muller, Jessé C Costa, Clecio R Bom, Matheus Klatt, Elisangela L Faria, Marcelo P de Albuquerque, and Marcio P de Albuquerque.  Deep pre-trained fwi: where supervised learning meets the physics-informed neural networks.  Geophysical Journal International, 235(1):119–134, 2023.
  * [39] Oleg Ovcharenko, Vladimir Kazei, Mahesh Kalita, Daniel Peter, and Tariq Alkhalifah.  Deep learning for low-frequency extrapolation from multioffset seismic data.  Geophysics, 84(6):R989–R1001, 2019.
  * [40] Min Jun Park and Mauricio D Sacchi.  Automatic velocity analysis using convolutional neural network and transfer learning.  Geophysics, 85(1):V33–V43, 2020.
  * [41] R Gerhard Pratt.  Seismic waveform inversion in the frequency domain, part 1: Theory and verification in a physical scale model.  Geophysics, 64(3):888–901, 1999.
  * [42] R Gerhard Pratt, Changsoo Shin, and GJ Hick.  Gauss–newton and full newton methods in frequency–space seismic waveform inversion.  Geophysical journal international, 133(2):341–362, 1998.
  * [43] Maziar Raissi, Paris Perdikaris, and George E Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational Physics, 378:686–707, 2019.
  * [44] Majid Rasht-Behesht, Christian Huber, Khemraj Shukla, and George Em Karniadakis.  Physics-informed neural networks (pinns) for wave propagation and full waveform inversions.  Journal of Geophysical Research: Solid Earth, 127(5):e2021JB023120, 2022.
  * [45] Alan Richardson.  Seismic full-waveform inversion using deep learning tools and techniques.  arXiv preprint arXiv:1801.07232, 2018.
  * [46] Leonid I Rudin, Stanley Osher, and Emad Fatemi.  Nonlinear total variation based noise removal algorithms.  Physica D: nonlinear phenomena, 60(1-4):259–268, 1992.
  * [47] S. Bernard, V. Monteiller, D. Komatitsch, and P. Lasaygues.  Ultrasonic computed tomography based on full-waveform inversion for bone quantitative imaging.  Physics in Medicine & Biology, 62, 2017.
  * [48] Yash Sanghvi, Yaswanth Kalepu, and Uday K Khankhoje.  Embedding deep learning in inverse scattering problems.  IEEE Transactions on Computational Imaging, 6:46–56, 2019.
  * [49] Chao Song, Tariq Alkhalifah, and Umair Bin Waheed.  Solving the frequency-domain acoustic vti wave equation using physics-informed neural networks.  Geophysical Journal International, 225(2):846–859, 2021.
  * [50] Song, Chao and Alkhalifah, Tariq and Waheed, Umair Bin.  A versatile framework to solve the helmholtz equation using physics-informed neural networks.  Geophysical Journal International, 228(3):1750–1762, 2022.
  * [51] Antonio Stanziola, Simon R Arridge, Ben T Cox, and Bradley E Treeby.  A helmholtz equation solver using unsupervised learning: Application to transcranial ultrasound.  Journal of computational physics, 441:110430, 2021.
  * [52] Hongyu Sun and Laurent Demanet.  Extrapolated full-waveform inversion with deep learning.  Geophysics, 85(3):R275–R288, 2020.
  * [53] Jian Sun, Kristopher Innanen, Tianze Zhang, and Daniel Trad.  Implicit seismic full waveform inversion with deep neural representation.  Journal of Geophysical Research: Solid Earth, 128(3):e2022JB025964, 2023.
  * [54] Christian Szegedy, Wei Liu, Yangqing Jia, Pierre Sermanet, Scott Reed, Dragomir Anguelov, Dumitru Erhan, Vincent Vanhoucke, and Andrew Rabinovich.  Going deeper with convolutions.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 1–9, 2015.
  * [55] Fabricio Arend Torres, Marcello Massimo Negri, Monika Nagy-Huber, Maxim Samarin, and Volker Roth.  Mesh-free eulerian physics-informed neural networks.  arXiv preprint arXiv:2206.01545, 2022.
  * [56] Eran Treister and Eldad Haber.  Full waveform inversion guided by travel time tomography.  SIAM Journal on Scientific Computing, 39(5):S587–S609, 2017.
  * [57] Jean Virieux and Stéphane Operto.  An overview of full-waveform inversion in exploration geophysics.  Geophysics, 74(6):WCC1–WCC26, 2009.
  * [58] UB Waheed, Ehsan Haghighat, Tariq Alkhalifah, Chao Song, and Qi Hao.  Eikonal solution using physics-informed neural networks.  In EAGE 2020 annual conference & exhibition online, pages 1–5. European Association of Geoscientists & Engineers, 2020.
  * [59] R Wang, C Bao, and L Qiu.  Seismic waveform inversion with source manipulation.  In 82nd EAGE Annual Conference & Exhibition, pages 1–5. European Association of Geoscientists & Engineers, 2021.
  * [60] Yue Wu and Youzuo Lin.  InversionNet: An efficient and accurate data-driven full waveform inversion.  IEEE Transactions on Computational Imaging, 6:419–433, 2019.


# Accurate Hyperfine Tensors for Solid State Quantum Applications: Case of the
NV Center in Diamond

István Takács  Department of Physics of Complex Systems, Eötvös Loránd
University, Egyetem tér 1-3, H-1053 Budapest, Hungary  MTA–ELTE Lendület
”Momentum” NewQubit Research Group, Pázmány Péter, Sétány 1/A, 1117 Budapest,
Hungary  Viktor Ivády<EMAIL_ADDRESS>Department of Physics of
Complex Systems, Eötvös Loránd University, Egyetem tér 1-3, H-1053 Budapest,
Hungary  MTA–ELTE Lendület ”Momentum” NewQubit Research Group, Pázmány Péter,
Sétány 1/A, 1117 Budapest, Hungary  Department of Physics, Chemistry and
Biology, Linköping University, SE-581 83 Linköping, Sweden

(August 30, 2024)

###### Abstract

The decoherence of point defect qubits is often governed by the electron spin-
nuclear spin hyperfine interaction that can be parameterized by using ab inito
calculations in principle. So far most of the theoretical works have focused
on the hyperfine interaction of the closest nuclear spins, while the accuracy
of the predictions for distinct nuclear spins is barely discussed. We
demonstrate for the case of the NV center in diamond that the absolute
relative error of the computed hyperfine parameters can exceed 100% in VASP
for weakly coupled nuclear spins. To overcome this issue, we implement an
alternative method and report on significantly improved hyperfine values with
$\mathcal{O}$(1%) relative mean error at all distances. The provided accurate
hyperfine data for the NV center enables high-precision simulation of NV
quantum nodes for quantum information processing and positioning of nuclear
spins by comparing experimental and theoretical hyperfine data.

Point defects have been widely used to control the optical and electronic
properties of semiconductors. Recently, the magnetic properties of these
materials have also been tailored by paramagnetic defects giving rise to
various microscopic and mesoscopic magnetic phenomena. At low defect
counteractions, controllable few-spin systems can be realized that has led to
the development of point defect quantum bitsJelezko  _et al._ (2001); Weber
_et al._ (2010) (qubits) and quantum nodesWehner  _et al._ (2018); Bradley
_et al._ ; Abobeih  _et al._ . In contrast to other qubit implementations,
point defect qubits in wide-bandgap semiconductors are highly coherent and
robust even at elevated temperatures.Weber  _et al._ (2010); Wolfowicz  _et
al._ (2021) Such optically addressable spin qubits, realized for example by
the NV center in diamondDoherty  _et al._ (2013), the silicon vacancyWidmann
_et al._ (2015); Ivády  _et al._ (2017) in silicon carbide (SiC), and
divacancy related defects in SiCKoehl  _et al._ (2011); Ivády  _et al._
(2019), can possess as long coherence time as 1 ms at room
temperature.Balasubramanian  _et al._ (2009); Christle  _et al._ (2015);
Widmann  _et al._ (2015) Research activities in this area have gained
considerable momentum over the last decades and point defect-based quantum
devices have become leading contenders in several areas of quantum
technologies, such as quantum sensing and quantum internet.Awschalom  _et al._
(2018)

The coherence of spin qubits in light element semiconductors is often limited
by spin-spin interactions with paramagnetic defects and nuclei. In high-purity
samples, the magnetic environment of a spin qubit is defined by the
surrounding nuclear spin bath.Zhao  _et al._ (2012); Seo  _et al._ (2016)
Point defect spins interact with nuclear spins through the hyperfine coupling
that depends on the spatial distribution of the defect’s spin density and the
position of the nuclear spins. The hyperfine spin Hamiltonian term is
parameterized by the hyperfine tensor, whose elements can be measured by
various magnetic resonance techniques and calculated by using first-principles
electronic structure methods. Conventionally, electron spin resonance (ESR)
has been used to determine the hyperfine tensor for the closest nuclear spins
giving rise to $\sim$10-300 MHz hyperfine splitting of the nuclear spin
sublevels. The high controllability and the long coherence time of point
defect qubits have enabled more sophisticated nuclear spin detection
techniques to be developed. Optically detected magnetic resonance (ODMR)
measurement of individual point defect qubits allowed the detection of nuclear
spins 2.5-7 Å distances from the NV center with hyperfine splitting ranging
from 430 kHz to 14 MHz.Smeltzer  _et al._ ; Dréau  _et al._ (2012) Dynamic
decoupling techniques can be used to boost further the sensitivity of the
measurements.Taminiau  _et al._ (2012); Wang  _et al._ (2016); Zopes  _et al._
(2018); Abobeih  _et al._ (2019); van de Stolpe  _et al._ (2023) Using such
techniques, nuclear spins as distant as 30 Å could be detected with hyperfine
splitting of $\sim 1$ kHz.Abobeih  _et al._ (2019); van de Stolpe  _et al._
(2023) These developments have opened new directions in magnetic resonance and
magnetic resonance imaging in nanometer scales.Taminiau  _et al._ (2012); Wang
_et al._ (2016); Zopes  _et al._ (2018); Abobeih  _et al._ (2019); van de
Stolpe  _et al._ (2023) In addition, the characterized nearby nuclear spins
can be utilized as additional highly coherent quantum resources for quantum
computation and quantum internet.Bradley  _et al._ ; Abobeih  _et al._

Hyperfine coupling tensors can be calculated using different first principles
methods, such as density functional theory (DFT) and wave functions-based
methods. The accuracy of the computed hyperfine tensors for close nuclear
spins is remarkable. For example, considering paramagnetic point defects in
semiconductors a mean absolute relative error of 4.7% has been reported.Szász
_et al._ (2013) Since the hyperfine structure of the point defect’s spin
sublevels is unique like a bar code, one can compare the measured and computed
hyperfine tensors to unambiguously identify paramagnetic defects in
semiconductors, see for instance Refs. Son  _et al._ ; Isoya  _et al._ ; Ivády
_et al._ (2017).

As demonstrated recently by a supercell-size-scaling test in Ref. Ivády  _et
al._ (2018), the remarkable accuracy of the computed hyperfine parameters is
only limited for the closest nuclear spin in 1-5 Å distances from the defect.
The relative error sharply increases for nuclear spins located at large
distances. As discussed in Ref. Ivády  _et al._ (2018), this is presumably due
to the periodic boundary condition and related finite size effects. Using a
finite cluster model, such as the C510H252 cluster used in Ref. Nizovtsev  _et
al._ to calculate the hyperfine coupling tensors at 1.5-8 Å distances from the
defect, errors related to the periodic boundary condition can be eliminated.
Another limitation of first principles calculations is the model size that
maximizes the number of lattice sites that can be considered.

In this Letter, we first demonstrate the inaccuracy of the numerical hyperfine
parameters obtained with the industry standard VASP codeKresse and Hafner
(1994); Kresse and Furthmüller (1996). To resolve the underlying
methodological issues we introduce a real-space integration method and the use
of a large support lattice for considering nuclear spins outside the
boundaries of the supercell. To benchmark our method, we carry out large-scale
calculations for the NV center in diamond using different exchange-correlation
functionals and compare the numerical results with available experimental
hyperfine data sets. From the comparison, we conclude that the HSE06Heyd  _et
al._ (2003, 2006) functional with 0.2 mixing parameter performs best for the
NV center in diamond resulting in a mean absolute percentage error of 1.7% for
nuclear spins 6-30 nm distances from the NV center. This is a significant
improvement compared to previous theoretical predictions obtained by using
VASP. We show that the residual errors are likely related to the inaccurate
calculation of the Fermi contact term. High accuracy hyperfine tensors of
$\approx 10^{4}$ lattice sites as well as volumetric hyperfine data with
$<0.1$ Å spatial resolution are published together with this Letter at Ref.
111Theoretical hyperfine tensors for all lattice sites within 30 Å distance
from the NV center are available at https://ivadygroup.elte.hu/hyperfine and
ready to be used for modeling NV center quantum nodes and positioning nuclear
spins around the NV center in diamond in nano-NMR measurements.

Figure 1: Spin density (red lobs) of the NV center viewed from [111]
direction.

The hyperfine interaction describes the weak coupling between the electron and
the nuclear spins. Integrating out the spatial degrees of freedom, the
corresponding hyperfine spin Hamiltonian term can be written as

$H_{\text{HF}}=\mathbf{S}\mathbf{\mathcal{A}}\mathbf{I}\text{,}$ (1)

where $\mathbf{\mathcal{A}}$ is the hyperfine tensor and $\mathbf{S}$ and
$\mathbf{I}$ are the electron spin and the nuclear spin operator vectors with
$S=1$ and $I=1/2$ quantum numbers for the NV center and adjacent 13C nuclear
spins of the diamond lattice, respectively. The two dominant contributions to
the hyperfine tensor are the Fermi contact interaction,
$\mathcal{A}_{\text{FC}}$, and the magnetic dipole-dipole coupling of the
electron and nuclear spins, $\mathcal{A}_{\text{SS}}$.

Elements of the hyperfine tensor can be calculated given the spin density of
the electron spin $\sigma(\mathbf{r})$, see Fig. 1 for the NV center, and the
position $\mathbf{R}_{J}$ of nuclear spin $J$ as Szász  _et al._ (2013)

$\mathcal{A}_{ij}^{J}=\frac{1}{2S}\gamma_{J}\gamma_{e}\hbar^{2}\Bigg{[}\frac{8\pi}{3}\int\delta(\mathbf{r}-\mathbf{R}_{J})\sigma(\mathbf{r})d\mathbf{r}+\mathcal{W}_{ij}(\mathbf{R}_{J})\Bigg{]}\text{
,}$ (2)

where $\gamma_{J}$ and $\gamma_{e}$ are the gyromagnetic ratio of nucleus $J$
and the electron, respectively. The first term on the r.h.s. of Eq. (2)
accounts for the Fermi contact interactions, where the Dirac delta function
$\delta(\mathbf{r}-\mathbf{R}_{J})$ takes the value of the spin density at the
position of the nucleus whose spatial distribution is neglected. The largest
contribution to the Fermi contact term originates from electronic states of
$s$ orbit character exhibiting non-zero probability at the nucleus site. The
second term on the r.h.s. of Eq. (2) accounts for the the dipole-dipole
interaction, where the integral $\mathcal{W}_{ij}$ can be expressed as

$\mathcal{W}_{ij}(\mathbf{R})=\int\Bigg{(}\frac{3(\mathbf{r}-\mathbf{R})_{i}(\mathbf{r}\mathbf{R})_{j}}{|\mathbf{r}-\mathbf{R}|^{5}}-\frac{\delta_{ij}}{|\mathbf{r}-\mathbf{R}|^{3}}\Bigg{)}\sigma(\mathbf{r})d\mathbf{r}\text{.}$
(3)

First principles electronic structure codes for solid state physics often use
periodic boundary conditions, plane wave basis sets, and pseudopotentials to
describe valance states of periodic lattices. For example, the industry
standard VASP software package Kresse and Hafner (1994); Kresse and
Furthmüller (1996) calculates the hyperfine tensor using the method of P. E.
Blöchl Blöchl (2000) while taking core polarization contributions into account
Szász  _et al._ (2013) in periodic boundary conditions. Employing the
projector augmented wave (PAW) method Blöchl (1994), the total spin density of
the system is composed of three partsBlöchl (2000); Szász  _et al._ (2013)

$\sigma=\tilde{\sigma}+\sigma^{1}-\tilde{\sigma}^{1}\text{,}$ (4)

where $\sigma^{1}$ and $\tilde{\sigma}^{1}$ are the atomic core-centered true
and pseudo-spin densities, respectively, and $\tilde{\sigma}$ is the total
spin density of the valance electrons calculated using pseudo-potentials. With
this differentiation, the Fermi contact interaction can be written as

$\displaystyle\frac{8\pi}{3}\int\delta(\mathbf{r}-\mathbf{R_{J}})\sigma(\mathbf{r})d\mathbf{r}=\frac{8\pi}{3}\Bigg{[}\sum_{\mathbf{G}}\tilde{\sigma}(\mathbf{G})e^{i\mathbf{G}\mathbf{R}}+$
(5)
$\displaystyle+\int\delta_{T}(r)\sigma^{1}_{s\mathbf{R}}(r)dr-\int\delta_{T}(r)\tilde{\sigma}^{1}_{s\mathbf{R}}(r)dr\Bigg{]}\text{,}$
(6)

where $\sigma^{1}_{s\mathbf{R}}$ and $\tilde{\sigma}^{1}_{s\mathbf{R}}$ are
the $s$-like contributions to the true and pseudo-core-centered spin
densities, respectively. Spin polarization of the core electrons can be
calculated within the frozen valence approximation.Yazyev  _et al._ (2005);
Szász  _et al._ (2013) The computed core polarization is added to $\sigma^{1}$
and the corresponding hyperfine contribution $\mathcal{A}_{\text{1c}}$ is
determined. With a similar line of thought, the dipole-dipole interaction can
also be expressed as

$\mathcal{W}_{ij}(\mathbf{R})=\tilde{\mathcal{W}}_{ij}(\mathbf{R})+\mathcal{W}_{ij}^{1}(\mathbf{R})-\tilde{\mathcal{W}}_{ij}^{1}(\mathbf{R})\text{,}$
(7)

where the valance electrons’ contribution $\tilde{\mathcal{W}}_{ij}$ is
obtained from the pseudo-spin density $\tilde{\sigma}$ as

$\tilde{\mathcal{W}}_{ij}(\mathbf{R})=-4\pi\sum_{\mathbf{G}}\Bigg{(}\frac{G_{i}G_{j}}{G^{2}}-\frac{\delta_{ij}}{3}\Bigg{)}\tilde{\sigma}(\mathbf{G})e^{i\mathbf{G}\mathbf{R}}\text{,}$
(8)

and the one-center contributions to the dipole-dipole term can be obtained
from $d$-like contribution to the one-center spin density.Blöchl (2000); Szász
_et al._ (2013)

First, we use this method as implemented in VASP to calculate hyperfine
tensors of all sites in a 512-atom and a 1728-atom supercell of diamond
including a single NV center in the middle of the supercell, see Fig. 1. For
the calculations, we use HSE06 exchange-correlation functional, 500 eV cutoff
energy of the plane-wave basis set, $\Gamma$-point sampling of the Brillouin
zone, and high convergence criteria. The structure of the defect is optimized
as far as the largest force is smaller than 10-3 eV/Å. To compare our results
with the experimental values, see Fig. 2(a), we either compute the hyperfine
splitting as $A_{z}=\sqrt{A_{xz}^{2}+A_{yz}^{2}+A_{zz}^{2}}$ Dréau  _et al._
(2012), where $A_{ij}$ are the elements of the hyperfine tensor, which is
compared with experimental values in data set I. taken from Ref. Dréau  _et
al._ (2012), or compare the $A_{zz}$ hyperfine tensor element with
experimental $A_{zz}$ values obtained from high precision measurements in Ref.
Taminiau  _et al._ (2012) (data set II.) and Ref. van de Stolpe  _et al._
(2023) (data set III.). Note that the latter data set is a refined and
extended version of the data published in Ref. Abobeih  _et al._ (2019). From
the reported 50 nuclear spins, we consider only those for which both the
position and the hyperfine parameter are reported with high accuracy, see Fig.
2(a).

Figure 2: Comparison of experimental and theoretical hyperfine parameters. a)
Experimental data sets of hyperfine parameters reported in Ref. Dréau  _et
al._ (2012) (data set I.), in Ref. Taminiau  _et al._ (2012) (data set II.),
and in Ref. van de Stolpe  _et al._ (2023) (data set III.). Gray columns
depict the measured hyperfine parameters, $A_{z}$ for data set I. and $A_{zz}$
for data sets II. and III. The red bars depict the absolute error of the
measurements. b) Absolute relative error of the calculated hyperfine values
showing increasing error due to finite-size effects. Value 1 on the $y$-scale
corresponds to 100% absolute percentage error. When no bars are depicted for
certain elements of the data sets, the corresponding nuclei site is outside
the supercell used for the calculations.

In Fig. 2(b), we depict the absolute relative error of the computed hyperfine
splitting values for data set I.-III. As can be seen, the absolute relative
error of the theoretical values rapidly increases with decreasing value of the
experimental hyperfine parameter and below $\sim$1 MHz wildly fluctuates. It
is also clear that the use of a large supercell does not improve the
comparison with the experiments. Finally, we note that for making the
comparison between theory and experiment we needed to position the nuclear
spins measured in data sets I. and II. For this, we used our improved
hyperfine tensor computation discussed next.

Finite-size correction of hyperfine tensors computed in periodic supercell
models has not been thoroughly investigated before. There are several possible
sources of finite-size effects that can lead to non-physical interaction of
nuclear spins and period defect structures. For instance, the spin density of
a defect and its periodic replicas may overlap in small supercell models
giving rise to overestimated Fermi contact interaction terms and perturbed
dipole-dipole interaction terms for certain nuclei sites. This error is,
however, assumed to decay exponentially with the size of the supercell since
localized defect states decay exponentially. More difficult-to-handle finite-
size effects arise from the long-range dipole-dipole interaction that decays
with the third power of the distance of the spins. In periodic supercell
models, nuclear spins interact with a lattice of defects. In contrast to the
Coulomb interaction, the dipole-dipole interaction does not diverge, although
the interaction strength and the hyperfine tensor’s principal axis may be
considerably perturbed. A nuclear spin halfway between the defect and one of
its periodic replicas interacts with two spin densities with approximately the
same coupling strength that can give rise to $\mathcal{O}(100\%)$ error
explaining our observations depicted in Fig. 2(b). In the following, we remedy
these finite-size effects.

The errors are derived from the long-ranged interaction term, i.e. the pseudo-
dipole-dipole integral $\tilde{\mathcal{W}}_{ij}$ defined in Eq. (8). The
periodicity of the spin density and thus the finite size effects are encoded
into the pseudo spin density $\tilde{\sigma}(\mathbf{G})$. The point defect
and its replicas cannot be separated in Fourier space; however, in real space,
this can be easily achieved by simply limiting the range of integration. To
overcome the finite-size dependence of the dipolar hyperfine interaction term,
we utilize this strategy and combine it with the PAW method. To this end, we
calculate $\tilde{\mathcal{W}}_{ij}$ as

$\tilde{\mathcal{W}}_{ij}(\mathbf{R^{+}})=\\!\\!\\!\bigintss\displaylimits^{\overline{\Omega}_{\text{SC}}}\\!\\!\\!\\!\Bigg{(}\frac{3(\mathbf{r}-\mathbf{R^{+}})_{i}(\mathbf{r}-\mathbf{R^{+}})_{j}}{|\mathbf{r}-\mathbf{R^{+}}|^{5}}-\frac{\delta_{ij}}{|\mathbf{r}-\mathbf{R^{+}}|^{3}}\Bigg{)}\sigma(\mathbf{r})d\mathbf{r}\text{,}$
(9)

where the extended position coordinates $\mathbf{R^{+}}$ include the lattice
sites within the supercell ($\mathbf{R}$) and outside the supercell provided
by a support lattice within a sphere of 30 Å radius centered on the NV center.
The support lattice is aligned with the supercell, although it does not
contain the atomic sites of the supercell. This way the hyperfine interaction
calculation is not limited by the size of the supercell. Note that for the
integration we use the full spin density $\sigma(\mathbf{r})$ expressed on a
fine grid and not the pseudo spin density in contrast to Eq. (8). For nuclear
spins contained within the supercell, the real-space integration is carried
out over the supercell except a sphere of $r_{\text{PAW}}$ radius centered at
the nuclear spin position $\mathbf{R}$, i.e.
$\overline{\Omega}_{\text{SC}}=\Omega_{\text{SC}}-\Omega_{\text{PAW}}(\mathbf{R})$.
The total dipole-dipole interaction term is defined as
$\mathcal{W}_{ij}(\mathbf{R})=\tilde{\mathcal{W}}_{ij}(\mathbf{R})+\mathcal{W}_{ij}^{1}(\mathbf{R})$
in contrast to Eq. (7). The Fermi contact interaction with a core polarization
contribution is obtained using the VASP’s implementationSzász  _et al._
(2013). For lattice sites outside the supercell’s boundaries, the spin density
is considered to be zero at the nuclear spin site, i.e. both the
$\mathcal{W}_{ij}^{1}(\mathbf{R^{+}})$ and the Fermi contact terms are
approximated to be zero. The sole non-zero term for these sites is given by
Eq. (9), where the integration is carried over the full supercell volume. With
these modifications, the computed tensors account for the case when the
nuclear spins interact with an isolated defect and not with a lattice of
defects. We anticipate that the leading finite-size effects will be removed
within our approach.

To compute the hyperfine tensors for a large number of lattice sites,
$\sim$20000 sites within a sphere $R^{+}_{\text{cut}}=30$ Å around the NV
center, we implement the real space integration method in an in-house code
that post-processes the VASP output files. The ground state calculations of
the NV center are carried out in 512-atom and 1728-atom supercells using the
experimental lattice parameter of 3.567 Å. We use both the semi-local PBE
Perdew  _et al._ (1996) and various forms of the HSE06 hybrid functional,
which are labeled as HSE($\alpha$), where $\alpha$ is the mixing parameters,
e.g. HSE(0.25) = HSE06. The ground state spin density used in the calculation
of Eq. (9) is defined on a real space grid of 0.036 Å spacing. The positions
of the nuclear spins for data set III. are taken from Ref. van de Stolpe  _et
al._ (2023). For data sets I. and II. we use the following strategy for
positioning the nuclear spins: Considering an experimental hyperfine value, we
look for the closest theoretical value in our data set. Depending on the error
bar of the measurement and the calculations, we can position nuclear spin up
to symmetrically equivalent sites with this method.

Figure 3: Absolute relative error (ARE) of the hyperfine splitting values
calculated with our improved integration method. Gray columns depict the ARE
of the experimental data, red thin columns with circles depict the ARE of the
computed hyperfine values obtained with HSE06 functional, and light blue line
with squares depict the ARE of the theoretical values obtained with PBE
functional.

The hyperfine splitting values computed with our method are significantly
improved compared to the currently available implementationSzász  _et al._
(2013). The relative error of the theoretical values reduce from
$\mathcal{O}(100\%)$ to $\mathcal{O}(1\%)$, see Fig. 2(b) and Fig. 3. For
different data sets, we obtain different mean absolute percentage errors
(MAPE). For data set I. we obtain a MAPE of 3.8%, while the experimental data
exhibit an averaged relative error margin of 1.2%Dréau  _et al._ (2012). The
obtained MAPE of the improved theoretical values for data set II. is 1.5%,
which is better than the relative error margin of the experimental data of
3.3%Taminiau  _et al._ (2012). This suggests that the theoretical values are
”overfitted” for data set II., i.e. multiple matching hyperfine values were
found within the error margin of the experimental data. By selecting the
closest ones, we could obtain a MAPE smaller than the experimental error bar.
The most accurate hyperfine values, with a relative error margin of $6\times
10^{-3}\%$, together with nuclear spin positions are provided for data set
III. in Ref van de Stolpe  _et al._ (2023). These high-precision measurements
allow us to make a reliable assessment of the error bar of the theoretical
values. Considering the 29 most accurately measured and positioned nuclear
spins, we obtain a MAPE of 1.79%. In Fig. (3), we also depicted the ARE of
hyperfine values obtained with the PBE functional. As can be seen, the
HSE(0.25) functional consistently improves on the PBE values, especially for
nuclear spins found close to the NV center. We also note that the use of
accurate DFT spin density has high relevance. Considering point spin density
approximation, i.e. $\sigma(\mathbf{r})=\delta(\mathbf{r}-\mathbf{r}_{0})$
where $\mathbf{r}_{0}$ is the center of the NV center, we obtain a MARE of 76%
for data set III. see Supplementary Information.

Figure 4:  Absolute relative error of the theoretical values as a function of
the Fermi contact contribution to the computed hyperfine parameter. The
theoretical values are obtained by using HSE06 functional and 1728 atom
supercell. Green, blue, and red columns depict the ARE values corresponding to
data sets I., II., and III., respectively.

Finally, we investigate the source of residual errors and possible further
measures to improve the computed hyperfine tensors. For this study, we use
data set III. First of all, it is notable that the mean signed relative error
(MSRE) of -0.12% is significantly smaller than the MARE of 1.79%, suggesting
that there are no large systematic errors and the discrepancies may be related
to numeral uncertainties. Next, we plot the MARE as a function of the Fermi
contact term, see Fig. (4), and as a function of the distance, see
Supplementary Information. There seems to be a correlation between the largest
errors and the value of the Fermi contact term. Due to the distance dependence
of the Fermi contact term, the MARE seems to decay with the distance of the
nuclear spins, see Supplementary Information. It should be noted; however,
that large errors are also obtained for a few nuclear spins that are beyond
the boundaries of the supercell, where we explicitly neglect the Fermi contact
terms. Here, the neglect of the Fermi contact term may be the source of the
error. To attempt reducing the errors, we tune the mixing parameter $\alpha$
of the HSE functional and study the variation of the MARE obtained for data
set III., see Supplementary Information. By reducing the mixing parameter from
0.25 to 0.2, we obtain a slightly decreased MARE of 1.69%, although an
enlarged MSRE of -0.32%. Overall, it seems like the hyperfine values may be
further improved enhancing numerical accuracy for the core/Fermi contact
contribution, for example using all-electron methods accounting for the core
electrons and larger supercells.

In summary, we demonstrated in this Letter that high-accuracy, finite-size
effect-free hyperfine tensors can be calculated using an improved integration
method. Compared to the industry standard VASP code, we could achieve a
$\sim$100 fold reduction of the MARE of the theoretical values. Having more
experimental data of high precision and further improved numerical accuracy in
larger supercells may help to obtain superior theoretical hyperfine data
compared to what has been presented here. The obtained and potentially
improved future data is available online under Ref. [33]. The provided data
can be used for high-precision simulation of the NV center-nuclear spin few-
body quantum systems as well as positioning nuclear spin around the NV center.

## Acknowledgments

This work was supported by the National Research, Development and Innovation
Office of Hungary (NKFIH) within the Quantum Information National Laboratory
of Hungary (Grant No. 2022-2.1.1-NL-2022-00004) and within projects FK 135496
and FK 145395. V.I. acknowledges support from the Knut and Alice Wallenberg
Foundation through the WBSQD2 project (Grant No. 2018.0071). The computations
were enabled by resources provided by the National Academic Infrastructure for
Supercomputing in Sweden (NAISS) at the Swedish National Infrastructure for
Computing (SNIC) at Tetralith, partially funded by the Swedish Research
Council through grant agreement no. 2022-06725.

## References

  * Jelezko  _et al._ (2001) F. Jelezko, C. Tietz, A. Gruber, I. Popa, A. Nizovtsev, S. Kilin, and J. Wrachtrup, Single Molecules 2, 255 (2001).
  * Weber  _et al._ (2010) J. R. Weber, W. F. Koehl, J. B. Varley, A. Janotti, B. B. Buckley, C. G. Van de Walle, and D. D. Awschalom, Proc. Natl. Acad. Sci. 107, 8513 (2010).
  * Wehner  _et al._ (2018) S. Wehner, D. Elkouss, and R. Hanson, Science 362 (2018), 10.1126/science.aam9288.
  * (4) C. Bradley, J. Randall, M. Abobeih, R. Berrevoets, M. Degen, M. Bakker, M. Markham, D. Twitchen, and T. Taminiau, Physical Review X 9, 031045.
  * (5) M. H. Abobeih, Y. Wang, J. Randall, S. J. H. Loenen, C. E. Bradley, M. Markham, D. J. Twitchen, B. M. Terhal, and T. H. Taminiau, Nature 606, 884.
  * Wolfowicz  _et al._ (2021) G. Wolfowicz, F. J. Heremans, C. P. Anderson, S. Kanai, H. Seo, A. Gali, G. Galli, and D. D. Awschalom, Nature Reviews Materials , 1 (2021).
  * Doherty  _et al._ (2013) M. W. Doherty, N. B. Manson, P. Delaney, F. Jelezko, J. Wrachtrup, and L. C. L. Hollenberg, Physics Reports 528, 1 (2013).
  * Widmann  _et al._ (2015) M. Widmann, S.-Y. Lee, T. Rendler, N. T. Son, H. Fedder, S. Paik, L.-P. Yang, N. Zhao, S. Yang, I. Booker, A. Denisenko, M. Jamali, S. A. Momenzadeh, I. Gerhardt, T. Ohshima, A. Gali, E. Janzén, and J. Wrachtrup, Nature Materials 14, 164 (2015).
  * Ivády  _et al._ (2017) V. Ivády, J. Davidsson, N. T. Son, T. Ohshima, I. A. Abrikosov, and A. Gali, Physical Review B 96, 161114 (2017).
  * Koehl  _et al._ (2011) W. F. Koehl, B. B. Buckley, F. J. Heremans, G. Calusine, and D. D. Awschalom, Nature 479, 84 (2011).
  * Ivády  _et al._ (2019) V. Ivády, J. Davidsson, N. Delegan, A. L. Falk, P. V. Klimov, S. J. Whiteley, S. O. Hruszkewycz, M. V. Holt, F. J. Heremans, N. T. Son, D. D. Awschalom, I. A. Abrikosov, and A. Gali, Nature Communications 10, 5607 (2019).
  * Balasubramanian  _et al._ (2009) G. Balasubramanian, P. Neumann, D. Twitchen, M. Markham, R. Kolesov, N. Mizuochi, J. Isoya, J. Achard, J. Beck, J. Tissler, V. Jacques, P. R. Hemmer, F. Jelezko, and J. Wrachtrup, Nature Materials 8, 383 (2009).
  * Christle  _et al._ (2015) D. J. Christle, A. L. Falk, P. Andrich, P. V. Klimov, J. U. Hassan, N. T. Son, E. Janzén, T. Ohshima, and D. D. Awschalom, Nature Materials 14, 160 (2015).
  * Awschalom  _et al._ (2018) D. D. Awschalom, R. Hanson, J. Wrachtrup, and B. B. Zhou, Nature Photonics 12, 516 (2018).
  * Zhao  _et al._ (2012) N. Zhao, S.-W. Ho, and R.-B. Liu, Physical Review B 85, 115303 (2012).
  * Seo  _et al._ (2016) H. Seo, A. L. Falk, P. V. Klimov, K. C. Miao, G. Galli, and D. D. Awschalom, Nature Communications 7, 12935 (2016).
  * (17) B. Smeltzer, L. Childress, and A. Gali, New Journal of Physics 13, 025021.
  * Dréau  _et al._ (2012) A. Dréau, J.-R. Maze, M. Lesik, J.-F. Roch, and V. Jacques, Physical Review B 85, 134107 (2012).
  * Taminiau  _et al._ (2012) T. H. Taminiau, J. J. T. Wagenaar, T. van der Sar, F. Jelezko, V. V. Dobrovitski, and R. Hanson, Physical Review Letters 109, 137602 (2012).
  * Wang  _et al._ (2016) Z.-Y. Wang, J. F. Haase, J. Casanova, and M. B. Plenio, Phys. Rev. B 93, 174104 (2016).
  * Zopes  _et al._ (2018) J. Zopes, K. Herb, K. S. Cujia, and C. L. Degen, Phys. Rev. Lett. 121, 170801 (2018).
  * Abobeih  _et al._ (2019) M. H. Abobeih, J. Randall, C. E. Bradley, H. P. Bartling, M. A. Bakker, M. J. Degen, M. Markham, D. J. Twitchen, and T. H. Taminiau, Nature 576, 411 (2019).
  * van de Stolpe  _et al._ (2023) G. L. van de Stolpe, D. P. Kwiatkowski, C. E. Bradley, J. Randall, S. A. Breitweiser, L. C. Bassett, M. Markham, D. J. Twitchen, and T. H. Taminiau, “Mapping a 50-spin-qubit network through correlated sensing,”  (2023), arXiv:2307.06939 .
  * Szász  _et al._ (2013) K. Szász, T. Hornos, M. Marsman, and A. Gali, Physical Review B 88, 075202 (2013).
  * (25) N. T. Son, P. Carlsson, J. ul Hassan, E. Janzén, T. Umeda, J. Isoya, A. Gali, M. Bockstedte, N. Morishita, T. Ohshima, and H. Itoh, Phys. Rev. Lett. 96, 055501.
  * (26) J. Isoya, T. Umeda, N. Mizuochi, N. T. Son, E. Janzén, and T. Ohshima, physica status solidi (b) 245, 1298.
  * Ivády  _et al._ (2018) V. Ivády, I. A. Abrikosov, and A. Gali, npj Computational Materials 4, 1 (2018).
  * (28) A. P. Nizovtsev, S. Y. Kilin, A. L. Pushkarchuk, V. A. Pushkarchuk, S. A. Kuten, O. A. Zhikol, S. Schmitt, T. Unden, and F. Jelezko, New Journal of Physics 20, 023022.
  * Kresse and Hafner (1994) G. Kresse and J. Hafner, Phys. Rev. B 49, 14251 (1994).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Heyd  _et al._ (2003) J. Heyd, G. E. Scuseria, and M. Ernzerhof,  118, 8207 (2003).
  * Heyd  _et al._ (2006) J. Heyd, G. E. Scuseria, and M. Ernzerhof,  124, 219906 (2006).
  * Note (1) Theoretical hyperfine tensors for all lattice sites within 30 Å distance from the NV center are available at https://ivadygroup.elte.hu/hyperfine.
  * Blöchl (2000) P. E. Blöchl, Physical Review B 62, 6158 (2000).
  * Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Yazyev  _et al._ (2005) O. V. Yazyev, I. Tavernelli, L. Helm, and U. Röthlisberger, Physical Review B 71, 115110 (2005).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).