# Coherence dynamics in low-energy nuclear fusion

Iain Lee  Alexis Diaz-Torres<EMAIL_ADDRESS>
###### Abstract

Low-energy nuclear fusion reactions have been described using a dynamical
coupled-channels density matrix method, based on the theory of open quantum
systems. For the first time, this has been combined with an energy projection
method, permitting the calculation of energy resolved fusion probabilities.
The results are benchmarked against calculations using stationary Schrödinger
dynamics and show excellent agreement. Calculations of entropy, energy
dissipation and coherence were conducted, demonstrating the capability of this
method. It is evident that the presence of quantum decoherence does not affect
fusion probability. This framework provides a basis for quantum thermodynamic
studies using thermal environments.

###### keywords:

Coupled-channels density matrix , open quantum systems , nuclear fusion ,
quantum tunnelling , entropy , decoherence

††journal: Physics Letters B

organization=Department of Physics, University of Surrey,addressline=Faculty
of Engineering and Physical Sciences, city=Guildford, postcode=GU2 7XH,
country=United Kingdom

##  1 Introduction

In stars, nuclear fusion takes place within a plasma — a hot ionised gas where
atoms are stripped apart into electrons and positively charged nuclei. For
fusion to occur, the Coulomb barrier between a projectile and target must be
overcome or tunnelled through. The Coulomb barrier is formed from the action
of two opposite nucleus-nucleus interactions: the short-range, attractive
nuclear interaction and the long-range, repulsive Coulomb interaction. There
is scarse literature detailing how stellar plasmas affect nuclear fusion
reactions, apart from studies on plasma screening [1]. Interactions on the
boundary of nuclear and atomic physics are theorised to excite the atomic
nucleus, and in a dense plasma environment, these excitations may impact
fusion reaction rates [2, 3, 4]. We present a method for describing low-energy
nuclear fusion of heavy ions, with the capability to include effects from
environments that influence these fusion reactions.

The coupled-channels model is combined with an open quantum systems approach
to model nuclear fusion. An open quantum system involves an environment that
interacts with the quantum system of interest (or reduced system) [5]. This is
necessary to describe physical systems, and it is of great interest across
disciplines [6] since the medium in which these systems evolve is often
neglected in calculations. In low-energy nuclear physics, the theory of open
quantum systems has been applied to understanding various dissipative quantum
phenomena, such as deep-inelastic heavy-ion collisions [7] and heavy-ion
capture reactions [8]. For nuclear fusion reactions, the inclusion of
environments may lead to increased population of the intrinsic nuclear excited
states, loss of quantum coherence (decoherence) and energy dissipation,
therefore affecting the fusion probability [9]. As a result of decoherence and
dissipation to an environment, we cannot work with pure states and instead we
use mixed states. These require the use of a density matrix, instead of a
wavefunction, to propagate the dynamics in time. To do this, we use a coupled-
channels density matrix (CCDM) approach based on the Lindblad master equation
[10]. The CCDM approach was introduced in chemical physics [11] and was later
adapted to the field of nuclear physics, where it was used to investigate
dissipative quantum dynamics in nuclear collisions [12, 13].

A key feature of this work is the ability to compare our dynamical (time-
dependent) CCDM technique with stationary (time-independent) methods, such as
the model implemented in the CCFULL code [14], which uses the time-independent
Schrödinger equation (TISE). We do this by using the window operator [15] as
an energy spectral decomposition method to calculate energy-resolved fusion
probabilities. This has been used in literature on wavefunctions [16, 17, 18],
but for the first time, we have applied the window operator to a final
(asymptotic) density matrix. We also present pioneering work on the dynamics
of quantum coherence and entropy production in nuclear fusion, as the CCDM
approach allows us to study these important concepts of the emerging field of
quantum thermodynamics [19].

For the sake of simplicity and as a test case, we currently consider only a
single internal vibrational excitation mode of the 144Sm target that interacts
with an inert 16O projectile. This is the test reaction in the CCFULL code
[14]. This paper serves as a proof of concept for the novel CCDM technique and
the preliminary results provide a benchmark for this approach compared to the
standard, time-independent coupled-channels method [20, 21].

##  2 Theoretical framework

In the present model calculations, there is a projectile nucleus (16O) and a
target nucleus (144Sm) that are initially separated at a large distance so
that virtually no interactions between them are present (due to nuclear and
Coulombic effects) and they reside in their ground states. The collision
described is central and it is assumed that there is no change in orbital
angular momentum $(L=0)$. A “fusion” environment is used for describing
dissipation out of the two-body direct nuclear reaction (reduced system).

###  2.1 Reduced-system Hamiltonian, fusion environment and initial density
matrix

The reduced-system Hamiltonian in the overall centre-of-mass reference frame,
$\hat{H}_{S}$, is the same as the Hamiltonian implemented in the CCFULL code,
and realistic parameters of the model can be found in Ref. [14]:

$\hat{H}_{S}=\hat{T}(r)+\hat{U}(r)+\hat{h}(\xi)+\hat{V}(r,\xi),$ (1)

where $\hat{T}(r)$ and $\hat{U}(r)$ are the kinetic energy operator and the
monopole, total real interaction potential (nuclear + Coulomb parts) between
the two nuclei in their ground states, respectively. $\hat{h}(\xi)$ is the
intrinsic Hamiltonian of the nuclei, whose eigenstates and eigenenergies are
$\ket{\alpha}$ and $e_{\alpha}$, respectively. The total real coupling
potential, $\hat{V}(r,\xi)$, determines how the radial motion affects the
population of the internal energy spectrum of the interacting nuclei. The
model values of $U(r)$ and $V(r,\xi)$ for the 16O $+$ 144Sm collision are
presented in Fig. 1 to show the position and strength of these potentials.
Plotting $V(r,\xi)$ shows the localisation of the coupling potential between
the ground-state and the $3^{-}$ vibrational excited state ($1.81$ MeV) of
144Sm. This coupling is negligible until a position near the Coulomb barrier
is reached.

Figure 1:  The potentials for the test collision 16O $+$ 144Sm as a function
of the internuclear radius. The thick red solid line is the total, bare
nucleus-nucleus potential $U(r)$ and the dashed line is the total coupling
potential, $V(r,\xi)$, between the ground-state and the $3^{-}$ vibrational
excited state ($1.81$ MeV) of 144Sm. The thin black solid line is the decay
function $\Gamma(r)$, which is a Fermi function that removes positional
probability from the reduced-system density matrix due to compound nucleus
formation.

The fusion channel has often been treated using a short-range imaginary
potential in other coupled-channels fusion models, since it is not explicitly
included in those model calculations [22]. We include the fusion channel in
our calculations, and call it the fusion environment which represents the high
energy-density compound-nucleus states. One can imagine that the fusion
environment acts like a detector, performing a continuous position
measurement, effectively measuring the distance between two nuclei within a
certain region of space where fusion takes place, causing loss of information
and increasing entropy production. In this dynamical theory, we extend the
Hilbert space of the system to include the fusion environment, which is
modelled by an auxiliary state. This idea was first introduced in Ref. [23] to
describe particle decay. The coupling of the reduced system to the fusion
environment is irreversible, and can be described by specific Lindblad
operators related to an imaginary potential, $-i\Gamma(r)$ [24]. The
$\Gamma(r)$ function in Fig. 1 determines the radial strength of the
absorption to the fusion environment. By recording the information that leaves
the reduced system with a fusion environment, we are able to quantify the
irreversible effects on the collision dynamics due to energy dissipation,
entropy production and loss of coherence.

To describe the initial radial motion of the nuclei, we have tested two
different initial wave packets with the same Gaussian envelope but two
different boosts, $\mathcal{B}(k_{0}r)$, namely a plane wave ($e^{-ik_{0}r}$)
and an incoming Coulomb wave, $H^{-}_{L=0}(k_{0}r)$ [25]:

$\psi(r,r_{0},\sigma_{0},k_{0})=\mathcal{N}^{-1}\,\exp[-\frac{(r-r_{0})^{2}}{2\sigma_{0}^{2}}]\mathcal{B}(k_{0}r),$
(2)

where $\mathcal{N}$ is a normalisation constant, $r_{0}$ is the initial,
central position of the wave packet, $r$ is a radial grid position,
$\sigma_{0}$ is the spatial dispersion, and $k_{0}$ is the average wave
number, which depends on the average incident energy $E_{0}$, $r_{0}$ and
$\sigma_{0}$ and is found by solving $E_{0}=\bra{\psi}\hat{H}_{S}\ket{\psi}$.
Eq. (2) with a plane-wave boost is a Gaussian wave packet, while this is a
Coulomb wave packet [26] when a Coulomb-wave boost is used. These wave packets
were used to construct the radial part of the initial density matrix.

The density matrix created from either of the wavefunctions above would be
correct if the target and projectile nuclei were to remain only in their
ground states. However, for every position on the radial grid, the target
and/or projectile could be in an excited state. Therefore the density operator
is a tensor and is formed from a mixed basis of states: the radial position
and the intrinsic energy states, $\ket{r}$ and $\ket{\alpha}$, respectively:

$\hat{\rho}=\sum_{\alpha,\beta,r,s}\ket{r}\ket{\alpha}\rho^{rs}_{\alpha\beta}(t)\bra{\beta}\bra{s}.$
(3)

In Eq. (3), the auxiliary state associated with the fusion environment is also
included. The density operator must be Hermitian and positive semi-definite
because the diagonal elements of the density matrix are the coefficients
$\rho_{\alpha\alpha}^{rr}(t)$, which have the physical meaning of
probabilities. Initially, when the nuclei are far apart, only the matrix
elements of the density matrix related to the ground state of the nuclei are
nonzero:

$\rho^{rs}_{11}(t=0)=\ket{r}\,\bra{s}.$ (4)

###  2.2 Time propagation of the coupled-channels density matrix and the
dynamics of quantum coherence

Inserting Eq. (3) into the Lindblad master equation for the density operator
and projecting this equation onto both the radial grid position and intrinsic
energy states, the equations of motion describing the time evolution of the
coupled-channels density matrix, $\rho^{rs}_{\alpha\beta}(t)$, are obtained
[12, 13]. These are separated into matrix elements for the reduced system:

$\dot{\rho}^{rs}_{\alpha\beta}=-\frac{i}{\hslash}\left[\rho^{rs}_{\alpha\beta}\>(e_{\alpha}-e_{\beta})+\sum^{M}_{t=1}(T^{rt}\>\rho^{ts}_{\alpha\beta}-\rho^{rt}_{\alpha\beta}\>T^{ts})\right.\\\
+\left.\rho^{rs}_{\alpha\beta}\>(U^{rr}-U^{ss})+\sum^{N}_{\mu=1}(V^{rr}_{\alpha\mu}\>\rho^{rs}_{\mu\beta}-\rho^{rs}_{\alpha\mu}\>V^{ss}_{\mu\beta})\right]\\\
+\delta_{\alpha\beta}\sum_{\nu=1}^{B}\sqrt{\Gamma^{rr}_{\alpha\nu}}\,\rho^{rs}_{\nu\nu}\sqrt{\Gamma^{ss}_{\alpha\nu}}-\>\frac{1}{2}\sum_{\nu=1}^{B}(\Gamma^{rr}_{\nu\alpha}+\Gamma^{ss}_{\nu\beta})\rho^{rs}_{\alpha\beta},$
(5)

and matrix elements involving the auxiliary state(s):

$\dot{\rho}^{rs}_{kl}=\delta_{kl}\sum_{\nu=1}^{B}\sqrt{\Gamma^{rr}_{k\nu}}\,\rho^{rs}_{\nu\nu}\sqrt{\Gamma^{ss}_{k\nu}}-\frac{1}{2}\sum_{\nu=1}^{B}(\Gamma^{rr}_{\nu
k}+\Gamma^{ss}_{\nu l})\rho^{rs}_{kl},$ (6)

where either $k$ or $l$ is equal to the index of any auxiliary state and
$\Gamma^{rr}_{\nu\nu}=\sum_{\mu\neq\nu}\Gamma^{rr}_{\mu\nu}$ [13]. The
auxiliary state in Eq. (6) allows one to store the information transferred out
of the reduced system. In the present model calculations, $\Gamma^{rr}_{31}$
and $\Gamma^{rr}_{32}$ describe the irreversible coupling between the reduced
system’s states (denoted by 1 and 2) and the auxiliary state (denoted by 3).
These matrix elements are equal to the function $\Gamma(r)$ in Fig. 1, whilst
$\Gamma^{rr}_{21}=\Gamma^{rr}_{12}=0$. In Eq. (5), the coupling matrix
$V^{rr}_{\alpha\beta}\equiv\bra{\alpha}V(r,\xi)\ket{\beta}$. Eqs. (5) and (6)
are directly integrated using the Faber polynomial integrator [27] with the
initial condition given by Eq. (4) and a time step of $1\times 10^{-22}$ s.
The initial radial centroid of the wave packet was $r_{0}=70$ fm, and its
typical spatial dispersion was $\sigma_{0}=10$ fm. The radial grid
($r=1.5–150$ fm) was evenly spaced with $M=1024$ points. The number of energy
states of the reduced system was $N=2$, and the total number of intrinsic
energy states was $B=3$.

The dynamics of quantum coherence in the reduced system can be described by
the time-dependent purity, $\bar{\mathcal{P}}(t)$, and the von Neumann
entropy, $\mathcal{S}(t)$, of the reduced-system density matrix,
$\rho_{S}(t)\equiv\\{\rho^{rs}_{\alpha\beta}(t)\\}$ from Eq. (5):

$\bar{\mathcal{P}}(t)=\textnormal{Tr}[\bar{\rho}_{S}^{2}(t)],$ (7)

$\mathcal{S}(t)=-\textnormal{Tr}\left\\{\rho_{S}(t)\,\textnormal{ln}\,[\rho_{S}(t)]\right\\}=-\sum_{j}\eta_{j}(t)\,\textnormal{ln}\,[\eta_{j}(t)],$
(8)

where $\bar{\rho}_{S}(t)$ is the diagonal-removed reduced system density
matrix [28, 29] (i.e.,
$\bar{\rho}_{S}(t)\equiv\\{\rho^{rs}_{\alpha\beta}(t)\\}$ with
$\rho^{rr}_{\alpha\alpha}(t)=0$), and $\eta_{j}(t)$ are the eigenvalues of
$\rho_{S}(t)$ [19]. Before using these equations, the trace of $\rho_{S}(t)$
should be normalised to unity. A constant value of purity indicates a
completely pure (coherent) state, showing that the interaction with the
environment has not caused a mixed state. Conversely, a decrease of purity
over time indicates decoherence caused by the effect of the environment.

###  2.3 Energy projection of the time-propagated density matrix for fusion
probability calculations

We model the initial projectile-target radial motion as a wave packet
containing a distribution of energies, with average energy $E_{0}$. To
calculate the fusion probability of the wave packet, the standard approach is
to sum the fusion probability contributions from all energies within the wave
packet, obtaining the total fusion probability for a fusion reaction. However,
this causes the loss of small probabilities from energies well-below the
Coulomb barrier, since the higher energies dominate the fusion probability.
Therefore, a method is needed to resolve the energies of the wave packet and
calculate energy-resolved fusion probabilities to fully understand these low-
energy contributions. To this end, we employ the window operator,
$\hat{\Delta}(E_{k},n,\epsilon)$ [15]:

$\hat{\Delta}(E_{k},n,\epsilon)=\frac{\epsilon^{2^{n}}}{(\hat{H}_{S}-E_{k})^{2^{n}}+\epsilon^{2^{n}}},$
(9)

where $E_{k}$ is the incident energy of interest, $\epsilon$ is an energy
resolution parameter that determines the width of the energy window, and $n$
is a parameter that determines the shape of the energy window. For a large
$n$, the window function is increasingly rectangular.

The amount of probability of the reduced-system density matrix, $\rho_{S}$, in
an energy window of width $2\epsilon$ around $E_{k}$ is:

$\mathcal{P}(E_{k},n,\epsilon)=\textnormal{Tr}\big{[}\hat{\Delta}(E_{k},n,\epsilon)\,\rho_{S}\big{]}.$
(10)

By expanding Eq.(10) with $n=2$, the successive linear equations to obtain an
energy-resolved density matrix, $\rho^{\prime}_{k}$, can be written as:

$(\hat{H}_{S}-E_{k}+i^{\frac{3}{2}}\epsilon)(\hat{H}_{S}-E_{k}-i^{\frac{3}{2}}\epsilon)(\hat{H}_{S}-E_{k}+\sqrt{i}\epsilon)(\hat{H}_{S}-E_{k}-\sqrt{i}\epsilon)\rho^{\prime}_{k}=\epsilon^{4}\,\rho_{S},$
(11)

and Eq. (10) can be rewritten as:

$\mathcal{P}(E_{k},n=2,\epsilon)=\text{Tr}(\rho^{\prime}_{k}).$ (12)

Eqs. (11) and (12) have been used to calculate the fusion probability,
$P_{fus}$, for a range of below-barrier to above-barrier incident energies,
$E_{k}$:

$P_{fus}(E_{k})=1-\frac{\mathcal{P}_{final}(E_{k})}{\mathcal{P}_{initial}(E_{k})},$
(13)

where $\epsilon$ was equal to $0.35$ MeV. In Eq. (13), $\mathcal{P}_{initial}$
and $\mathcal{P}_{final}$ refer to the energy spectrum associated with the
initial and final reduced-system density matrices, respectively. Having
applied the window operator, we can now compare our results of $P_{fus}$ from
density matrix calculations with those from TISE calculations.

##  3 Results

(a) Subfigure 1 list of figures text

(b) Subfigure 2 list of figures text

Figure 2: Radial position probability as a function of internuclear radius and
time for a head-on collision of 16O + 144Sm with a mean energy of $60$ MeV.
The radial probability decreases and increases for (a) the elastic and (b)
inelastic channels respectively, as the nuclei approach their Coulomb barrier
($r\approx 10$ fm). For visualisation, when the mean radius is larger than
$20$ fm, the time step is $3\times 10^{-22}$ s.

(a)

(b)

Figure 3: The dynamics of energy dissipation and quantum coherence in the
collision scenario shown in Fig. 2. (a) Energy loss (thick solid red line) and
entropy production (thin solid black line) increase as the nuclei approach
their Coulomb barrier, and reach a plateau when the nuclei re-separate. The
dotted black line shows the average internuclear radius of the reduced system.
(b) The purity of the reduced-sytem density matrix decreases (decoherence)
when the nuclei are near the Coulomb barrier. There is coherence resurgence
(purity increases and entropy decreases) when the nuclei re-separate as the
interaction with the fusion environment diminishes. Asymptotically, the
elastic and inelastic channels evolve in a coherent superposition.

One useful feature of this method is the ability to follow the dynamics in
time and understand what is happening at each time step of the propagation.
Fig. 2 shows the population of the radial grid basis states over time for an
initial Coulomb wave packet with an average energy $E_{0}=60$ MeV. The change
of population of both the ground state and the $3^{-}$ excited state of 144Sm
due to the radial coupling between these states can be observed.

Fig. 3 displays the dynamics of (i) energy dissipation (thick solid red line),
entropy production (thin solid black line) [panel (a)], and (ii) quantum
coherence (solid line) [panel (b)] for the collision scenario shown in Fig. 2.
The energy loss is determined by the change of the average energy of the
reduced system relative to its initial value, i.e.,
$E_{0}-\textnormal{Tr}[\hat{H}_{S}\rho_{S}(t)]$. It is interesting to observe
that the transient, strong interaction of the reduced system with the fusion
environment makes the dynamics dissipative and decoherent (the entropy
increases and the purity of the reduced-system density matrix decreases). When
the nuclei re-separate, there is a revival of quantum coherence and,
asymptotically, the elastic and inelastic channels move in a coherent
superposition. Does the transient decoherent phase of the collision affect the
fusion probability?

Figure 4 shows the fusion probabilities calculated using optimal conditions
for the initial Gaussian and Coulomb wave packets for the test case collision
16O + 144Sm. The height of the Coulomb barrier between these two nuclei is
$V_{B}=61.1$ MeV. The Gaussian wave packet calculations clearly provide the
best results in terms of energy-resolved fusion probability at deep sub-
barrier energies ($E_{c.m.}/V_{B}<0.9$), reaching an order of magnitude of
$10^{-8}$ (purple triangles). Having accurate values of fusion probability for
low-energy fusion reactions allow us to extend the scope of our calculations
and will allow us to validate interesting effects at these low fusion
probability regions. This could be the case when, for example, using a plasma
environment.

Figure 4:  The energy resolved fusion probabilities from the CCDM method for
different initial wave packets for a head-on collision of 16O + 144Sm with
various mean energies. These probabilities are compared with those from CCFULL
[14]. The CCFULL results (red solid line) are reproduced by the CCDM method
(symbols).

The Coulomb wave packets can still achieve values of $\sim 10^{-6}$ (green
solid squares), however the broadness of these wave packets in momentum space
allows a better global description of energy values below, near and above the
barrier ($E_{c.m.}/V_{B}>0.9$), as also shown in Ref. [26]. This could be
useful as less computational time is needed if the lowest penetration values
are not so important. The physically meaningful values of the CCDM fusion
probability are the invariant values when the parameters of the initial wave
packet are changed [16, 17, 18]. The probabilities that are out of this trend
are determined with a very small, initial weighting factor in the window
function and are numerically inaccurate.

The very good agreement between the CCDM and CCFULL fusion calculations
demonstrates the reliability of the CCDM approach. Since CCFULL does not
include decoherence, these results indicate that the transient quantum
decoherence caused by the fusion environment localised in the nucleus-nucleus
potential pocket does not affect the fusion probability. Quantum decoherence
is a dynamical phenomenon [30] which cannot be described by a static coupled-
channels model such as CCFULL. Dynamical and static coupled-channels
calculations based on the Schrödinger equation provide the same (asymptotic)
observables, as shown in Ref. [26]. The same is expected to happen using a
semi-classical dynamical coupled-channels model [31] at above-barrier incident
energies. In such semi-classical calculations, without complex potentials, the
occupation probability of direct reaction channels oscillates in time at
internuclear radii around the distance of minimal approach (e.g., see Fig. 2
in Ref. [31]). These oscillations are due to the reversible couplings among
the reaction channels which have not yet achieved equilibrium in their
occupation probabilities. The reversible couplings do not affect the coherent
superposition of intrinsic energy states. In the present work, the change of
purity in Fig. 3(b) is caused by the effect of the irreversible coupling
between the fusion environment and the reduced system. Without the fusion
environment, the dynamics of the CCDM model would be Hamiltonian, preserving
energy, entropy and the purity of the reduced density matrix.

##  4 Conclusions

We have successfully shown that the CCDM approach can be used for describing
low-energy nuclear fusion reactions. The window operator has been applied to a
density matrix for energy resolving purposes, which in this case, is crucial
for validating the treatment of fusion using the CCDM method. It is shown that
the CCDM and TISE calculations of fusion are in very good agreement with each
other. Unexpectedly, the transient decoherence within the fusion pocket does
not affect the fusion probability [9]. The CCDM approach has the advantage of
allowing us to study quantities throughout the time propagation of the density
matrix, such as coherence dynamics, energy dissipation and entropy production.
Most importantly, we have built a foundation for the study of environmental
effects on a fusion reaction, via the open quantum systems methodology. The
next step is to introduce a physical, plasma environment which is present in
all stellar fusion reactions. This will be implemented in a similar manner as
the fusion environment, using specifically chosen Lindblad operators.
Achieving accurate fusion probabilities well-below the Coulomb barrier will
allow us to extend our study since plasma interactions are more likely to
affect sub-barrier fusion. Further questions, such as the experimentally
observed hindrance of fusion at sub-barrier energies [32] and dissipative
multi-nucleon transfer reactions at sub-barrier energies [33], may also be
addressed.

Acknowledgement. This work has received funding from the Leverhulme Trust (UK)
under Grant No. RPG-2019-325.

## References

  * Yakovlev [1989] D.G. Yakovlev, D.A. Shalybkov, _Degenerate cores of white dwarfs and envelopes of neutron stars: Thermodynamics and plasma screening in thermonuclear reactions_ , Astrophys. Space Phys. 7 (1989) 311.
  * Harston [1999] M. Harston, J.F. Chemin, _Mechanisms of nuclear excitation in plasmas_ , Phys. Rev. C 59 (1999) 2462.
  * Adriana [2006] A. Pálffy, W. Scheid, Z. Harman, _Theory of nuclear excitation by electron capture for heavy ions_ , Phys. Rev. A 73 (2006) 012715.
  * Gosselin [2009] G. Gosselin, N. Pillet, V. Méot, P. Morel, A.Ya. Dzyublik, _Nuclear transition induced by low-energy unscreened electron inelastic scattering_ , Phys. Rev. C 79 (2009) 014604.
  * breuer [2002] H. Breuer, F. Petruccione, _The Theory of Open Quantum Systems_ , Oxford University Press, 2002.
  * rotter [2015] I. Rotter, J.P. Bird, _A review of progress in the physics of open quantum systems: theory and experiment_ , Rep. Prog. Phys. 78 (2015) 114001.
  * sandulescua [1987] A. Sandulescu, H.Scutaru, _Open quantum systems and the damping of collective modes in deep inelastic collisions_ , Ann. Phys. (NY) 173 (1987) 277.
  * vazgen [2016] V.V. Sargsyan, Z. Kanokov, G.G. Adamian, N.V. Antonenko, _Application of the theory of open quantum systems to nuclear physics problems_ , Phys. Part. Nucl. 47 (2016) 157.
  * Dasgupta [2007] M. Dasgupta, D.J. Hinde, A. Diaz-Torres, B. Bouriquet, C.I. Low, G. Milburn, J.O. Newton, _Beyond the coherent coupled channels description of nuclear fusion_ , Phys. Rev. Lett. 99 (2007) 192701.
  * Lindblad1976 [1976] G. Lindblad, _On the generators of quantum dynamical semigroups_ , Comm. Math. Phys. 48 (1976) 119.
  * pesce [1997] L. Pesce, P. Saalfrank, _“Free” nuclear density propagation in two dimensions: the coupled-channel density matrix method and its application to inelastic molecule-surface scattering_ , Chem. Phys. 219 (1997) 43.
  * alexis2008 [2008] A. Diaz-Torres, D.J. Hinde, M. Dasgupta, G.J. Milburn, J.A. Tostevin, _Dissipative quantum dynamics in low-energy collisions of complex nuclei_ , Phys. Rev. C 78 (2008)064604.
  * alexis2010 [2010] A. Diaz-Torres, _Coupled-channels density-matrix approach to low-energy nuclear collision dynamics: A technique for quantifying quantum decoherence effects on reaction observables_ , Phys. Rev. C 82 (2010) 054617.
  * ccfull [1999] K. Hagino, N. Rowley, A. Kruppa, _A programme for coupled-channel calculations with all order couplings for heavy-ion fusion reactions_ , Comp. Phys. Comm. 123 (1999) 143.
  * schafer [1991] K. Schafer, _The energy analysis of time-dependent numerical wave functions_ , Comp. Phys. Comm. 63 (1991) 427.
  * boselli [2015] M. Boselli, A. Diaz-Torres, _Quantifying low-energy fusion dynamics of weakly bound nuclei from a time-dependent quantum perspective_ , Phys. Rev. C 92 (2015) 044610.
  * alexis [2018] A. Diaz-Torres, M. Wiescher, _Characterizing the astrophysical S factor for 12C + 12C fusion with wave-packet dynamics_, Phys. Rev. C 97 (2018) 055802.
  * terry [2019] T. Vockerodt, A. Diaz-Torres, _Describing heavy-ion fusion with quantum coupled-channels wave-packet dynamics_ , Phys. Rev. C 100 (2019) 034606.
  * kosloff [2019] R. Kosloff, _Quantum thermodynamics and open-systems modeling_ , J. Chem. Phys. 150 (2019) 204105.
  * Balantekin [1998] A.B. Balantekin, N. Takigawa, _Quantum tunneling in nuclear fusion_ , Rev. Mod. Phys. 70 (1998) 77.
  * Hagino [2012] K. Hagino, N. Takigawa, _Subbarrier fusion reactions and many-particle quantum tunneling_ , Prog. Theor. Phys. 128 (2012) 1061.
  * thompson [2001] I.J. Thompson, _Methods of direct reaction theories_ , in _Scattering_ , E.R. Pike and P.C. Sabatier (eds), Academic Press, 2001, pp. 1360-1372.
  * bertlmann [2006] R.A. Bertlmann, W. Grimus, B.C. Hiesmayr, _Open-quantum-system formulation of particle decay_ , Phys. Rev. A 73 (2006) 054101.
  * thompson [2004] I.J. Thompson, A. Diaz-Torres, _Modelling effects of halo breakup on fusion_ , Prog. Theor. Phys. Suppl. 154 (2004) 69.
  * thompson [1986] A.R. Barnett, _COULFG: Coulomb and Bessel functions and their derivatives, for real arguments, by Steed’s method_ , Comp. Phys. Comm. 27 (1982) 147.
  * terry [2021] T. Vockerodt, A. Diaz-Torres, _Calculating the S-matrix of low-energy heavy-ion collisions using quantum coupled-channels wave-packet dynamics_ , Phys. Rev. C 104 (2021) 064601.
  * huisinga [1999] W. Huisinga, L. Pesce, R. Kosloff, P. Saalfrank, _Faber and Newton polynomial integrators for open-system density matrix propagation_ , J. Chem. Phys. 110 (1999) 5538.
  * Saalfrank [1995] C. Scheurer, P. Saalfrank, _Density matrix model for hydrogen transfer in the benzoic acid dimer_ , Chem. Phys. Lett. 245 (1995) 201.
  * lenton [2021] J. Lenton, I. Lee, A. Diaz-Torres, _Quantum dynamics of a nucleon in the Fermi accelerator_ , Ann. Phys. (NY) 434 (2021) 168624.
  * schlosshauer [2007] M. Schlosshauer, _Decoherence and The Quantum-to-Classical Transition_ , Springer Verlag, 2007.
  * Laura [2021] L. Moschini, A. Diaz-Torres, _Tracing the dynamical interplay of low-energy reaction processes of exotic nuclei using a two-center molecular continuum_ , Phys. Lett. B 820 (2021) 136513.
  * hindrance [2021] C.L. Jiang, B.B. Back, K.E. Rehm, K. Hagino, G. Montagnoli, A.M. Stefanini, _Heavy-ion fusion reactions at extreme sub-barrier energies_ , Eur. Phys. J. A 57 (2021) 235.
  * mnt [2020] G.G. Adamian, N.V. Antonenko, A. Diaz-Torres, S. Heinz, _How to extend the chart of nuclides?_ , Eur. Phys. J. A 56 (2020) 47.


# Accelerating innovation with software abstractions for scalable
computational geophysics

Mathias Louboutin*,1, Philipp Witte2, Ali Siahkoohi1, Gabrio Rizzuti3,
Ziyi Yin1, Rafael Orozco1 and Felix J. Herrmann1
1 Georgia Institute of Technology, 2 Microsoft Research, 3 Utrecht University

##  1 Summary

We present the SLIM open-source software framework for computational
geophysics, and more generally, inverse problems based on the wave-equation
(e.g., medical ultrasound). We developed a software environment aimed at
scalable research and development by designing multiple layers of
abstractions. This environment allows the researchers to easily formulate
their problem in an abstract fashion, while still being able to exploit the
latest developments in high-performance computing. We illustrate and
demonstrate the benefits of our software design on many geophysical
applications, including seismic inversion and physics-informed machine
learning for geophysics(e.g., loop unrolled imaging, uncertainty
quantification), all while facilitating the integration of external software.

##  2 Introduction

Software development for exploration geophysics has been traditionally driven
by performance. Though this has resulted in very efficient code, the usability
and portability of these codes has been compromised as a result of the absence
of user-level design. Abstractions and user interfaces at a high level have
been developed in recent decades to facilitate research and development. A
number of such interfaces are available, ranging from programming languages,
e.g., Python [1], Julia [2], and Matlab, to domain-specific languages (DSLs),
including RVL [3], Firedrake [4], and Devito [5, 6]. Additionally, there has
been a community initiative towards open-source codes and reproducibility, led
by Madagascar [7]. These efforts have brought attention to the need for user
abstraction to easily translate mathematics into code without having to
manually refactor thousands of lines of code.

Motivated by this background, we introduce our fully open-source software SLIM
framework based on high-level abstractions and separation of concerns. Our
software design relies on three principles: (1) A high-level abstraction that
represents the mathematical problem at hand; (2) Scalability with vertical
integration of high-performance software and DSLs; (3) Inter-operability
through the adoption of language standards (object oriented, multiple
dispatch, inheritence, …). In the following, we will detail these three
points, highlighting their importance and implementation. We will follow with
some illustrations of the integration with external software to demonstrate
how research and development can be eased and potentially accelerated by these
design principles.

##  3 Software design

The fundamental objective of a scientific software framework is to provide
users with an interface that allows them to express their own scientific
problems. In particular, the interface should provide abstractions to define
the problem as closely as possible to its mathematical formulation. With this
design philosophy, research and development turnaround time can be drastically
reduced with quick prototyping and testing of new ideas and algorithms. Our
software framework is grounded in this idea and led to the development of
legacy software [8, 9] adopted by industry for real-world application. Based
on these ideas and modern programming paradigms, we designed a high-level
framework that encapsulates the mathematical definition of geoscientific
problems at every level. These different levels are handled through the
vertical integration of domain-specific languages (DSLs). At the lowest level,
Devito provides a symbolic DSL for the definition of wave-equation and a just-
in-time compiler to target the available hardware with near-optimal
performance. On top of Devito, we developed JUDI, a linear algebra DSL for
wave-equation based modeling and inversion. Finally, we created a machine
learning framework that allows us to integrate Devito propagators and JUDI
Linear Operators into machine learning (ML) frameworks (PyTorch and Flux)
opening the door to ML-augmented geophysical inverse problems and uncertainty
quantification (UQ). Finally, throughout these different layers, we committed
to follow language standards allowing easy interfacing and integration with
external software and frameworks.

##  4 Wave-equation based inversion

Linear operators are at the core of applied geophysics. Some relevant examples
include data filtering tools such as Fourier-based f-k filters, the Radon
transform, NMO corrections to traditional imaging operator such as Kirchhoff
migration, post-stack migration and wave-equation based inversion where the
discrete wave-equation is a linear operator. One of the first frameworks for
abstract matrix-free linear algebra was sPOT[10], later extended to wave-
equation based inversion and distributed computing for wavefield
reconstruction [11]. While its adoption was limited by its implementation in
Matlab, such user oriented abstraction laid the ground for modern abstracted
frameworks such as JOLI[12] in Julia and pyLops in Python [13].

Similarly, wave-equation based inversion can be trivially formulated as simple
least-squares optimization problem for the non-linear wave-equation
represented as an abstract matrix $\mathbf{A}(\mathbf{m})$ parametrized by the
model parameter $\mathbf{m}$:

$\begin{split}\underset{\mathbf{m}}{\operatorname{minimize}}\hskip
5.69046pt\Phi(\mathbf{m})=\sum_{i=1}^{n_{s}}||\mathbf{P}_{r}\mathbf{A}(\mathbf{m})^{-1}\mathbf{P}_{s}^{\top}\mathbf{q}_{i}-\mathbf{d}_{i}||^{2}_{2},\\\
\end{split}$

Building on modern software solutions JUDI [14, 15] provides a high-level
interface allowing to define a solver for this problem with a few lines of
code. For example, FWI using Gauss-Newton updates at each iteration can be
summarized in as little as five lines (Listing 1).

⬇

1# Gauss-Newton method

2F = judiModeling(model, srcGeom, recGeom)

3J = judiJacobain(F, q)

4for j=1:maxiter

5 d_pred = F*q

6 fhistory_GN[j] = .5f0*norm(d_pred - d_obs)^2

7 # Gauss-Newton update

8 p = lsqr(J, d_pred - d_obs)

9 model0.m .= proj(model0.m .- reshape(p, model0.n))

10end

Listing 1: FWI with Gauss-Newton updates using JUDI. The complete example is
available and reproducible in the JUDI github repository.

In addition to offering a high-level, mathematical interface to the wave-
equation and related linear operators, JUDI builds on multiple layers of
abstractions while still providing computational performance that is
comparable to the state-of-the-art[6]. JUDI integrates Devito as a backend for
the actual wave-equation solves. Devito is a stencil DSL [5, 6] that offers a
symbolic user-level interface to the underlying numerical solver. This
symbolic interface allows the user to easily and mathematically define generic
partial differential equations (PDEs), such as the acoustic, tilted-
transverse-isotropic (TTI), or elastic wave-equations. The underlying software
then provides a code generation framework than can generate highly-performant
code for numerous architectures, (CPU, GPUs, ARM, POWER). Its ease of use and
high-level interface have contributed to its adoption in the oil-and-gas
industry for production and research at scale [16]. We summarize the vertical
integration of the different technologies (compiler, Devito, JUDI, …) in
Figure 1, which highlights the different levels of abstraction. The 2007 BP
TTI dataset, for example, can be setup and processed with RTM within days. The
results are shown in Figure 2.

Figure 1: Schematics of the vertical integration in JUDI that enables at scale
inversion through layers of abstractions. Figure 2: 2D RTM on the 2007 BP TTI
model with a marine acquisition. This RTM was run on a GPU without moving off
to the CPU at any time using randomized trace estimation as an extension of
JUDI [17, 18].

###  4.1 Interoperability

While a proper separation of concerns can lead to highly usable, portable, and
performant software, another important aspect to consider is interoperability,
namely the ability to integrate software from different sources and
organizations. We now describe how the design of our software framework allows
for easy interfacing with external frameworks. One major drawback of
proprietary and low-level software is the limited capability to interface and
interact with external software. Through open-source code development, and
committing to language-specific standards (such as proper usage of types,
hierarchy, inheritance, etc.) in Julia and Python, we enable seamless
interoperability with external software and frameworks. This greatly increases
the ability to perform research, which involves being able to easily test and
combine different ideas. For instance, we combine core Julia packages, our
numerical optimization toolbox SlimOptim.jl [19], a constrained optimization
software SetIntersectionProjection.jl [20], and COFII’s wave-equation
propagators [16], in order to setup a full waveform inversion (FWI) exercise.
We were able to perform FWI, in parallel, on the Marmousi-II model by taking
advantage of the best technology available. Additionally, this demonstrative
example, and in general JUDI (or COFII), can trivially be deployed in the
cloud using once again dedicated software abstractions [16, 21].

⬇

1function objective(F, velocity, d_obs)

2 J = jacobian(F, velocity)

3 d_pred = F * velocity # Forward modeling with COFFI

4 G = J’ * ({d_pred.- d_obs) # Gradient with COFII

5 f = .5f0*norm(d_pred .- d_obs)^2 # Misfit

6 return f, G

7end

8# Projection on TV+bounds, SetIntersectionProjection.jl

9prj = setup_projection(...)

10# Projected Quasi-Newton (l-BFGS) with SlimOptim.jl

11sol = pqn(x->objective(F, x, d_obs), vec(v2), prj);

Listing 2: Projected Quasi-Newton FWI with multiple software packages. The
complete example can be found here

Since different software modules can be pieced together with minimal effort,
it greatly aids the development and testing of new research. This flexibility
allows, for instance, the integration of machine learning into conventional
algorithm pipelines in a plug-and-play fashion. In the following section, we
provide some notable examples of this approach.

##  5 Machine learning

Thanks to recent theoretical and practical advances, deep learning has seen a
wide adoption in various geophysical applications ranging from processing
(e.g., denoising, multiple removal) to inverse problems and seismic
interpretation [e.g., 22, 23, 24]. One of the core challenges in the adoption
of deep learning for seismic applications is the integration of existing codes
(e.g. seismic propagators) into deep-learning frameworks based on automatic
differentiation (AD). To be able to differentiate through networks that
contain both standard PyTorch/Tensorflow layers, as well as third party
functions, such as a forward modeling kernel, the third party code must have
manually implemented gradients and must be properly integrated into the deep
learning framework. Our high-level linear algebra abstraction framework, which
sits on top of our Devito-based propagators, facilitates this integration, as
deep learning frameworks like PyTorch or Flux are already able to
backpropagate through a linear operator (namely by applying its adjoint to the
data residual). This allows us to implement deep neural networks in Flux, in
which we can combine standard deep learning layers with external third party
functions, e.g. convolutional layers with migration/demigration operators.
This capability has enabled various research projects in our group, including
surface-related multiple elimination [25], dispersion attenuation [26], and
ML-augmented imaging and uncertainty quantification [27, 28]. We show in the
following two examples of machine learning for exploration geophysics that
take advantage of these high-level abstractions to interface PDE solvers
(JUDI, Devito) with deep learning frameworks.

###  5.1 Seismic imaging with deep priors and uncertainty quantification

Since we can integrate Devito’s highly optimized wave-equation solvers with
the PyTorch deep learning library, we propose to use deep priors [29] to
regularize seismic imaging, where we reparameterize the seismic image as the
output of an untrained convolutional neural network (CNN). This approach acts
as a regularization in the image space [29, 30, 31], which exploits the
inductive bias of the CNN in representing images without noisy artifacts [29].
We perform Bayesian inference via a gradient-based Markov chain Monte Carlo
(MCMC) algorithm [32], where each iteration requires differentiating the
action of the linearized Born scattering operator on the CNN output with
respect to the CNN’s weights. In order to have access to the automatic
differentiation utilities of PyTorch, we expose Devito’s matrix-free
implementations for the migration operator and its adjoint to PyTorch via
Devito4PyTorch [33]. This is exemplified in Listing 3. Figure 3 summarizes the
results of seismic imaging and uncertainty quantification with this approach,
which are borrowed from Siahkoohi et al. [34].

⬇

1def deep_prior_imaging(d_obs, n, maxiter=1000):

2 # PyTorch Integrated Devito Jacobian.

3 J = ForwardBornLayer(model, geometry)

4 # Initialize deep prior CNN with random input.

5 G = deep_prior_net()

6 z = torch.randn([1, 1, n[1], n[2]])

7 # Setup pSGLD MCMC sampler.

8 optim, samples = pSGLD(G.parameters(), lr=0.001), []

9 for i in range(maxitr):

10 # Compute predicted perturbation model and data.

11 x = G(z)

12 d_pred = J(x)

13 # Compute the negative log-posterior.

14 nlp = n_log_posterior(d_pred, d_obs, G.parameters())

15 # Compute the gradient w.r.t. the CNN weights,

16 nlp.backward()

17 # Sample with pSGLD.

18 optim.step()

19 samples.append(x.detach().numpy())

20 return samples

Listing 3: Seismic imaging and uncertainty quantification with PyTorch and
Devito as a wave-equation solver.

(a)

(b)

(c)

Figure 3: Imaging a 2D subset of the Parihaka [35] dataset with deep priors
[34]. (a) LSRTM without any regularization. (b) The conditional (posterior)
mean estimate, using the deep prior as regularization. (c) Normalized
pointwise standard deviation.

###  5.2 Loop-unrolled seismic imaging

Another set of applications that are enabled by the proper integration of
abstract seismic modeling operators into deep learning frameworks are loop-
unrolled optimization algorithms such as the learned primal-dual
reconstruction [36, 37]. These are special types of networks that follows the
general structure of gradient-based optimization algorithms, in which
conventional gradients are augmented by additional neural network layers.
Using our integration of JUDI operators into Flux via the JUDI4Flux package
[38], we apply the loop-unrolled network architecture from [39] to seismic
imaging. Every iteration of the loop-unrolled algorithm consists of computing
the (conventional) gradient of the LS-RTM objective function $g$, using the
forward and adjoint linearized Born scattering operator, followed by the
application of a shallow CNN (Listing 4). The network takes a single
simultaneous (super-) shot record as the input and predicts an LS-RTM (or
true) image. We train the network in a supervised fashion, in which we
minimize the misfit between the predicted image and the true image (whereas in
reality, one would train with available LS-RTM images). Figure 4 shows the
output of the loop-unrolled gradient descent algorithm, Listing 4, after
training the network for 2000 iterations on a training dataset of 2,000 data-
image pairs. The example is not meant to represent a realistic scenario (as
the true images which were used in the training process are obviously
unknown), but serve as a proof of concept on how to augment physical models
with data-driven approaches. The complete example and additional variations
are available at LoopUnrolledSeismicImaging.

⬇

1function loop_unrolled_lsrtm(d_obs, n; maxiter=10)

2 x = randn(Float32, n[1], n[2], 1, 1)

3 nx, ny, nc, nb = size(x)

4 s = randn(Float32, nx, ny, 5, nb) # memory term

5 for j=1:maxiter

6 g = adjoint(J)*(J*vec(x) - d_obs)

7 g = reshape(Flux.normalise(g), nx, ny, 1, 1)

8 u = cat(x, g, s, dims=3)

9 u = bn1(conv1(u)); u = bn2(conv2(u)); u = bn3(conv3(u));

10 s = relu.(u[:, :, 2:6, :])

11 dx = u[:, :, 1:1, :]

12 x += dx

13 end

14 return vec(x)

15end

Listing 4: Example of a physics-augmented neural network for seismic imaging.
The network consists of 10 iterations of a loop-unrolled gradient descent
algorithm, in which the conventional LS-RTM gradient $g$ is augmented through
convolution layers. The input into the network is the observed seismic data
and the output is the predicted image.  Figure 4: Loop-unrolled LSRTM on a
test 2D slice of the 2D overthrust model. (a) True image, (b) RTM of
simultaneous shot and (c) loop-unrolling of.

###  5.3 End-to-end inversion for geological carbon storage monitoring

Finally, we show that we can easily build a framework for seismic monitoring
of geological carbon storage that integrates PDE solvers, deep learning models
and physical constraints. These monitoring problems rely on three different
types of physics [40]: fluid-flow physics, rock physics and wave physics.
Given time-lapse seismic data collected over multiple years, we jointly invert
for the rocks intrinsic permeability which can be used to recover the CO2
concentration. By leveraging our high-level abstractions and AD, we can
directly differentiate through all three physics solvers which map
permeability to seismic data and minimize a fully coupled data misfit
objective function. At each iteration, given an estimate of the permeability
$K$, we model seismic data where the subsurface velocity is translated from
the solution of the fluid-flow simulation. We can then compute the standard
data misfit and use automatic differentiation to compute the permeability
update. Because each step is abstracted, we can easily swap each physical
solver with a different one, such as replacing the fluid-flow solver by a
trained Fourier Neural Operator (FNO) for computational efficiency [41]. We
show in Figure 5 that we can recover the permeability from seismic
measurements using a pre-trained FNO in place of a fluid-flow solver.

⬇

1opt = RMSprop()

2for j=1:maxiter

3 theta = params(K)

4 grads = gradient(theta) do

5 c = S(K); v = R(c); d_pred = F(v)

6 return 0.5f0 * norm(d_pred-d_obs)^2f0

7 end

8 for p in theta

9 update!(opt, p, grads[p])

10 end

11end

Listing 5: Example of an end-to-end coupled inversion for seismic monitoring
of geological carbon storage. A pre-trained FNO $S$ is used as a surrogate for
fluid-flow simulation. In each iteration, we calculating the seismic data
misfit, compute the gradient with respect to input permeability, $K$, via AD,
and update it according to an optimizer $opt$.  Figure 5: End-to-end inversion
of CO2 concentration from seismic measurements. In order to obtain this end-
to-end result, the three operators representing the fluid-flow modeling, rock
property modeling, wave modeling, and their derivatives are integrated by
combining AD and manually defined gradients.

##  6 Conclusions

Through this paper, we have introduced a software design philosophy aimed at
enabling research and development at scale in geoscience, with the potential
to generalize to any inverse problem such as medical imaging. We demonstrated
through carefully chosen example that high-level abstractions allow to express
complex problem in a clear and representative way without incurring additional
computational costs. Our software framework already allows for a wide range of
application, including industry scale inversion and cutting-edge physics-
informed deep learning for geophysics. We intend to build additional
capabilities to tackle modern computing environment such as Cloud computing,
task dedicated accelerators (i.e TPUs) and non-convex optimization.

##  7 Acknowledgment

This research was carried out with the support of Georgia Research Alliance
and partners of the ML4Seismic Center. We thank J. Washbourne at CVX for the
fruitful discussions on open source software development and interoperability.

## References

  * van Rossum and Drake [2009] Guido van Rossum and Fred L. Drake.  _Python 3 Reference Manual_.  CreateSpace, Scotts Valley, CA, 2009.  ISBN 1441412697.
  * Bezanson et al. [2017] Jeff Bezanson, Alan Edelman, Stefan Karpinski, and Viral B Shah.  Julia: A fresh approach to numerical computing.  _SIAM Review_ , 59(1):65–98, 2017.  doi: 10.1137/141000671.  URL https://epubs.siam.org/doi/10.1137/141000671.
  * Padula et al. [2009] Anthony D. Padula, Shannon D. Scott, and William W. Symes.  A software framework for abstract expression of coordinate-free linear algebra and optimization algorithms.  _ACM Trans. Math. Softw._ , 36(2), apr 2009.  ISSN 0098-3500.  doi: ricevector.  URL https://doi.org/ricevector.
  * Rathgeber et al. [2016] Florian Rathgeber, David A. Ham, Lawrence Mitchell, Michael Lange, Fabio Luporini, Andrew T. T. Mcrae, Gheorghe-Teodor Bercea, Graham R. Markall, and Paul H. J. Kelly.  Firedrake: Automating the finite element method by composing abstractions.  _ACM Trans. Math. Softw._ , 43(3), dec 2016.  ISSN 0098-3500.  doi: 10.1145/2998441.  URL https://doi.org/10.1145/2998441.
  * Louboutin et al. [2019] M. Louboutin, M. Lange, F. Luporini, N. Kukreja, P. A. Witte, F. J. Herrmann, P. Velesko, and G. J. Gorman.  Devito (v3.1.0): an embedded domain-specific language for finite differences and geophysical exploration.  _Geoscientific Model Development_ , 12(3):1165–1187, 2019.  doi: 10.5194/gmd-12-1165-2019.  URL https://www.geosci-model-dev.net/12/1165/2019/.
  * Luporini et al. [2020] Fabio Luporini, Mathias Louboutin, Michael Lange, Navjot Kukreja, Philipp Witte, Jan Hückelheim, Charles Yount, Paul H. J. Kelly, Felix J. Herrmann, and Gerard J. Gorman.  Architecture and performance of devito, a system for automated stencil computation.  _ACM Trans. Math. Softw._ , 46(1), apr 2020.  ISSN 0098-3500.  doi: 10.1145/3374916.  URL https://doi.org/10.1145/3374916.
  * Madagascar [2017] Madagascar, 2017.  URL http://www.ahay.org/wiki/Main_Page.
  * Lin and Herrmann [2015] Tim T.Y. Lin and Felix J. Herrmann.  The student-driven hpc environment at slim.  In _Inaugural Full-Waveform Inversion Workshop_ , 08-09 2015.  (Natal, Brazil).
  * Silva and Herrmann [2019] Curt Da Silva and Felix J. Herrmann.  A unified 2d/3d large scale software environment for nonlinear inverse problems.  _ACM Transactions on Mathematical Software_ , 2019.  URL https://slim.gatech.edu/Publications/Public/Journals/ACMTOMS/2019/dasilva2017uls/dasilva2017uls.html.
  * van den Berg and Friedlander [2009] Ewout van den Berg and Michael P. Friedlander.  Spot: A linear-operator toolbox for matlab.  In _SCAIM_ , University of British Columbia, 2009. SCAIM Seminar, SCAIM Seminar.
  * Kumar [2010] Nameet Kumar.  Parallelizing operations with ease using parallel spot, 2010.  URL https://slim.gatech.edu/Publications/Public/Conferences/SINBAD/2010/Fall/kumar2010SINBADpoe/kumar2010SINBADpoe_pres.pdf.
  * Modzelewski et al. [2022] Henryk Modzelewski, Mathias Louboutin, Ziyi (Francis) Yin, and Rafael Orozco.  slimgroup/joli.jl: v0.7.16, January 2022.  URL https://doi.org/10.5281/zenodo.5933224.
  * Ravasi and Vasconcelos [2020] Matteo Ravasi and Ivan Vasconcelos.  Pylops—a linear-operator python library for scalable algebra and optimization.  _SoftwareX_ , 11:100361, 2020.  ISSN 2352-7110.  doi: https://doi.org/10.1016/j.softx.2019.100361.  URL https://www.sciencedirect.com/science/article/pii/S2352711019301086.
  * Witte et al. [2019a] Philipp A. Witte, Mathias Louboutin, Navjot Kukreja, Fabio Luporini, Michael Lange, Gerard J. Gorman, and Felix J. Herrmann.  A large-scale framework for symbolic implementations of seismic inversion algorithms in julia.  _Geophysics_ , 84(3):F57–F71, 03 2019a.  doi: 10.1190/geo2018-0174.1.  URL https://slim.gatech.edu/Publications/Public/Journals/Geophysics/2019/witte2018alf/witte2018alf.pdf.  (Geophysics).
  * Louboutin et al. [2022a] Mathias Louboutin, Philipp Witte, Ziyi Yin, Henryk Modzelewski, and Carlos da Costa.  slimgroup/judi.jl: v2.6.4, January 2022a.  URL https://doi.org/10.5281/zenodo.5893940.
  * Washbourne et al. [2021] John Washbourne, Sam Kaplan, Miguel Merino, Uwe Albertin, Anusha Sekar, Chris Manuel, Scott Mishra, Matthew Chenette, and Alex Loddoch.  _Chevron optimization framework for imaging and inversion (COFII) – An open source and cloud friendly Julia language framework for seismic modeling and inversion_ , pages 792–796.  2021\.  doi: 10.1190/segam2021-3594362.1.  URL https://library.seg.org/doi/abs/10.1190/segam2021-3594362.1.
  * Louboutin and Herrmann [2021a] Mathias Louboutin and Felix J. Herrmann.  Ultra-low memory seismic inversion with randomized trace estimation.  In _SEG Technical Program Expanded Abstracts_ , pages 787–791, 09 2021a.  doi: 10.1190/segam2021-3584072.1.  URL https://slim.gatech.edu/Publications/Public/Conferences/SEG/2021/louboutin2021SEGulm/louboutinp.html.  (IMAGE, Denver).
  * Louboutin [2021] Mathias Louboutin.  slimgroup/timeprobeseismic.jl, October 2021.  URL https://doi.org/10.5281/zenodo.5605055.
  * Louboutin and Herrmann [2021b] Mathias Louboutin and Felix Herrmann.  slimgroup/slimoptim.jl: v0.1.6, August 2021b.  URL https://doi.org/10.5281/zenodo.5196691.
  * Peters and Herrmann [2019] Bas Peters and Felix J Herrmann.  Algorithms and software for projections onto intersections of convex and non-convex sets with applications to inverse problems.  _arXiv preprint arXiv:1902.09699_ , 2019.
  * Louboutin et al. [2022b] Mathias Louboutin, Ziyi (Francis) Yin, and Felix Herrmann.  slimgroup/judi4cloud.jl: First public release, March 2022b.  URL https://doi.org/10.5281/zenodo.6386831.
  * Zhang et al. [2018] Guoyin Zhang, Zhizhang Wang, and Yangkang Chen.  Deep learning for seismic lithology prediction.  _Geophysical Journal International_ , 215(2):1368–1387, 2018.
  * Wang et al. [2019] Benfeng Wang, Ning Zhang, Wenkai Lu, and Jialin Wang.  Deep-learning-based seismic data interpolation: A preliminary result.  _Geophysics_ , 84(1):V11–V20, 2019.
  * Yang and Ma [2019] Fangshu Yang and Jianwei Ma.  Deep-learning inversion: A next-generation seismic velocity model building method.  _Geophysics_ , 84(4):R583–R599, 2019.
  * Siahkoohi et al. [2019a] Ali Siahkoohi, Dirk J. Verschuur, and Felix J. Herrmann.  Surface-related multiple elimination with deep learning.  In _SEG Technical Program Expanded Abstracts_ , pages 4629–4634, 09 2019a.  doi: 10.1190/segam2019-3216723.1.  URL https://slim.gatech.edu/Publications/Public/Conferences/SEG/2019/siahkoohi2019SEGsrm/siahkoohi2019SEGsrm.html.
  * Siahkoohi et al. [2019b] Ali Siahkoohi, Mathias Louboutin, and Felix J. Herrmann.  The importance of transfer learning in seismic modeling and imaging.  _Geophysics_ , 2019b.  doi: 10.1190/geo2019-0056.1.  URL https://slim.gatech.edu/Publications/Public/Journals/Geophysics/2019/siahkoohi2019itl/siahkoohi2019itl.html.
  * Siahkoohi et al. [2019c] Ali Siahkoohi, Mathias Louboutin, and Felix J. Herrmann.  Neural network augmented wave-equation simulation.  Technical Report TR-CSE-2019-1, Georgia Institute of Technology, 09 2019c.  URL https://slim.gatech.edu/Publications/Public/TechReport/2019/siahkoohi2019TRnna/siahkoohi2019TRnna.html.
  * Siahkoohi et al. [2020] Ali Siahkoohi, Gabrio Rizzuti, and Felix J. Herrmann.  A deep-learning based bayesian approach to seismic imaging and uncertainty quantification.  In _EAGE Annual Conference Proceedings_ , 1 2020.  URL https://slim.gatech.edu/Publications/Public/Conferences/EAGE/2020/siahkoohi2020EAGEdlb/siahkoohi2020EAGEdlb.html.
  * Lempitsky et al. [2018] V. Lempitsky, A. Vedaldi, and D. Ulyanov.  Deep Image Prior.  In _2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 9446–9454, June 2018.  doi: 10.1109/CVPR.2018.00984.
  * Cheng et al. [2019] Zezhou Cheng, Matheus Gadelha, Subhransu Maji, and Daniel Sheldon.  A Bayesian Perspective on the Deep Image Prior.  In _The IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 5443–5451, June 2019.
  * Dittmer et al. [2020] Sören Dittmer, Tobias Kluth, Peter Maass, and Daniel Otero Baguer.  Regularization by architecture: A deep prior approach for inverse problems.  _Journal of Mathematical Imaging and Vision_ , 62(3):456–470, 2020.
  * Welling and Teh [2011] Max Welling and Yee Whye Teh.  Bayesian Learning via Stochastic Gradient Langevin Dynamics.  In _Proceedings of the 28th International Conference on Machine Learning_ , ICML’11, pages 681–688, Madison, WI, USA, 2011. Omnipress.  ISBN 9781450306195.  doi: 10.5555/3104482.3104568.  URL https://dl.acm.org/doi/abs/10.5555/3104482.3104568.
  * Siahkoohi and Louboutin [2021] Ali Siahkoohi and Mathias Louboutin.  slimgroup/devito4pytorch: Initial release, March 2021.  URL https://doi.org/10.5281/zenodo.4610087.
  * Siahkoohi et al. [2021] Ali Siahkoohi, Gabrio Rizzuti, and Felix J Herrmann.  Deep bayesian inference for seismic imaging with tasks.  _arXiv preprint arXiv:2110.04825_ , 2021.
  * WesternGeco. [2012] WesternGeco.  Parihaka 3D PSTM Final Processing Report.  Technical Report New Zealand Petroleum Report 4582, New Zealand Petroleum & Minerals, Wellington, 2012.
  * Andrychowicz et al. [2016] Marcin Andrychowicz, Misha Denil, Sergio Gomez, Matthew W Hoffman, David Pfau, Tom Schaul, Brendan Shillingford, and Nando De Freitas.  Learning to learn by gradient descent by gradient descent.  _Advances in neural information processing systems_ , 29, 2016.
  * Adler and Öktem [2018] Jonas Adler and Ozan Öktem.  Learned primal-dual reconstruction.  _IEEE transactions on medical imaging_ , 37(6):1322–1332, 2018.
  * Witte et al. [2019b] Philipp A. Witte, Mathias Loubouting, and Felix J. Herrmann.  Judi4flux: Seismic modeling for deep learning, December 2019b.  URL https://doi.org/10.5281/zenodo.4301018.
  * Adler and Öktem [2017] Jonas Adler and Ozan Öktem.  Solving ill-posed inverse problems using iterative deep neural networks.  _Inverse Problems_ , 33(12):124007, 2017.
  * Li et al. [2020a] Dongzhuo Li, Kailai Xu, Jerry M Harris, and Eric Darve.  Coupled time-lapse full-waveform inversion for subsurface flow problems using intrusive automatic differentiation.  _Water Resources Research_ , 56(8):e2019WR027032, 2020a.
  * Li et al. [2020b] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar.  Fourier neural operator for parametric partial differential equations, 2020b.