# Supersolid Stacks in Antidipolar Bose-Einstein Condensates

K. Mukherjee  M. Nilsson Tengstrand  T. Arnone Cardinale  S. M. Reimann
Mathematical Physics and NanoLund, LTH, Lund University, Box 118, 22100 Lund,
Sweden

###### Abstract

We theoretically investigate a novel supersolid structure taking the form of
stacked, disk-shaped superfluid droplets connected via a dilute superfluid, in
an antidipolar condensate. A phase diagram is determined for varying the
particle number and scattering length, identifying the regions of a regular
dipolar superfluid, supersolid stacks, and isolated stacked disk-shaped
droplets in an experimentally realizable trapping potential. The collective
Bogoliubov excitation spectrum across the superfluid-supersolid phase
transition is studied, and the transition point is found to be associated with
the breaking of the degeneracy of the two lowest-lying modes. The dynamical
generation of the supersolid stacks is also investigated by ramping down the
scattering length across the phase transition. Moreover, we have studied the
impact of vortex-line penetration on the phase transition. We have found that
the presence of a vortex line causes the supersolid region to move towards
weaker contact interactions. Our detailed numerical simulations highlight that
an antidipolar condensate can create such supersolid stacks within an
experimentally reachable parameter regime.

##  I Introduction

Dipolar Bose-Einstein condensates offer an intriguing platform Baranov (2008);
Lahaye  _et al._ (2009); Baranov  _et al._ (2012); Böttcher  _et al._ (2020);
Chomaz  _et al._ (2022) to study the supersolid state of matter Böttcher  _et
al._ (2019a); Tanzi  _et al._ (2019a); Chomaz  _et al._ (2019); Natale  _et
al._ (2019); Tanzi  _et al._ (2019b); Guo  _et al._ (2019); Hertkorn  _et al._
(2019, 2021a), displaying both diagonal and off-diagonal long-range order
Gross (1957); *Gross1958; Yang (1962); Andreev and Lifshits (1969); Chester
(1970); Leggett (1970); Pomeau and Rica (1994); Boninsegni and Prokof’ev
(2012). In a supersolid, translational symmetry is spontaneously broken,
leading to density modulations while partially maintaining superfluid
properties. The concept was originally introduced and debated over an
extensive period of time in the context of 4He Kim and Chan (2004a);
*Kim2004b; Balibar (2010); Kim and Chan (2012); Boninsegni and Prokof’ev
(2012); Chan  _et al._ (2013), however, only recently ultra-cold atoms have
emerged as a viable alternative setup Henkel  _et al._ (2010); Cinti  _et al._
(2010); Saccani  _et al._ (2011); Léonard  _et al._ (2017a); *Leonard2017b;
Lin  _et al._ (2011); Li  _et al._ (2016, 2017). Experiments involving highly
magnetic dysprosium Kadau  _et al._ (2016); Ferrier-Barbut  _et al._ (2016);
Schmitt  _et al._ (2016); Böttcher  _et al._ (2019b) and erbium atoms Chomaz
_et al._ (2016, 2018) uncovered the emergence of supersolidity in (quasi-)one
Tanzi  _et al._ (2019a); Böttcher  _et al._ (2019b); Chomaz  _et al._ (2018);
Sohmen  _et al._ (2021) as well as two spatial dimensions Schmidt  _et al._
(2021); Biagioni  _et al._ (2022); Bland  _et al._ (2022a). The underlying
mechanism is attributed to the interplay between inter-particle interactions
and quantum fluctuations Lima and Pelster (2011); Wächtler and Santos (2016a);
*Wachtler2016b; Bisset  _et al._ (2016), a mechanism similar to the one seen
in binary Bose gases forming droplets Petrov (2015); Petrov and Astrakharchik
(2016); Jørgensen  _et al._ (2018) which also have been realized
experimentally Cabrera  _et al._ (2017); Semeghini  _et al._ (2018); Skov  _et
al._ (2021). Owing to the unique and fascinating properties of dipolar
droplets and supersolids, there has been a recent surge in research focused on
understanding these systems Bisset and Blakie (2015); Xi and Saito (2016);
Blakie (2016); Saito (2016); Bisset  _et al._ (2016); Macia  _et al._ (2016);
Baillie  _et al._ (2017); Edler  _et al._ (2017); Baillie and Blakie (2018);
Roccuzzo and Ancilotto (2019); Zhang  _et al._ (2019); Blakie  _et al._
(2020); Mishra  _et al._ (2020); Chomaz (2020); Poli  _et al._ (2021);
Hertkorn  _et al._ (2021b, a); Zhang  _et al._ (2021); Young-S. and Adhikari
(2022); Ghosh  _et al._ (2022); Gallemí and Santos (2022); Schmidt  _et al._
(2022); Halder  _et al._ (2022a); Tengstrand  _et al._ (2021); Roccuzzo  _et
al._ (2022), including the exploration of out-of-equilibrium dynamics Sohmen
_et al._ (2021); Ilzhöfer  _et al._ (2021); Norcia  _et al._ (2022); Mukherjee
and Reimann (2022), vortices Roccuzzo  _et al._ (2020); Gallemí  _et al._
(2020); Šindik  _et al._ (2022), and extensions to dipolar mixtures Smith  _et
al._ (2021); Bisset  _et al._ (2021); Li  _et al._ (2022); Scheiermann  _et
al._ (2023); Bland  _et al._ (2022b); Halder  _et al._ (2022b).

The bulk of current and recent work set focus on dipole-dipole interactions
under a fixed magnetic field, oriented along a particular direction, say, the
$z$-axis. However, by rotating the polarizing magnetic field (a technique also
demonstrated in experiments Tang  _et al._ (2018)) it is possible to
manipulate the dipole-dipole interaction Giovanazzi  _et al._ (2002); Baillie
and Blakie (2020). A time-averaged DDI can be utilized when the rotation
frequency exceeds the trap frequencies but is still smaller than the Larmor
frequency. The magnitude and polarity of this DDI is determined by the angle
between the dipole and the $z$-axis Prasad  _et al._ (2019). By adjusting the
tilt angle, the effective interaction can be modulated from dipolar to
antidipolar regimes where the dipolar interaction is reversed. In this case,
head-to-tail antidipoles repel, while side-by-side antidipoles attract, which
is just opposite to the behavior observed in the ordinary dipolar regime. In a
single trapped antidipolar condensate, stacks of disk-like droplets may form
that can even be supersolid when the disks connect by a dilute superfluid.
These structures are quite different from the usual linear arrays of elongated
filaments found in conventional dipolar condensates, as mentioned early-on in
the PhD work of Wenzel Wenzel (2022). Stacked droplets or supersolids however
have, to the best of our knowledge, not yet gained much attention, despite the
in principle realistic experimental scenario, and the many interesting future
prospects concerning for example studies of vorticity in such systems. Similar
stacked supersolid structures were previously only discussed for a different
setting of dipolar and non-dipolar mixtures Kirkby  _et al._ (2023), where the
immiscibility of the components stabilized the system. Other theoretical Ghosh
_et al._ (2022); Bilitewski  _et al._ (2023) or experimental studies Natale
_et al._ (2022); Du  _et al._ (2023) with focus on layered structures of
dipolar gases required an optical lattice potential Natale  _et al._ (2022);
Bilitewski  _et al._ (2023); Du  _et al._ (2023) or an electric dipole moment
Ghosh  _et al._ (2022) for their fabrication.

Here, we set focus on the antidipolar single-component condensates where
remarkably, the rather unique stack structures can arise exclusively from the
combined effects of long-range and contact forces Wenzel (2022). The layers
can be completely isolated from each other, or be connected by a superfluid
link, forming a novel supersolid that we in the following refer to as a
“supersolid stack”. We determine the ground-state phase diagram depicting the
superfluid (SF) phase, supersolid stack (SSS), and isolated stacked droplets
(SDL) as a function of particle number and scattering length, with a fixed
antidipolar length. After demarcating the explicit phase boundaries, we
proceed to analyze the low-lying collective excitation spectra across the
transition using a linear stability analysis based on the Bogoliubov de-Gennes
(BdG) approach Wilson  _et al._ (2009a, b); Lu  _et al._ (2010); Blakie  _et
al._ (2012); Martin and Blakie (2012); Hertkorn  _et al._ (2019). A crucial
finding is that the degeneracy of the two lowest modes is broken at the
transition point, serving as an indicator of the SF-SSS phase transition.
Additionally, we provide evidence that the generation of layered structures
can be achieved dynamically from a superfluid state by performing an
interaction quench. Such a quenching technique for supersolid state generation
is frequently used for regular dipolar condensates Tanzi  _et al._ (2019a);
Böttcher  _et al._ (2019a); Bland  _et al._ (2022a). Quenching the interaction
across the phase boundary results in density oscillations that we have found
to be associated with patterns in the low-lying excitation spectra of the
system. Finally, we investigate the impact of imposing a vortex line through
the stacks. Our analysis reveals that a higher-charge vortex shifts the
supersolid phase towards lower scattering lengths relative to the system
without a vortex. Intriguingly, this suggests that the presence of a vortex
affects the superfluid connection and reduces the density of the crystal
structure.

The remainder of the paper is structured as follows. Section II describes the
extended Gross-Pitaevskii (eGP) equation and the BdG approach. In Sec. III, we
present the ground state properties of an antidipolar condensate, delineating
the phase diagram in Sec. III.1, and describing the collective excitation
properties in Sec. III.2. The dynamical formation of a layered supersolid
structure is addressed in Sec. IV.1, and collective oscillations in Sec. IV.2.
The effect of a vortex line on supersolidity is discussed in Sec. V. We
provide a summary of our findings, along with future perspectives, in Sec. VI.
In Appendix A, we discuss the impact of a three-body loss term. Appendix B
provides some details of our numerical simulations.

##  II MODEL AND METHODS

Figure 1: (Color online) Density isosurfaces representing (a) the regular
superfluid and (b)-(c) the supersolid state formed by layers of four (b) and
three (c) disks for a few selected particle numbers $N$ and scattering lengths
$a$ (see the legends). The density isosurfaces are taken at $50\%$, $16\%$,
and $2.5\%$ of the maximum density. The system is confined in a harmonic
trapping potential with frequencies
$(\omega_{x},\omega_{y},\omega_{z})/(2\pi)=(100,100,50)$, and dipolar length
$a_{\rm dd}=-65.5a_{0}$.

###  II.1 Extended Gross-Pitaevskii Equation

We consider a dipolar BEC (dBEC) composed of atoms with mass $M$ and magnetic
dipole moment $\mu_{m}$. The dBEC is confined within a three-dimensional
harmonic trapping potential of the form
$V(\mathbf{r})=M(\omega_{x}^{2}x^{2}+\omega_{y}^{2}y^{2}+\omega_{z}^{2}z^{2})/2$.
The atomic dipoles are aligned with a rotating uniform magnetic field of
strength $\mathbf{B}(t)$ along $\mathbf{e}(t)$, where
$\mathbf{e}(t)=\mathbf{B}(t)/|\mathbf{B}(t)|$ Prasad  _et al._ (2019); Baillie
and Blakie (2020). If the external magnetic field rapidly rotates at an angle
$\phi$ relative to the $z$-axis and with a much higher frequency than the trap
frequency, it is possible to perform a time-averaging process on the dipole-
dipole interaction (DDI). The time-averaged DDI is given by Giovanazzi  _et
al._ (2002)

$\displaystyle
U_{\mathrm{dd}}(\textbf{r},t)=\frac{\mu_{0}\mu^{2}_{m}}{4\pi}\left(\frac{3\cos^{2}\phi-1}{2}\right)\left(\frac{1-3\cos^{2}\theta}{|\mathbf{r}|^{3}}\right),$
(1)

where $\theta$ denotes the angle between $\mathbf{r}$ and the $z$-axis. The
DDI is highly dependent upon the orientation of the constituent dipoles. The
interaction vanishes entirely at the so-called magic angle, denoted by
$\phi_{m}\approx 54.7^{\circ}$. Specifically, for $\phi<\phi_{m}$, the
preferred orientation of the dipoles is head-to-tail, while for
$\phi>\phi_{m}$, an antidipolar configuration with a side-by-side arrangement
of dipoles becomes energetically favorable. At zero temperature, the system is
well described by the eGP equation Lima and Pelster (2011); Wächtler and
Santos (2016a); *Wachtler2016b; Chomaz  _et al._ (2016)

$\displaystyle
i\hbar\partial_{t}\psi(\textbf{r},t)=\bigg{[}-\frac{\hbar^{2}}{2M}\nabla^{2}+V(\textbf{r})+g\absolutevalue{\psi(\textbf{r},t)}^{2}+\frac{3}{4\pi}g_{\rm
dd}\times$ $\displaystyle\int
d\textbf{r}^{\prime}\frac{1-3\cos^{2}\theta}{\absolutevalue{\textbf{r}-\textbf{r}^{\prime}}^{3}}\absolutevalue{\psi(\textbf{r}^{\prime},t)}^{2}+\gamma(\epsilon_{\mathrm{dd}})\absolutevalue{\psi(\textbf{r},t)}^{3}\bigg{]}\psi(\mathbf{r},t).$
(2)

Here, the short-range repulsive contact interaction, $g=4\pi\hbar^{2}a/M$, is
determined by the $s$-wave scattering length $a$. The dipolar interaction
coefficient is $g_{\rm dd}=4\pi\hbar^{2}a_{\rm dd}/M$ with $a_{\rm
dd}=\mu_{0}\mu^{2}_{m}M(3\cos^{2}\phi-1)/24\pi\hbar^{2}$ being the so-called
dipolar length. The final term in Eq. (II.1) is given by the repulsive Lee-
Huang-Yang (LHY) correction with $\gamma(\epsilon_{\rm
dd})=\frac{32}{3}g\sqrt{\frac{a^{3}}{\pi}}\left(1+\frac{3}{2}\epsilon_{\rm
dd}^{2}\right)$ Lima and Pelster (2011, 2012), where the dimensionless
parameter $\epsilon_{\rm dd}=a_{\mathrm{dd}}/a$ quantifies the relative
strength of the DDI compared to the contact interaction. The solution of Eq.
(II.1) is obtained by employing the split-step Crank-Nicholson method Crank
and Nicolson (1947); Antoine  _et al._ (2013) in imaginary time to determine
the initial ground state, and in real time to monitor the system’s dynamics.
The behavior of the system can be classified into different phases depending
on the absolute value of the parameter $\epsilon_{\mathrm{dd}}$. When
$|\epsilon_{\mathrm{dd}}|$ is sufficiently small, the system exhibits a
superfluid phase. However, for larger values of $|\epsilon_{\mathrm{dd}}|$,
the supersolid phase is favored within a specific range of values of
$|\epsilon_{\mathrm{dd}}|$. Beyond this range, the system transitions into the
isolated droplet phase.

In the following, we consider a dBEC of 164Dy atoms with a magnetic moment of
$\mu_{m}=9.93\mu_{B}$, where $\mu_{B}$ is the Bohr magneton. The system is
examined in the maximally antidipolar regime, with $\phi=\pi/2$. This
particular configuration leads to an overall factor of $-1/2$ in Eq. (1)
compared to non-rotating dipoles corresponding to $\phi=0^{\circ}$.
Consequently, the dipolar length becomes $a_{\rm dd}=-65.5a_{0}$, where
$a_{0}$ is the Bohr radius. The frequencies of the harmonic potential used in
this work are $\omega_{x}/(2\pi)=100$ Hz, $\omega_{y}/(2\pi)=100$ Hz, and
$\omega_{z}/(2\pi)=50$ Hz, resulting in an elongated geometry along the
$z$-axis.

###  II.2 BdG Approach

In order to unveil the collective excitation spectrum of the system, we
perform a Bogoliubov-de Gennes analysis Pethick and Smith (2002); Stringari
and Pitaevskii (2003). To do so, we perturb the equilibrium solution
$\psi_{0}$ using the following ansatz:

$\displaystyle\psi(\mathbf{r},t)=\left\\{\psi_{0}(\mathbf{r})+\epsilon\left[u(\mathbf{r})e^{-i\Omega
t}+v^{\ast}(\mathbf{r})e^{i\Omega t}\right]\right\\}$ $\displaystyle\times
e^{-i\mu t/\hbar}.$ (3)

The parameter $\epsilon$ represents a small-amplitude perturbation, while
$\mu$ is the chemical potential. The eigenfrequencies and eigenfunctions,
represented by $\Omega$ and $(u,v^{*})^{T}$, respectively, are the solutions
of the eigenvalue problem resulting from the substitution of Eqs. (II.2) into
Eqs. (II.1) with terms retained up to the first order in $\epsilon$.
Specifically,

$\begin{pmatrix}H_{\mathrm{s}}-\mu+X&&X\\\
-X&&-(H_{\mathrm{s}}-\mu+X)\end{pmatrix}\begin{pmatrix}u\\\
v\end{pmatrix}=\hbar\Omega\begin{pmatrix}u\\\ v\end{pmatrix},$ (4)

where $H_{s}=-\frac{\hbar^{2}}{2M}\nabla^{2}+V(\textbf{r})$ is the single
particle Hamiltonian, and the operator $X$ is defined by its action on $q=u,v$
according to

$\displaystyle Xq(\mathbf{r})=$ $\displaystyle\int
d\mathbf{r}^{\prime}U_{\mathrm{dd}}(\mathbf{r}-\mathbf{r}^{\prime})\psi^{\ast}_{0}(\mathbf{r}^{\prime})\psi_{0}(\mathbf{r})q(\mathbf{r}^{\prime})$
(5)
$\displaystyle+g|\psi_{0}(\mathbf{r})|^{2}q(\mathbf{r})+\frac{3\gamma}{2}|\psi_{0}(\mathbf{r})|^{3}q(\mathbf{r}).$

Through the variable transformations $f=(u+v)/\sqrt{2}$ and
$g=(u-v)/\sqrt{2}$, Eqs. (4) can be reduced to two equations with only half
the dimensions of the original ones. The resulting reduced equations are then
solved by using standard diagonalization methods.

Figure 2: (Color online) (a) Phases of the antidipolar condensate in terms of
the particle number $N$ and scattering length $a$. The thick blue line
separates the supersolid disk phases (light grey) from the unmodulated
superfluid phase (light blue), which occurs when the contrast $\mathcal{C}=0$.
On the other hand, the thin red line divides the supersolid disks from
isolated disks (light pink), with reference to contrast $\mathcal{C}>0.99$.
The triangular markers indicate a phase transition that manifests via a
density minimum near the center of the trap, while the circular markers denote
the same via a density maximum. The bottom panels show the integrated density
profiles $n_{\rm 1D}$ along the $z$-direction for (b) $N=7\times 10^{4}$ and
(c) $N=10^{5}$. All other parameters are the same as in Fig. 1. Figure 3:
(Color online) Ground state contrast $\mathcal{C}$ (blue circular markers) and
superfluid fraction $f_{s}$ (red triangular markers) as function of the
scattering length $a$ for two different particle number (a) $N=1.3\times
10^{5}$ and (b) $N=10^{5}$. The insets show the integrated density profiles
close to the supersolid and superfluid transitions at scattering length
$a=105.45a_{0}$ for $N=1.3\times 10^{5}$ [(a)], and $a=104.24a_{0}$ and
$a=104.27a_{0}$ for $N=10^{5}$[(b)]. All other parameters are the same as in
Fig. 1.  Figure 4: (Color online) (a) Frequencies $\Omega/(2\pi)$ of the eight
lowest excitation modes as a function of the scattering length $a$. The two
lower panels show the function $f_{z}(z)=\int dxdyf(\mathbf{r})$ for the three
lowest modes in the (b) supersolid phase for $a=104a_{0}$, and in the (c)
superfluid phase for $a=104.7a_{0}$, normalized to the maximal value
$f_{z}^{\mathrm{max}}$. The system consists of $N=10^{5}$ particles, and all
other parameters are the same as in Fig. 1. Figure 5: (Color online) (a)-(b)
Dynamical generation of supersolid states consisting of stacks of four disks
via an interaction quench. The interaction quench has been performed by
linearly reducing the scattering length from $a=110a_{0}$ to $a=105a_{0}$ in
$15/\omega_{z}\approx 47.75\rm ms$. The density isosurfaces, taken at $20\%$
and $1.5\%$ of the maximum densities, are shown at different instants of time
(see the legends). The system is composed of $N=1.3\times 10^{5}$ particles
with dipolar length $a_{\rm dd}=-65.5a_{0}$ and is confined in a potential
with frequencies $(\omega_{x},\omega_{y},\omega_{z})/(2\pi)=(100,100,50)\rm
Hz$.

##  III Ground State properties

Our investigation focuses on exploring the ground-state properties of the
antidipolar condensate, with the primary objective of identifying various
density structures, determining their region of existence in a phase diagram,
and mapping the underlying collective excitation properties during the phase
transition. Furthermore, we explore the impact of the number of disks on the
sharpness of the phase transition.

###  III.1 Identification of Different Phases

Figure 1 displays various ground-state density profiles featuring the
unmodulated and modulated phases. For large scattering lengths where the
contact interaction dominates the dipolar interaction, a typical superfluid
density profile elongated along the $z$-direction is formed, as shown in Fig.
1(a) for $N=1.3\times 10^{5}$ particles at $a=110a_{0}$. As the relative
strength of the dipolar interaction $|\varepsilon_{\mathrm{dd}}|$ increases,
side-by-side configurations become more favored, leading to their arrangement
in the radial plane. However, the presence of a tight radial confinement
breaks the condensate into multiple segments along the $z$-direction as the
configuration with lowest energy, resulting in the formation of stacked disks.
Examples of modulated density structures for two different particle numbers
are presented in Figs. 1(b)-(c), which demonstrate that the individual disks
are connected by a dilute background density, substantiating a SSS phase. Note
that the number and distribution of disks along the $z$-direction strongly
depend on the total number of particles.

Next, we determine the ground state phase diagram, considering different atom
numbers $N$ and scattering lengths $a$ for a fixed $a_{\rm dd}=-65.5a_{0}$. To
differentiate between distinct emergent phases, we introduce a contrast
measure defined as Böttcher  _et al._ (2019a); Bland  _et al._ (2022b)

$\displaystyle\mathcal{C}=\frac{n_{\rm max}-n_{\rm min}}{n_{\rm max}+n_{\rm
min}},$ (6)

where $n_{\rm max}$ and $n_{\rm min}$ are the neighboring maximum and minimum
densities, respectively. A SF phase lacks periodic density modulation,
resulting in $\mathcal{C}=0$, while a completely isolated density modulated
phase corresponds to $\mathcal{C}=1$. Thus, in the SSS phase, the following
relation, $0<\mathcal{C}<1$, holds. The thick blue line in Fig. 2 indicates
the transition from a regular SF to SSS. Notably, the critical scattering
length at which the transition occurs increases with the particle number $N$.
Additionally, the number and positions of disks are sensitive to $N$ (and
hence density). For instance, an even number of disk-shaped droplets can form
along the $z$-direction with a dilute superfluid at $z=0$ (see the triangular
markers). This formation occurs for particle number between $1.03\times
10^{5}<N<1.04\times 10^{5}$. However, for a lower number of particles, one
disk is present at $z=0$, leading to an odd number of droplets (see the
circular markers). Nevertheless, further decreasing $a$ results in the gradual
suppression of the density between individual disks, ultimately rendering them
completely isolated from each other. The onset of the SDL phase is indicated
in Fig. 2 by the thin red line. This transition also manifests by a density
maximum (circular markers) or minimum (triangular markers) at $z=0$, depending
on particle number. In our current setup, we notice that an odd number of
disks with one situated at the center form for particle numbers between
$9\times 10^{4}<N<9.1\times 10^{4}$. Typical examples for isolated three- and
four-disk stacks are demonstrated for $N=7\times 10^{4}$ and $N=10^{5}$
particles, respectively, in Figs. 2(b)-(c), showing integrated density
profiles defined as $n_{\rm 1D}(z)=\int dxdy\absolutevalue{\psi(x,y,z)}^{2}$.
This observable, which we also utilize later, is experimentally detectable,
e.g. via in-situ imaging Hertkorn  _et al._ (2021a); Chomaz  _et al._ (2022).

To gain a deeper understanding of the nature of the transition, we investigate
the contrast $\mathcal{C}$ as a function of the scattering length $a$ for two
different particle numbers; see Fig. 3. In addition, we use Leggett’s upper
bound to estimate the superfluid fraction $f_{s}$ Leggett (1970, 1998) in the
central region of length $2L$111The length $2L$ spans the central region
between the first two side-minima for three disks and the two first-side
maxima for four disks, along the $z$-axis Scheiermann  _et al._ (2023)..

$\displaystyle f_{s}=\frac{(2L)^{2}}{N}\left[\int dz\frac{1}{\int
dxdy\absolutevalue{\psi}^{2}}\right]^{-1}.$ (7)

For $N=1.3\times 10^{5}$ particles, we observe the continuous formation of
disks on both sides of $z=0$ as $\mathcal{C}(f_{s})$ gradually approaches zero
(one), indicating the formation of modulated density structures with a strong
superfluid connection between them. As the scattering length $a$ decreases, we
observe that $\mathcal{C}$($f_{s}$) increases (decreases) due to a pronounced
density modulation. At $a=102.1a_{0}$, the disks become completely isolated
from each other and $\mathcal{C}>0.99$. Interestingly, for $N=10^{5}$
particles, the superfluid to supersolid transition is much sharper, but it
still remains continuous. We observe a disk formation at $z=0$ for this
transition, unlike in the case of $N=1.3\times 10^{5}$ particles. The contrast
drops from $\mathcal{C}=0.252$ at $a=104.24a_{0}$ to $\mathcal{C}=0$ at
$a=104.27a_{0}$. Our results are consistent with the findings of Ref. Zhang
_et al._ (2019), where the supersolid phase transition is argued to be
continuous at larger particle numbers closer to the thermodynamic limit due to
insignificant kinetic energy contributions. A close inspection also reveals
that the superfluid state just before the transition to the supersolid state
does not exhibit the typical parabolic density. Nonetheless, they are still
classified as superfluid since they demonstrate a non-periodic density
modulation. For instance, for $N=1.3\times 10^{5}$, $n_{\rm 1D}$ has a
tendency to develop flat top [Fig. 3] profile, while for $N=10^{5}$, a sharp
peak appears [Fig. 3(b)]. Turning to the SDL side of the curves, we observe
that the contrast $C$ gradually approaches zero as the system enters the
isolated droplet regime for both $N=1.3\times 10^{5}$ and $N=10^{5}$
particles. Notably, our numerical simulations show that the superfluid
fraction, as calculated by Eq. (7), does not completely vanish during the SSS-
SDL transition. Additionally, $f_{s}$ gradually decreases and becomes
$f_{s}<0.1$ when $\mathcal{C}>0.99$, indicating the transition to the SDL
phase. This observation is consistent with the findings of Ref. Kirkby  _et
al._ (2023).

###  III.2 Collective Excitation Spectra

Having identified the ground-state phase diagram and showcased some density
profiles, we turn to study the behavior of the collective excitations. The
low-lying collective excitation spectra for the supersolid phase in regular
dipolar condensate have been reported in Refs. Chomaz  _et al._ (2018); Natale
_et al._ (2019); Hertkorn  _et al._ (2019); Schmidt  _et al._ (2021); Hertkorn
_et al._ (2021a, c); Bühler  _et al._ (2022). To accomplish this for the
antidipolar condensate, we fix the particle number at $N=10^{5}$ and compute
the spectrum over a range of values of the $s$-wave scattering length by
solving the BdG Eqs. (II.2). The obtained results are displayed in Fig. 4.
Specifically, we have examined the low-lying excitations that are primarily
related to excitations in the $z$-direction. The frequencies of the eight
lowest modes are presented in Fig. 4(a). We notice that as the scattering
length is decreased from $a=105a_{0}$, the two lowest (degenerate) modes
eventually split into two around $a\sim 104.26a_{0}$, marking the transition
from the superfluid to the supersolid phase. It is interesting to note that
the point at which the degeneracy breaks coincides with the point where the
contrast drops to zero, as shown in Fig.3(b). This implies that the breaking
of degeneracy can be used as a precise definition of the transition point from
the superfluid to the supersolid phase. This behavior is analogous to what has
been observed for a regular dipolar supersolid in an elongated harmonic trap
Hertkorn  _et al._ (2019). To investigate the nature of these modes further,
we note that to the lowest non-vanishing order in $\varepsilon$, the density
can be written as

$n(\mathbf{r},t)=\psi_{0}^{2}(\mathbf{r})+2\sqrt{2}\varepsilon
f(\mathbf{r})\psi_{0}(\mathbf{r})\cos(\Omega t).$ (8)

The time-evolution of the density associated with a mode is consequently
characterized by the function $f$, and we display the integrated version
$f_{z}(z)=\int dxdyf(\mathbf{r})$ of the three lowest modes in Figs. 4(b)-(c).
The modes in the SSS phase at $a=104a_{0}$ are shown in Fig. 4(b), while the
modes in the superfluid phase at $a=104.7a_{0}$ are shown in Fig. 4(c).
Interestingly, the two lowest modes are very similar in both phases, although
their effects are quite different due to their different spatial
distributions. For the supersolid, the lowest mode corresponds to a center of
mass excitation where the high-density regions move back and forth. Thus it is
a Nambu-Goldstone mode associated with the spontaneous breaking of translation
symmetry Nambu (2009). The second-lowest mode is a Higgs (amplitude) mode
Pekker and Varma (2015), where the crystalline regions periodically increase
and decrease in density in a fashion opposite to that of the superfluid
background. The third-lowest mode is the dipole mode, which has a constant
frequency $\Omega/(2\pi)=50$ Hz, independent of the interaction strength Kohn
(1961). The dipole mode remains decoupled from other modes, thereby serving as
a gauge for numerical accuracy. In the superfluid phase [Figs. 4(c)], the
first two modes, $\Omega_{1}$ and $\Omega_{2}$, have the same frequency, and
the third mode is the dipole mode. This dipole mode changes character
significantly as the system transitions from a superfluid to a supersolid. The
function $f_{z}$ of the dipole mode in the supersolid phase has a very
prominent dipole-like structure that is composed of a single positive and
negative region, see the red dashed curve in Fig. 4(c). This mode in the
supersolid phase [red dashed line, $\Omega_{3}$, in the Fig. 4(b)] acquires a
behaviour that appears to be a mix of its lowest mode [solid black line,
$\Omega_{1}$, in Fig. 4(b)] and the dipole mode in the superfluid phase
[$\Omega_{3}$ in Fig. 4(c)], possessing multiple local minima and maxima. This
mode is often attributed to in-phase (between the crystal and the superfluid
background) centre-of-mass oscillation, referred to as the in-phase Goldstone
mode Guo  _et al._ (2019).

##  IV Quench Dynamics across the phase transition

After explicating the supersolid state and its underlying collective
excitation spectrum within the ground-state phase diagram, we now shift our
focus to studying the non-equilibrium dynamics of our system. To trigger these
dynamics, we will vary the scattering length from an initial value of
$a=a_{i}$ to a final value of $a=a_{f}$. A particular focus will be on
examining how the layers of disks can be generated dynamically. We will also
illustrate the emergence of excitations that manifests as density oscillations
and how the excitation frequencies relate to those calculated via the BdG
analysis at $a=a_{f}$.

###  IV.1 Dynamical Generation of Supersolid Disks

We present our results for a system comprising $N=1.3\times 10^{5}$ particles
by showing the 3D density isosurfaces in Figs. 5(a)-(f). Our investigation
begins with an initial state prepared at $a_{i}=110a_{0}$. We then introduce a
very small amplitude noise [Figs. 5(a)] to the ground state, followed by a
gradual decrease of the scattering length to $a_{f}=105a_{0}$ over a time span
of $47.75$ ms222This time span is chosen such that we can adiabatically
produce the states which have the same number of droplets as those in the
ground state of the system. A sudden quench would typically produce more
droplets than the ground state Böttcher  _et al._ (2019a).. As the system
dynamically enters the SSS phase within the ground-state phase diagram [Fig.
2(a)], a modulational instability sets in, resulting in the onset of density
modulation. This can be observed by careful inspection in Fig. 5(b).
Subsequently, for $t>47.75$ ms, the system quickly forms four prominent
circular disks arranged in layers, connected by a dilute superfluid background
as depicted in Fig. 5(c). A more precise analysis of the temporally resolved
dynamics leading to the SSS phase can be done by invoking the integrated
density $n_{\rm 1D}(z,t)$. The time evolution of $n_{\rm 1D}$ for two
different particle numbers is shown in Fig. 6. Indeed, the formation of four
disks, with the central two having the highest densities from an initial non-
modulated state, is evident in Fig. 6(a), corresponding to $N=1.3\times
10^{5}$ particles. Similarly, an odd number of disks can be formed
dynamically. An example of the formation of a state with three disks is
illustrated in Fig. 6(b) for $N=7\times 10^{4}$. The final scattering length
$a_{f}=102a_{0}$ is achieved from an initial $a_{i}=110a_{0}$ in the same time
span of $47.75$ ms as before. The superfluid connection remains more robust
for $N=1.3\times 10^{5}$ particles when compared to $N=7\times 10^{4}$. We
remark that the value of contrast at the ground state corresponds to
$\mathcal{C}=0.356$ at $a=105a_{0}$ for $N=1.3\times 10^{5}$, and
$\mathcal{C}=0.54$ at $a=102a_{0}$ for $N=7\times 10^{4}$, implying a weaker
background for the latter, which also becomes evident in dynamically generated
SSS phase [see Figs. 6(a)-(b)]. Another important observation from Figs.
5(d)-(f) is that the thickness of each individual crystal changes during the
dynamics, implying particle flow among them as well as the emergence of
density oscillations. This particle flow causes a persistent alteration of the
density, radial width, and thickness [see Figs. 6(a)-(f)], thereby triggering
compressional dynamics, which we will expound on next.

Figure 6: (Color online) Generation of supersolid states for two different
particle numbers (a) $N=1.3\times 10^{5}$, and (b) $N=7\times 10^{4}$,
utilizing the dynamic ramp of the scattering length $a$ from an initial value
$a=110a_{0}$ to the final value $a=105a_{0}$[(a)], and $a=102a_{0}$[(b)]. The
colorbar represents the integrated density $n_{\rm 1D}$ in unit of $1000$
$\mu\mathrm{m}^{-1}$. The duration of dynamic ramp is $15/\omega_{z}\approx
47.75$ ms. All other parameters are the same as in Fig. 1.  Figure 7: (Color
online) (a) Time evolution of the integrated density difference with respect
to the ground state $\psi_{0}$ at the final scattering lengths is quantified
via $\Delta n(t)=(1/N)\int
d\mathbf{r}(\absolutevalue{\psi(\mathbf{r};t)}^{2}-\absolutevalue{\psi_{0}{(\mathbf{r})}}^{2})$
for two different particle numbers (see the legends). The dynamics is
triggered by ramping the scattering length from $a=110a_{0}$ to $a=105a_{0}$
($N=1.3\times 10^{5}$), and $a=102a_{0}(N=7\times 10^{4})$. The real and
imaginary parts of the collective excitation spectra $\Omega$ are presented in
(b) for $N=1.3\times 10^{5}$ and (c) for $N=7\times 10^{4}$, calculated by BdG
analysis at $a=105a_{0}$ and $a=102a_{0}$, respectively. The BdG frequencies
are shown in units of $2\pi\times$Hz.

###  IV.2 Collective Oscillation During Dynamics

Here we evaluate the extent to which the density of the system during the time
evolution remains similar to that of the ground state at the same scattering
length. This can be monitored by analyzing the observable $\Delta n(t)$, as
given by

$\displaystyle\Delta n(t)=(1/N)\int
d\mathbf{r}(\absolutevalue{\psi(\mathbf{r},t)}^{2}-\absolutevalue{\psi_{0}{(\mathbf{r})}}^{2}),$
(9)

where $\psi_{0}$ represents the ground state wavefunction at $a_{f}$. When the
value of $\Delta n(t)$ approaches zero, it indicates that the dynamically
generated supersolid state is very similar to the ground state. However, any
deviation from zero arises due to the collective excitations triggered during
the dynamics. The behavior of $\Delta n(t)$, which relates to the density
evolution displayed in Fig. 6(a)-(b), has been presented in Fig. 7(a). Note
that $\Delta n(t=0)$ is higher for $N=1.5\times 10^{5}$ than for $N=7\times
10^{4}$ due to a stronger superfluid connection in the final state,
$\psi_{0}$, of the former, making it much closer to the initial superfluid
ground state. As the scattering length is gradually ramped down to the final
value, $\Delta n(t)$ gradually decreases and eventually reaches a minimum, as
shown in Fig. 7(a). This minimum value indicates the generation of a density-
modulated state that closely resembles the ground state at the final
scattering length $a_{f}$. Subsequently, $\Delta n(t)$ exhibits an oscillatory
behaviour that involves multiple oscillation frequencies. The oscillation
amplitude is larger for $N=7\times 10^{4}$ than $N=1.3\times 10^{5}$ because
the change in scattering length $a_{f}-a_{i}$ is larger in the same time span,
making the quench less adiabatic in the former case. As mentioned earlier,
such oscillations in $\Delta n(t)$ stem from the continuous particle flow
between the crystal and background, leading to compressional dynamics. The two
dominant frequencies of oscillations are calculated by performing a Fourier
transformation of $\Delta n(t)$333To compute the frequencies, we have
propagated $\Delta n(t)$ until the time $T_{f}=1$ s, and as a result, the
maximum frequency error is estimated to be $1/T_{f}=1$ Hz., yielding
$\omega_{f_{1}}/(2\pi)\approx\rm 36\pm 0.5Hz$ and
$\omega_{f_{2}}/(2\pi)\approx\rm 72\pm 0.5Hz$ for $N=1.3\times 10^{5}$, and
$\omega_{f_{1}}/(2\pi)\approx\rm 28\pm 0.5Hz$ and
$\omega_{f_{2}}/(2\pi)\approx\rm 72\pm 0.5Hz$ for $N=7\times 10^{4}$. Next, we
check whether these frequencies are related to the low-lying excitations of
the ground state at the final scattering length, $a_{f}$. It should be noted
that the excitation spectra at scattering lengths other than the final value
may also be triggered during the quenching process. However, their
contribution to the oscillation of $\Delta n(t)$ is expected to be negligible,
as the system spends most of the time at $a_{f}$. The frequencies of the low-
lying collective excitation spectra calculated at $a=105a_{0}$ and
$a=102a_{0}$ are presented in Figs. 7(b) and (c) for $N=1.3\times 10^{5}$ and
$N=7\times 10^{4}$, respectively. Notably, the dominant oscillation frequency
in $\Delta n(t)$ is very close to the second frequency (excluding the zero
frequency) that corresponds to the Higgs amplitude mode, as depicted in Figs.
7(b)-(c). Additionally, the second dominant frequency close to $72\rm Hz$ is
also evident in the spectrum. Finally, it should be noted that we have
illustrated the intricate connection discussed above by considering these two
particle numbers as representative of three- and four-crystal states. However,
the discussion is generally valid for all states in the SSS phase, where the
non-equilibrium dynamics induced by the quench indeed emanates from triggering
the low-lying excitations of the system.

Figure 8: (Color online) Density isosurfaces within a supersolid state in the
presence of a vortex line with charge (a) $S=1$ and (b) $S=2$. (c) The
contrast $\mathcal{C}$ as a function of scattering length $a$ in the presence
($S=1,S=2$) or absence ($S=0$) of a vortex line. The density isosurfaces
represent $20\%$ and $2.5\%$ of the maximum density. The 164Dy system is
composed of $N=10^{5}$ particles with dipolar length $a_{\rm dd}=-65.5a_{0}$,
and confined in a traping potential with frequencies
$(\omega_{x},\omega_{y},\omega_{z})/(2\pi)=(100,100,50)$ Hz.

##  V Impact of a vortex line

Let us finally briefly illustrate the effect of a non-linear defect on the
formation of multiple-disk structures and the associated superfluid
background. As the disks are relatively flat structures, they offer an
opportunity to investigate quasi-2D physics, particularly in relation to
vortices. Hence, it is crucial to examine the influence of vortex line
penetration on the phase transition. To realize a vortex line of charge $S$
passing through the center of the disks, we perform the following
transformation

$\displaystyle\psi(x,y,z)=\psi(x,y,z)e^{iS\tan^{-1}(y/x)}$ (10)

to the wavefunction during the imaginary time evolution, which ensures a
vanishing density at the center.

Figures 8(a) and 8(b) depict the 3D isosurfaces representing the supersolid
state for $S=1$ and $S=2$, respectively, realized at $a=103a_{0}$ for
$N=10^{5}$ particles. The vortex line has a larger radial thickness in the
dilute superfluid background than in the localized crystal structure, which
varies between $S=1$ and $S=2$. This suggests that the superfluid connection
between the crystal changes due to the presence of a vortex. To investigate
this further, we explore whether imprinting a vortex line changes the
supersolid region, as shown in Fig. 8. The contrast $\mathcal{C}$ is plotted
as a function of the scattering length $a$ in both the presence and absence of
a vortex line. The results show that a higher-charge vortex shifts the
supersolid phase towards a lower scattering length compared to the vortex-free
system. However, the contrast asymptotically approaches each other towards the
isolated droplet phases, indicating that the SSS to SDL transition is not
affected by the presence of a vortex. Between the SF and SSS phase, the
contrast decreases sharply towards lower values for higher vortex charges. It
is energetically costly for a vortex to dig a hole within the highly localized
crystal compared to the dilute superfluid, and therefore it increases
superfluid connection, minimizing the crystal density.

In addition, we remark that we have analyzed the stability of vortex lines in
real-time dynamics. It is observed that a single charge vortex remains stable,
but the double-unit charged vortex line eventually bends and breaks. However,
the detailed discussion of their dynamics is beyond the scope of this
manuscript.

##  VI Conclusions

In this paper, we have reported a novel supersolid state formed in the form of
stacked disk-like droplets connected by a dilute superfluid in an antidipolar
condensate Wenzel (2022). Considering an elongated geometry with
experimentally relevant trapping frequencies we have presented a phase diagram
in the parameter space of particle number $N$ and scattering length $a$,
identifying the regions of existence for a superfluid, supersolid disks, and
isolated disks.

We have utilized a contrast measure and superfluid fraction to differentiate
between different emergent phases. As the scattering length decreases, the
contrast (superfluid fraction) increases (decreases) due to a pronounced
density modulation, leading to the formation of disk-shaped density
modulations. The critical scattering length for the transition from a regular
superfluid to a supersolid stack phase increases with the particle number. The
number and distribution of disks along the $z$-direction are strongly
dependent on the total number of particles.

Building upon the phase diagram, we have studied the collective excitation
spectrum across the superfluid-supersolid phase transition in the antidipolar
condensate. Specifically, we have computed the frequency of the eight lowest
modes and characterized the density of the first three modes. The breaking of
degeneracy of the two lowest-lying modes marks the transition point from the
superfluid to the supersolid stack phase. Additionally, the third lowest mode,
the so-called dipole mode, also changes significantly as the system undergoes
a phase transition. We have also examined the influence of vortex line
penetration on the phase transition, with the supersolid region shifting
towards weaker contact interaction upon imprinting a vortex line.

As a next step, we have studied the non-equilibrium dynamics triggered by
ramping down the scattering length across the superfluid and supersolid phase
transition. We have demonstrated the dynamic generation of four and three
disks connected by a dilute superfluid, considering two different particle
numbers. The particle flow between the disks triggers a collective oscillation
in the density. We have identified two dominant frequencies of oscillations,
which can also be found in the low-lying collective excitation spectra of the
ground state at the final scattering length.

The present work has opened up numerous promising research directions for
future endeavors. To further enhance our understanding, it would be
fascinating to investigate the response of the disks and background superfluid
to external rotation and dynamically observe how vortex lines penetrate
through the disk-like density modulations. Currently, vortices are being
sought actively in dipolar supersolids, and the recent creation of vortices in
an unmodulated dipolar condensate presents optimistic prospects in this area
Klaus  _et al._ (2022); Bland  _et al._ (2023). The 2D nature of the crystal
makes the antidipolar condensate an intriguing setup along that direction.
Further exploration will be undertaken to examine supersolid stacks and
isolated stacked droplets in the context of binary mixtures Arazo  _et al._
(2023). Direct formation of such a supersolid structure via evaporative
cooling along the line of Ref. Bland  _et al._ (2022a) and investigating the
thermal properties Sánchez-Baena  _et al._ (2023) would be highly interesting.

## Acknowledgements

This work was financially supported by the Knut and Alice Wallenberg
Foundation and the Swedish Research Council. Many valuable inputs from P.
Stürmer and L. Chergui are gratefully acknowledged.

##  Appendix A Three-Body Atom Loss

Figure 9: (Color online) (a) The time evolution of the integrated density
$n_{\rm 1D}$ for a quench from an unmodulated state to a supersolid state,
including three-body atom loss with a three-body coefficient $K_{3}=1.5\times
10^{-40}\text{~{}}\rm m^{6}/s$. The quench has been performed by ramping off
the scattering length from $a=110a_{0}$ to $a=102a_{0}$ in $48$ ms. (b) The
time evolution of the atom number $N$ for different loss coefficients (see the
legends), considering the same quench protocol. The system consists of
$N=70000$ 164Dy atoms confined in a harmonic trapping potential with
frequencies
$(\omega_{x},\omega_{y},\omega_{z})/(2\pi)=(100,100,50)\text{~{}}\rm Hz$, with
$a_{\rm dd}=-65.5a_{0}$. The colorbar represents the density in units of
$1000\text{~{}}\mu\mathrm{m}^{-1}$

In experiments, the lifetime of the supersolid state is limited by three-body
losses, which has been previously discussed in the context of dBECs with
positive dipolar lengths Tanzi  _et al._ (2019a); Chomaz  _et al._ (2019);
Böttcher  _et al._ (2019a). This type of losses is caused by highly localized
density arrays, making it crucial to understand how it affects the formation
and persistence of the SSS in anti-dBECs. A comprehensive understanding of
this phenomenon is necessary for the successful experimental realization of
SSS states. In simulations, the loss can be modeled by adding an imaginary
contribution, represented by $(-i\hbar
K_{3}/2)\absolutevalue{\psi(\textbf{r})}^{4}$, to Eq. (II.1), where $K_{3}$
denotes the three-body recombination rate.

To demonstrate our findings, we examine the formation of a supersolid state in
a system composed of $N=70000$ 164Dy atoms by reducing the scattering length
from $a=110a_{0}$ to $a=102a_{0}$ over a period of $47.78$ ms. Realistic
three-body losses are considered in simulations using the experimentally
measured loss coefficient, $K_{3}=1.5\times 10^{-40}\ \rm m^{6}/s$ Böttcher
_et al._ (2019a). Fig. 9(a) shows the time evolution of the integrated density
$n_{\rm 1D}(z)$. The maximum atom loss is $\approx 5.6\%$ of the initial
population in $250$ ms, as displayed by the blue line in Fig. 9(b). After
$t>47.75$ ms, the initially unmodulated state deforms into three disks.
Despite continuous population loss, see Fig. 9(b), a dilute background density
always connects the disks throughout the dynamics. We have also investigated
the population loss and supersolid state generation with larger three-body
loss coefficients, $K_{3}=7.5\times 10^{-40}\text{~{}}\rm m^{6}/s$ and
$K_{3}=1.5\times 10^{-39}\text{~{}}\rm m^{6}/s$, respectively. As expected,
larger $K_{3}$ leads to more atom losses, and $N(t)$ exhibits oscillatory
behavior [Fig. 9(b)]. This tendency is caused by the atom losses being
influenced by the severe changes in density distribution (via the term
$\absolutevalue{\psi(\mathbf{r})}^{4}$), which also undergoes periodic
expansion and contraction after the quench. Nonetheless, the supersolid state
remains robust in long-time dynamics even for $K_{3}=1.5\times
10^{-39}\text{~{}}\rm m^{6}/s$, as evidenced, for example, by the density
isosurface at time $t=240\text{~{}}\rm ms$.

##  Appendix B Computational Details

Here, we provide a detailed account of the numerical simulations used to
obtain the results described in the main text. The extended Gross-Pitaevskii
equation, Eq. (II.1), is cast into a dimensionless form in our simulations by
rescaling the length, the time in terms of the harmonic oscillator length
scale $l_{\rm osc}=\sqrt{\hbar/m\omega_{z}}=1.11\text{~{}}\rm\mu m$, and the
trap frequency $\omega_{z}$, respectively. Then the wavefunction is scaled
accordingly as $\psi(\mathbf{r}^{\prime},t)=\sqrt{l^{3}_{\rm
osc}/N}\psi(\mathbf{r},t)$. We solve the resulting dimensionless equation
using the split-step Crank-Nicholson method, see Ref.Crank and Nicolson
(1947). The stationary (lowest energy) states of the dBEC are obtained through
imaginary time propagation. Furthermore, we apply the transformation
$\frac{\psi(\mathbf{r},t)}{\norm{\psi(\mathbf{r},t)}}\rightarrow 1$ at each
imaginary time-step of this procedure. This preserves the normalization of the
wavefunction, while convergence is reached until relative deviations of the
wave function (at every grid point) and energy between consecutive time-steps
are smaller than $10^{-5}$ and $10^{-7}$, respectively. This solution is then
used as an initial state for the dynamical simulations, where the eGPE is
propagated in real time. It should be noted that calculating the stationary
state solution of Eq. (II.1) is an involved task due to many close-lying local
minima in the energy surface, which necessitates extensive sampling over many
different initial conditions to identify the most probable lowest-energy
solutions. Our simulations are carried out in a 3D box characterized by a grid
$n_{x}\times n_{y}\times n_{z}$ corresponding to $(128\times 128\times 256)$.
The employed spatial discretization (grid spacing) refers to
$\Delta_{x}=0.12l_{\rm osc}$, $\Delta_{y}=0.12l_{\rm osc}$, and
$\Delta_{z}=0.08l_{\rm osc}$, while the time step of the numerical integration
is $\Delta_{t}=10^{-5}/\omega_{z}$.

## References

  * Baranov (2008) M. Baranov, Physics Reports 464, 71 (2008).
  * Lahaye  _et al._ (2009) T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, and T. Pfau, Reports on Progress in Physics 72, 126401 (2009).
  * Baranov  _et al._ (2012) M. A. Baranov, M. Dalmonte, G. Pupillo, and P. Zoller, Chem. Rev. 112, 5012 (2012).
  * Böttcher  _et al._ (2020) F. Böttcher, J.-N. Schmidt, J. Hertkorn, K. S. H. Ng, S. D. Graham, M. Guo, T. Langen, and T. Pfau, Reports on Progress in Physics 84, 012403 (2020).
  * Chomaz  _et al._ (2022) L. Chomaz, I. Ferrier-Barbut, F. Ferlaino, B. Laburthe-Tolra, B. L. Lev, and T. Pfau, Reports on Progress in Physics 86, 026401 (2022).
  * Böttcher  _et al._ (2019a) F. Böttcher, J.-N. Schmidt, M. Wenzel, J. Hertkorn, M. Guo, T. Langen, and T. Pfau, Phys. Rev. X 9, 011051 (2019a).
  * Tanzi  _et al._ (2019a) L. Tanzi, E. Lucioni, F. Famà, J. Catani, A. Fioretti, C. Gabbanini, R. N. Bisset, L. Santos, and G. Modugno, Phys. Rev. Lett. 122, 130405 (2019a).
  * Chomaz  _et al._ (2019) L. Chomaz, D. Petter, P. Ilzhöfer, G. Natale, A. Trautmann, C. Politi, G. Durastante, R. M. W. van Bijnen, A. Patscheider, M. Sohmen, M. J. Mark, and F. Ferlaino, Phys. Rev. X 9, 021012 (2019).
  * Natale  _et al._ (2019) G. Natale, R. M. W. van Bijnen, A. Patscheider, D. Petter, M. J. Mark, L. Chomaz, and F. Ferlaino, Phys. Rev. Lett. 123, 050402 (2019).
  * Tanzi  _et al._ (2019b) L. Tanzi, S. M. Roccuzzo, E. Lucioni, F. Fama, A. Fioretti, C. Gabbanini, G. Modugno, A. Recati, and S. Stringari, Nature 574, 382–385 (2019b).
  * Guo  _et al._ (2019) M. Guo, F. Böttcher, J. Hertkorn, J.-N. Schmidt, M. Wenzel, H. P. Büchler, T. Langen, and T. Pfau, Nature 574, 386–389 (2019).
  * Hertkorn  _et al._ (2019) J. Hertkorn, F. Böttcher, M. Guo, J. N. Schmidt, T. Langen, H. P. Büchler, and T. Pfau, Phys. Rev. Lett. 123, 193002 (2019).
  * Hertkorn  _et al._ (2021a) J. Hertkorn, J.-N. Schmidt, M. Guo, F. Böttcher, K. S. H. Ng, S. D. Graham, P. Uerlings, H. P. Büchler, T. Langen, M. Zwierlein, and T. Pfau, Phys. Rev. Lett. 127, 155301 (2021a).
  * Gross (1957) E. P. Gross, Phys. Rev. 106, 161 (1957).
  * Gross (1958) E. P. Gross, Annals of Physics 4, 57 (1958).
  * Yang (1962) C. N. Yang, Rev. Mod. Phys. 34, 694 (1962).
  * Andreev and Lifshits (1969) A. Andreev and I. Lifshits, Zhur Eksper Teoret Fiziki 56, 2057 (1969).
  * Chester (1970) G. V. Chester, Phys. Rev. A 2, 256 (1970).
  * Leggett (1970) A. J. Leggett, Phys. Rev. Lett. 25, 1543 (1970).
  * Pomeau and Rica (1994) Y. Pomeau and S. Rica, Phys. Rev. Lett. 72, 2426 (1994).
  * Boninsegni and Prokof’ev (2012) M. Boninsegni and N. V. Prokof’ev, Rev. Mod. Phys. 84, 759 (2012).
  * Kim and Chan (2004a) E. Kim and M. H. W. Chan, Nature 427, 225 (2004a).
  * Kim and Chan (2004b) E. Kim and M. H. W. Chan, Science 305, 1941 (2004b).
  * Balibar (2010) S. Balibar, Nature 464, 176 (2010).
  * Kim and Chan (2012) D. Y. Kim and M. H. W. Chan, Phys. Rev. Lett. 109, 155301 (2012).
  * Chan  _et al._ (2013) M. H. W. Chan, R. B. Hallock, and L. Reatto, Journal of Low Temperature Physics 172, 317 (2013).
  * Henkel  _et al._ (2010) N. Henkel, R. Nath, and T. Pohl, Phys. Rev. Lett. 104, 195302 (2010).
  * Cinti  _et al._ (2010) F. Cinti, P. Jain, M. Boninsegni, A. Micheli, P. Zoller, and G. Pupillo, Phys. Rev. Lett. 105, 135301 (2010).
  * Saccani  _et al._ (2011) S. Saccani, S. Moroni, and M. Boninsegni, Phys. Rev. B 83, 092506 (2011).
  * Léonard  _et al._ (2017a) J. Léonard, A. Morales, P. Zupancic, T. Donner, and T. Esslinger, Science 358, 1415 (2017a).
  * Léonard  _et al._ (2017b) J. Léonard, A. Morales, P. Zupancic, T. Esslinger, and T. Donner, Nature 543, 87 (2017b).
  * Lin  _et al._ (2011) Y. J. Lin, K. Jiménez-García, and I. B. Spielman, Nature 471, 83 (2011).
  * Li  _et al._ (2016) J. Li, W. Huang, B. Shteynas, S. Burchesky, F. C. Top, E. Su, J. Lee, A. O. Jamison, and W. Ketterle, Phys. Rev. Lett. 117, 185301 (2016).
  * Li  _et al._ (2017) J.-R. Li, J. Lee, W. Huang, S. Burchesky, B. Shteynas, F. C. Top, A. O. Jamison, and W. Ketterle, Nature 543, 91–94 (2017).
  * Kadau  _et al._ (2016) H. Kadau, M. Schmitt, M. Wenzel, C. Wink, T. Maier, I. Ferrier-Barbut, and T. Pfau, Nature 530, 194–197 (2016).
  * Ferrier-Barbut  _et al._ (2016) I. Ferrier-Barbut, H. Kadau, M. Schmitt, M. Wenzel, and T. Pfau, Phys. Rev. Lett. 116, 215301 (2016).
  * Schmitt  _et al._ (2016) M. Schmitt, M. Wenzel, F. Böttcher, I. Ferrier-Barbut, and T. Pfau, Nature 539, 259 (2016).
  * Böttcher  _et al._ (2019b) F. Böttcher, M. Wenzel, J.-N. Schmidt, M. Guo, T. Langen, I. Ferrier-Barbut, T. Pfau, R. Bombín, J. Sánchez-Baena, J. Boronat, and F. Mazzanti, Phys. Rev. Research 1, 033088 (2019b).
  * Chomaz  _et al._ (2016) L. Chomaz, S. Baier, D. Petter, M. J. Mark, F. Wächtler, L. Santos, and F. Ferlaino, Phys. Rev. X 6, 041039 (2016).
  * Chomaz  _et al._ (2018) L. Chomaz, R. M. W. van Bijnen, D. Petter, G. Faraoni, S. Baier, J. H. Becher, M. J. Mark, F. Wächtler, L. Santos, and F. Ferlaino, Nature Phys. 14, 442–446 (2018).
  * Sohmen  _et al._ (2021) M. Sohmen, C. Politi, L. Klaus, L. Chomaz, M. J. Mark, M. A. Norcia, and F. Ferlaino, Phys. Rev. Lett. 126, 233401 (2021).
  * Schmidt  _et al._ (2021) J.-N. Schmidt, J. Hertkorn, M. Guo, F. Böttcher, M. Schmidt, K. S. H. Ng, S. D. Graham, T. Langen, M. Zwierlein, and T. Pfau, Phys. Rev. Lett. 126, 193002 (2021).
  * Biagioni  _et al._ (2022) G. Biagioni, N. Antolini, A. Alaña, M. Modugno, A. Fioretti, C. Gabbanini, L. Tanzi, and G. Modugno, Phys. Rev. X 12, 021019 (2022).
  * Bland  _et al._ (2022a) T. Bland, E. Poli, C. Politi, L. Klaus, M. A. Norcia, F. Ferlaino, L. Santos, and R. N. Bisset, Phys. Rev. Lett. 128, 195302 (2022a).
  * Lima and Pelster (2011) A. R. P. Lima and A. Pelster, Phys. Rev. A 84, 041604 (2011).
  * Wächtler and Santos (2016a) F. Wächtler and L. Santos, Phys. Rev. A 93, 061603 (2016a).
  * Wächtler and Santos (2016b) F. Wächtler and L. Santos, Phys. Rev. A 94, 043618 (2016b).
  * Bisset  _et al._ (2016) R. N. Bisset, R. M. Wilson, D. Baillie, and P. B. Blakie, Phys. Rev. A 94, 033619 (2016).
  * Petrov (2015) D. S. Petrov, Phys. Rev. Lett. 115, 155302 (2015).
  * Petrov and Astrakharchik (2016) D. S. Petrov and G. E. Astrakharchik, Phys. Rev. Lett. 117, 100401 (2016).
  * Jørgensen  _et al._ (2018) N. B. Jørgensen, G. M. Bruun, and J. J. Arlt, Phys. Rev. Lett. 121, 173403 (2018).
  * Cabrera  _et al._ (2017) C. R. Cabrera, L. Tanzi, J. Sanz, B. Naylor, P. Thomas, P. Cheiney, and L. Tarruell, Science 359, 301–304 (2017).
  * Semeghini  _et al._ (2018) G. Semeghini, G. Ferioli, L. Masi, C. Mazzinghi, L. Wolswijk, F. Minardi, M. Modugno, G. Modugno, M. Inguscio, and M. Fattori, Phys. Rev. Lett. 120, 235301 (2018).
  * Skov  _et al._ (2021) T. G. Skov, M. G. Skou, N. B. Jørgensen, and J. J. Arlt, Phys. Rev. Lett. 126, 230404 (2021).
  * Bisset and Blakie (2015) R. N. Bisset and P. B. Blakie, Phys. Rev. A 92, 061603 (2015).
  * Xi and Saito (2016) K.-T. Xi and H. Saito, Phys. Rev. A 93, 011604 (2016).
  * Blakie (2016) P. B. Blakie, Phys. Rev. A 93, 033644 (2016).
  * Saito (2016) H. Saito, Journal of the Physical Society of Japan 85, 053001 (2016).
  * Macia  _et al._ (2016) A. Macia, J. Sánchez-Baena, J. Boronat, and F. Mazzanti, Phys. Rev. Lett. 117, 205301 (2016).
  * Baillie  _et al._ (2017) D. Baillie, R. M. Wilson, and P. B. Blakie, Phys. Rev. Lett. 119, 255302 (2017).
  * Edler  _et al._ (2017) D. Edler, C. Mishra, F. Wächtler, R. Nath, S. Sinha, and L. Santos, Phys. Rev. Lett. 119, 050403 (2017).
  * Baillie and Blakie (2018) D. Baillie and P. B. Blakie, Phys. Rev. Lett. 121, 195301 (2018).
  * Roccuzzo and Ancilotto (2019) S. M. Roccuzzo and F. Ancilotto, Phys. Rev. A 99, 041601 (2019).
  * Zhang  _et al._ (2019) Y.-C. Zhang, F. Maucher, and T. Pohl, Phys. Rev. Lett. 123, 015301 (2019).
  * Blakie  _et al._ (2020) P. B. Blakie, D. Baillie, L. Chomaz, and F. Ferlaino, Phys. Rev. Research 2, 043318 (2020).
  * Mishra  _et al._ (2020) C. Mishra, L. Santos, and R. Nath, Phys. Rev. Lett. 124, 073402 (2020).
  * Chomaz (2020) L. Chomaz, Phys. Rev. A 102, 023333 (2020).
  * Poli  _et al._ (2021) E. Poli, T. Bland, C. Politi, L. Klaus, M. A. Norcia, F. Ferlaino, R. N. Bisset, and L. Santos, Phys. Rev. A 104, 063307 (2021).
  * Hertkorn  _et al._ (2021b) J. Hertkorn, J.-N. Schmidt, M. Guo, F. Böttcher, K. S. H. Ng, S. D. Graham, P. Uerlings, T. Langen, M. Zwierlein, and T. Pfau, Phys. Rev. Research 3, 033125 (2021b).
  * Zhang  _et al._ (2021) Y.-C. Zhang, T. Pohl, and F. Maucher, Phys. Rev. A 104, 013310 (2021).
  * Young-S. and Adhikari (2022) L. E. Young-S. and S. K. Adhikari, Phys. Rev. A 105, 033311 (2022).
  * Ghosh  _et al._ (2022) R. Ghosh, C. Mishra, L. Santos, and R. Nath, 2210.01093 (2022).
  * Gallemí and Santos (2022) A. Gallemí and L. Santos, Phys. Rev. A 106, 063301 (2022).
  * Schmidt  _et al._ (2022) M. Schmidt, L. Lassablière, G. Quéméner, and T. Langen, Phys. Rev. Research 4, 013235 (2022).
  * Halder  _et al._ (2022a) S. Halder, K. Mukherjee, S. I. Mistakidis, S. Das, P. G. Kevrekidis, P. K. Panigrahi, S. Majumder, and H. R. Sadeghpour, Phys. Rev. Research 4, 043124 (2022a).
  * Tengstrand  _et al._ (2021) M. N. Tengstrand, D. Boholm, R. Sachdeva, J. Bengtsson, and S. M. Reimann, Phys. Rev. A 103, 013313 (2021).
  * Roccuzzo  _et al._ (2022) S. M. Roccuzzo, A. Recati, and S. Stringari, Phys. Rev. A 105, 023316 (2022).
  * Ilzhöfer  _et al._ (2021) P. Ilzhöfer, M. Sohmen, G. Durastante, C. Politi, A. Trautmann, G. Natale, G. Morpurgo, T. Giamarchi, L. Chomaz, M. J. Mark, and F. Ferlaino, Nature Phys. 17, 356–361 (2021).
  * Norcia  _et al._ (2022) M. A. Norcia, E. Poli, C. Politi, L. Klaus, T. Bland, M. J. Mark, L. Santos, R. N. Bisset, and F. Ferlaino, Phys. Rev. Lett. 129, 040403 (2022).
  * Mukherjee and Reimann (2022) K. Mukherjee and S. M. Reimann, arXiv:2212.09658 (2022).
  * Roccuzzo  _et al._ (2020) S. M. Roccuzzo, A. Gallemí, A. Recati, and S. Stringari, Phys. Rev. Lett. 124, 045702 (2020).
  * Gallemí  _et al._ (2020) A. Gallemí, S. M. Roccuzzo, S. Stringari, and A. Recati, Phys. Rev. A 102, 023322 (2020).
  * Šindik  _et al._ (2022) M. Šindik, A. Recati, S. M. Roccuzzo, L. Santos, and S. Stringari, Phys. Rev. A 106, L061303 (2022).
  * Smith  _et al._ (2021) J. C. Smith, D. Baillie, and P. B. Blakie, Phys. Rev. Lett. 126, 025302 (2021).
  * Bisset  _et al._ (2021) R. N. Bisset, L. A. P. n. Ardila, and L. Santos, Phys. Rev. Lett. 126, 025301 (2021).
  * Li  _et al._ (2022) S. Li, U. N. Le, and H. Saito, Phys. Rev. A 105, L061302 (2022).
  * Scheiermann  _et al._ (2023) D. Scheiermann, L. A. P. n. Ardila, T. Bland, R. N. Bisset, and L. Santos, Phys. Rev. A 107, L021302 (2023).
  * Bland  _et al._ (2022b) T. Bland, E. Poli, L. A. P. n. Ardila, L. Santos, F. Ferlaino, and R. N. Bisset, Phys. Rev. A 106, 053322 (2022b).
  * Halder  _et al._ (2022b) S. Halder, S. Das, and S. Majumder, arXiv:2210.00598 (2022b).
  * Tang  _et al._ (2018) Y. Tang, W. Kao, K.-Y. Li, and B. L. Lev, Phys. Rev. Lett. 120, 230401 (2018).
  * Giovanazzi  _et al._ (2002) S. Giovanazzi, A. Görlitz, and T. Pfau, Phys. Rev. Lett. 89, 130401 (2002).
  * Baillie and Blakie (2020) D. Baillie and P. B. Blakie, Phys. Rev. A 101, 043606 (2020).
  * Prasad  _et al._ (2019) S. B. Prasad, T. Bland, B. C. Mulkerin, N. G. Parker, and A. M. Martin, Phys. Rev. Lett. 122, 050401 (2019).
  * Wenzel (2022) M. Wenzel, Macroscopic States of Dipolar Quantum Gases, Ph.D thesis, Stuttgart, Germany (2022).
  * Kirkby  _et al._ (2023) W. Kirkby, T. Bland, and B. R. N. Ferlaino, F., arXiv:2301.08007 (2023).
  * Bilitewski  _et al._ (2023) T. Bilitewski, G. A. Dominguez-Castro, D. Wellnitz, A. M. Rey, and L. Santos, arXiv:2302.09059 (2023).
  * Natale  _et al._ (2022) G. Natale, T. Bland, S. Gschwendtner, L. Lafforgue, D. S. Grün, A. Patscheider, M. J. Mark, and F. Ferlaino, Communications Physics (2022), 10.1038/s42005-022-01009-8.
  * Du  _et al._ (2023) L. Du, B. P., C. M., J. d. Hond, Y. Lu, and W. Ketterle, arXiv:2302.07209 (2023).
  * Wilson  _et al._ (2009a) R. M. Wilson, S. Ronen, and J. L. Bohn, Phys. Rev. A 79, 013621 (2009a).
  * Wilson  _et al._ (2009b) R. M. Wilson, S. Ronen, and J. L. Bohn, Phys. Rev. A 80, 023614 (2009b).
  * Lu  _et al._ (2010) H.-Y. Lu, H. Lu, J.-N. Zhang, R.-Z. Qiu, H. Pu, and S. Yi, Phys. Rev. A 82, 023622 (2010).
  * Blakie  _et al._ (2012) P. B. Blakie, D. Baillie, and R. N. Bisset, Phys. Rev. A 86, 021604 (2012).
  * Martin and Blakie (2012) A. D. Martin and P. B. Blakie, Phys. Rev. A 86, 053623 (2012).
  * Lima and Pelster (2012) A. R. P. Lima and A. Pelster, Phys. Rev. A 86, 063609 (2012).
  * Crank and Nicolson (1947) J. Crank and P. Nicolson, Math. Proc. Camb. Philos. Soc 43, 50–67 (1947).
  * Antoine  _et al._ (2013) X. Antoine, W. Bao, and C. Besse, Comput. Phys. Commun. 184, 2621 (2013).
  * Pethick and Smith (2002) C. J. Pethick and H. Smith,  _Bose-Einstein Condensation in Dilute Gases_ (Cambridge University Press, Cambridge, United Kingdom, 2002).
  * Stringari and Pitaevskii (2003) S. Stringari and L. Pitaevskii,  _Bose-Einstein Condensation_ (Oxford University Press, Oxford, United Kingdom, 2003).
  * Leggett (1998) A. J. Leggett, J. Stat. Phys. 93, 927 (1998).
  * Hertkorn  _et al._ (2021c) J. Hertkorn, J.-N. Schmidt, F. Böttcher, M. Guo, M. Schmidt, K. S. H. Ng, S. D. Graham, H. P. Büchler, T. Langen, M. Zwierlein, and T. Pfau, Phys. Rev. X 11, 011037 (2021c).
  * Bühler  _et al._ (2022) C. Bühler, T. Ilg, and H. P. Büchler, arXiv:2211.17251 (2022).
  * Nambu (2009) Y. Nambu, Rev. Mod. Phys. 81, 1015 (2009).
  * Pekker and Varma (2015) D. Pekker and C. Varma, Annual Review of Condensed Matter Physics 6, 269 (2015), https://doi.org/10.1146/annurev-conmatphys-031214-014350 .
  * Kohn (1961) W. Kohn, Phys. Rev. 123, 1242 (1961).
  * Klaus  _et al._ (2022) L. Klaus, T. Bland, E. Poli, C. Politi, G. Lamporesi, E. Casotti, R. N. Bisset, M. J. Mark, and F. Ferlaino, Nature Physics (2022), 10.1038/s41567-022-01793-8.
  * Bland  _et al._ (2023) T. Bland, G. Lamporesi, M. J. Mark, and F. Ferlaino, arXiv:2303.13263 (2023).
  * Arazo  _et al._ (2023) M. Arazo, A. Gallemí, M. Guilleumas, R. Mayol, and L. Santos, arXiv:2303.02087 (2023).
  * Sánchez-Baena  _et al._ (2023) J. Sánchez-Baena, C. Politi, F. Maucher, F. Ferlaino, and T. Pohl, Nature Communications 14 (2023), 10.1038/s41467-023-37207-3.


  * Citation

Mustafa, Ahmad, and Ghassan AlRegib. ”Explainable machine learning for
hydrocarbon prospect risking.” Second International Meeting for Applied
Geoscience & Energy. Society of Exploration Geophysicists and American
Association of Petroleum Geologists, 2022.

  * DOI

https://doi.org/10.1190/image2022-3752055.1

  * Review

Date of presentation: 31st August 2022

  * Data and Code

Data is proprietary and can not be shared.

  * Bib

@inproceedings{mustafa2022explainable,
title={Explainable machine learning for hydrocarbon prospect risking},
author={Mustafa, Ahmad and AlRegib, Ghassan},
booktitle={Second International Meeting for Applied Geoscience & Energy},
pages={1825–1829},
year={2022},
organization={Society of Exploration Geophysicists and American Association of
Petroleum Geologists} }

  * Contact

<EMAIL_ADDRESS>OR<EMAIL_ADDRESS>
https://ghassanalregib.info/

Example Mustafa & AlRegib Explainable Machine Learning for Hydrocarbon
Prospect Risking

# Explainable Machine Learning for Hydrocarbon Prospect Risking

Ahmad Mustafa11footnotemark: 1 and Ghassan AlRegib  Center for Energy and Geo
Processing (CeGP)  School of Electrical and Computer Engineering  Georgia
Institute of Technology  {amustafa9<EMAIL_ADDRESS>
##  1 Summary

Hydrocarbon prospect risking is a critical application in geophysics
predicting well outcomes from a variety of data including geological,
geophysical, and other information modalities. Traditional routines require
interpreters to go through a long process to arrive at the probability of
success of specific outcomes. AI has the capability to automate the process
but its adoption has been limited thus far owing to a lack of transparency in
the way complicated, black box models generate decisions. We demonstrate how
LIME—a model-agnostic explanation technique—can be used to inject trust in
model decisions by uncovering the model’s reasoning process for individual
predictions. It generates these explanations by fitting interpretable models
in the local neighborhood of specific datapoints being queried. On a dataset
of well outcomes and corresponding geophysical attribute data, we show how
LIME can induce trust in model’s decisions by revealing the decision-making
process to be aligned to domain knowledge. Further, it has the potential to
debug mispredictions made due to anomalous patterns in the data or faulty
training datasets.

##  2 Introduction

Progress in machine learning research has led to major advancements for
various computer vision and natural language processing tasks like image
classification, object detection and localization, natural scene
understanding, speech recognition, and machine translation, to name a few
(Szegedy et al.,, 2013; Krizhevsky et al.,, 2017; Rahnemoonfar et al.,, 2021;
Deng et al.,, 2013; Singh et al.,, 2017). It has also resulted in major
breakthroughs in the extremely critical field of medical image analysis,
helping in the detection of cellular structures, tissue segmentation, and in
diagnosing various diseases (Shen et al.,, 2017). In exploration geophysics,
it has now become a vital component in many seismic interpretation and
inversion tasks, as evidenced by the plethora of works in in salt body
delineation (Wang et al.,, 2015; Amin et al.,, 2017; Shafiq et al.,, 2017; Di
et al., 2018b, ), fault detection (Di and AlRegib,, 2019; Di et al., 2019a, ;
Di et al., 2019b, ), facies classification (Alaudah et al., 2019b, ; Alaudah
et al., 2019c, ), seismic attribute analysis (Long et al.,, 2018; Di et al.,
2018a, ; Alfarraj et al.,, 2018), structural similarity based seismic image
retrieval and segmentation (Alaudah et al., 2019a, ), and seismic inversion
(Alfarraj and AlRegib,, 2019; Mustafa et al.,, 2019; Mustafa and AlRegib,,
2020).

Critical to a more widespread adoption of AI in geophysical routines and
workflows is a deeper understanding of the model’s reasoning process as it
arrives at a particular decision. This is true of all areas heavy with
medical, legal, ethical, and monetary concerns. It is now common knowledge
that accuracy metrics alone seldom guarantee a model’s reliability under
unknown conditions. This can happen owing to several factors: a model could
overfit to training data but perform poorly on test data; a distributional
shift in the test data at run time could lead to significant degradation in
model performance, especially if it has been trained on data of a different
nature; or there could be spurious correlations present in the training data
leading to the model learning to predict on false causal features.

Hydrocarbon prospect risking Roden et al., (2012) is one such application
requiring interpreters to gather data from multiple sources and modalities to
process that and arrive at a decision: success or failure for a given oil
prospect. A single wrong prediction could mean the loss of several million
dollars. The DHI consortium aimed to gather best practices for the purposes of
risking prospects by comparing and correlating the outcomes on drilled
prospects from around the world to various features collected on those
prospects. Given the extremely large number of feature data collected and
interpreted for this purpose, the application stands to benefit greatly from
adoption of AI—subject to addressing the concerns raised above.

In this paper, we demonstrate how the adoption of a model-agnostic explanation
technique leads us to better understand the reasoning process behind a variety
of ML models as they arrive at decisions for prospect success or failure on a
dataset of well outcomes and input features. The method called LIME (Locally
Interpretable Model Explanations) works by fitting interpretable models in the
neighborhood of datapoints of interest (Ribeiro et al.,, 2016). LIME does not
assume any knowledge of the inner workings of the model being explained and
works only by querying the model at and around individual datapoints to
understand the local behavior of the decision boundary. As we demonstrate
later, this helps to quantify trust in individual model predictions by
correlating its reasoning process to prior geophysical knowledge and to debug
models in case of faulty behavior.

Our major contributions in this work are 1) training various machine learning
models for the task of well outcome prediction for hydrocarbon risking, 2)
using LIME to generate local explanations for various datapoints in the test
set, and 3) identifying and verifying important correlations between features
and well outcomes based on LIME explanations. To the best of our knowledge,
this is the first work of its kind applying explainable machine learning to a
geophysics application.

##  3 Theory

Given a black box machine learning model $f$ and a datapoint $x$, LIME
explains the prediction of $f$ over $x$ by fitting an interpretable model in
the local neighborhood of $x$. The local neighborhood is denoted by $\pi_{x}$.
The interpretable model, denoted by $g$, could be a logistic regression
classifier or a decision tree. For an explanation to be plausible, the fitted
model $g$ must be faithful to the decision boundary of the original model $f$
in the neighborhood of datapoint $x$. This property of the explanation is
termed _fidelity_ and is measured by a loss term (denoted $L(f,g,x)$)
penalizing the difference between the predictions over the datapoint $x$ for
both $f$ and $g$. To better explore the decision boundary of $f$ in the region
around $x$, LIME randomly samples new datapoints $x^{\prime}$ in the vicinity
of $x$ to obtain a new dataset ($Z$) containing both the original and sampled
datapoints. It then queries the model $f$ to obtain their respective
predictions. In the process to obtain the best fitted model $g$, the fidelity
term incorporates predictions from all the datapoints contained in $Z$, the
contribution from each being weighed according to its proximity to the
original data point of interest $x$. The fidelity term is then more
appropriately defined as $L(f,g,\pi_{x})$, where the latter refers to the fact
that $g$ is being fitted to a dataset containing the original and sampled
examples in the neighborhood $\pi_{x}$ and not just $x$ itself. Another
property controlled by LIME for the generated explanations is the
_interpretability_ of the fitted model $g$ in a term $\Omega(g)$. For the
family of logistic regression classifiers, this could be a term penalizing the
absolute norm of fitted coefficients (e.g., LASSO regression) to lead to more
interpretable insights in the sparse weights. The complete optimization
problem is then a sum of the objective functions measuring both fidelity and
interpretability as shown in

$\epsilon(x)=\operatorname*{arg\,min}_{g\in G}\quad L(f,g,\pi_{x})+\Omega(g),$
(1)

where $\epsilon(x)$ refers to the explanation generated for the datapoint $x$
as the weights of the model $g$ (from a family of models $G$) that best
minimizes the above cost function. A potential form $L$ could take is

$L(f,g,\pi_{x})=\sum_{z\in Z}\pi_{x}(z)(f(z)-g(z))^{2},$ (2)

where $z$ refers to the datapoints in set $Z$ (containing both the original
and sampled data) and $\pi_{x}$ is a weighting function that weighs the
contribution of each $z$ according to its distance from $x$. Points sampled
closer to $x$ would contribute more to the loss function as opposed to those
further away. An example of such a weighting function could be

$\pi_{x}(z)=\text{exp}(-\frac{||x-z||^{2}}{\sigma}),$ (3)

where $\sigma$ is a controlled hyperparameter. For more details, the reader is
encouraged to refer to the work by Ribeiro et al., (2016).

##  4 Methodology

The DHI consortium provides a dataset of drilled prospects around the world in
a variety of geologic settings. For each prospect, it collects—among other
things—data containing scores on how each prospect ranked on a variety of
seismic amplitude anomalies, termed DHIs (Direct Hydrocarbon Indicators). The
DHI risk assessment process as described in the work by Roden et al., (2012)
computes the final Pg (Probability of Geologic success) in a three-step
process incorporating 1) Initial Pg (estimated solely from geologic factors),
2) DHI Index (Pg estimated specifically from seismic amplitude anomalies
associated with DHIs), and 3) data quality scores for various geophysical
data. For all prospects pertaining to AVO class 3 geologic settings, we obtain
their outcomes measured as either success (1) or failure (0). Also
corresponding to each prospect, we obtain data on its initial Pg, DHI index,
final Pg, and 3 other attributes measuring data quality for various
geophysical modalities. These 6 attributes then form the feature data for each
prospect in a six dimensional vector ($x\in\mathbb{R}^{6}$) whereas the
prospect outcome forms the target class ($y\in[0,1]$). From the 258 such
instances available, we randomly divide such prospects into training and test
data with a 50/50 ratio. We then train a variety of machine learning models on
the training split including a logistic regression classifier, support vector
machines, and multilayer perceptrons. Afterwards, we use LIME as described in
the section above to generate explanations for various datapoints in the test
set. In the section below, we describe and elaborate on some of our findings.

##  5 Results

Figure 1 depicts the LIME explanations generated for a certain datapoint in
the test split for the logistic regression classifier. With the features as
selected above, the classifier achieved training and test set accuracies of
0.84 and 0.77, respectively. This tells us the classifier was able to
generalize to a reasonable extent. LIME is used to generate insights about the
model’s workings for specific data points. In this case, the classifier
strongly predicts the prospect to be a success (with probability 0.975) where
the ground truth is also success (class 1). In generating this prediction, it
weighs positively all the 6 features present. This makes sense from a
geophysical perspective since the better one scores on different criteria of
the prospect risking process, the more likely one’s chance of recovering
hydrocarbons. Moreover, the explanations reveal that scoring high on features
2 and 3 (DHI index and Final Pg, respectively) is much more important in
making this decision than for other features. This also makes sense from a
geophysical point of view since these attributes generally play a very
significant role in determining a particular prospect’s chances of success.
This finding is also corroborated by Figure 2 where we plot the feature
distributions for DHI Index, Final Pg, and Seismic data quality for both
successes and failures in the DHI data. As can be observed, the former two
features tend to be much more discriminatory compared to the last—something
also verified by LIME.

Figure 1: LIME explanations for logistic regression classifier. x-axis
indicates weight values fitted to interpretable LIME models for various
features. The values in parentheses next to feature names indicate their
respective values. Figure 2: Plotting feature distributions for various
geophysical markers in both success and failed prospects for DHI data.

Figure 3 further demonstrates the model’s inner workings for two different
test examples with negative ground-truths. In the first case (Figure 3a)), the
model strongly predicts the prospect to be a failure, weighing negatively all
features, but especially focusing on features 3 and 2 (Final Pg and DHI
Index). Once again, this corresponds very well with what we know from prior
geophysical knowledge. We also make a curious observation in Figure 3b); the
models predicts a weak negative and the LIME explanations reveal it does so
primarily on the back of weak values for features 2 and 3 but offsetting that
with stronger values for features 1, 6, and 5 (Initial Pg and data quality
scores), resulting in an overall weak negative.

Figure 3: LIME explanations for failure class datapoints in test set for two
cases: (a) strong negative and (b) a weak negative prediction.

We also generate LIME explanations similarly for an SVM model trained on a
50/50 split of the data between training and test sets. The results are
depicted in Figure 4. Once again, we achieve reasonably high training and test
accuracies of 0.84 and 0.75 respectively. Yet again, we see the SVM model fits
positively to all features (but especially to final Pg and DHI index) to
strongly predict positive a successful prospect. It underscores weak values in
the same features to strongly predict a negative a prospect with negative
ground truth. The bottom figure presents an interesting case study: the model
predicts a weak negative against a successful prospect owing to relatively
weak values for most features. It puts a positive weight on feature 1 (Initial
Pg) despite common sense dictating otherwise (i.e., low initial Pg to predict
a negative). For an interpreter, this would be a red flag for them to probe
further into the particular prospect to see if there could be anomalous
patterns present or to check for bugs in the model.

Figure 4: LIME explanations for SVM for a strong positive (top), a strong
negative prediction (middle), and a weak negative wrongly predicted (bottom).

Figure 5 generates explanations for a two layer MLP acting on a positive and
two negative examples from the test set. We obtain training and test
accuracies of 0.83 and 0.81, respectively. As before, we observe the model to
be learning to correlate primarily the DHI Index and Final Pg to the predicted
result, although sometimes other features can come into play to weaken the
prediction.

Figure 5: LIME explanations for an MLP for a strong positive (top) and two
negatives (middle and bottom).

##  6 Conclusion

Understanding the way complicated, black box models reason is central to
widespread adoption of machine learning methods in areas requiring critical
decision making having significant ethical, medical, legal, and monetary
impacts. One way to assess trust in ML models is via model-agnostic techniques
exploring their decision boundaries in local sub-regions around specific data
points of interest. LIME generates explanations by fitting interpretable
models to original and sampled data in these sub-regions. Applied to
hydrocarbon prospect risking, LIME explanations can be used to verify model
calibration by comparing its reasoning process to prior geophysical knowledge.
Accuracy metrics may be misleading for a variety of reasons including but not
limited to overfitting, data drift, data leakage, and spurious correlative
features, among other reasons. Therefore, having qualitative interpretable
explanations for model decisions over individual datapoints helps inject trust
in model predictions, and hence, the model as a whole.

By generating such explanations for different ML models for prospect risking,
we verified the model reasoning process to broadly match an interpreter’s by
looking for the presence or absence of relevant features to generate
predictions for prospect outcome. Our experiments revealed DHI Index and Final
Pg to matter the most when deciding the outcome on a given prospect. We also
observed how models at times made weak predictions owing to certain factors
pulling the explanations in opposite directions (e.g., Figure 3b)).

In addition to inducing trust in the model, the proposed workflow has the
potential to be used to identify important causal factors and correlations in
a much wider set of features than used for the experiments here. Moreover,
such explanations can be used to perform feature selection to choose only the
most relevant attributes for the task of prospect risking. By comparing
explanations across multiple models, one can select for the most reliable
model while still keeping testing/validation set accuracy in consideration.

## References

  * (1) Alaudah, Y., M. Alfarraj, and G. AlRegib, 2019a, Structure label prediction using similarity-based retrieval and weakly supervised label mapping: GEOPHYSICS, 84, V67–V79.
  * (2) Alaudah, Y., P. Michałowicz, M. Alfarraj, and G. AlRegib, 2019b, A machine-learning benchmark for facies classification: Interpretation, 7, SE175–SE187.
  * (3) Alaudah, Y., M. Soliman, and G. AlRegib, 2019c, Facies classification with weak and strong supervision: A comparative study: Presented at the SEG Technical Program Expanded Abstracts 2019, Society of Exploration Geophysicists.
  * Alfarraj et al., (2018) Alfarraj, M., Y. Alaudah, Z. Long, and G. AlRegib, 2018, Multiresolution analysis and learning for computational seismic interpretation: The Leading Edge, 37, 443–450.
  * Alfarraj and AlRegib, (2019) Alfarraj, M., and G. AlRegib, 2019, Semisupervised sequence modeling for elastic impedance inversion: Interpretation, 7, no. 3, SE237–SE249.
  * Amin et al., (2017) Amin, A., M. Deriche, M. A. Shafiq, Z. Wang, and G. AlRegib, 2017, Automated salt-dome detection using an attribute ranking framework with a dictionary-based classifier: Interpretation, 5, SJ61–SJ79.
  * Deng et al., (2013) Deng, L., G. Hinton, and B. Kingsbury, 2013, New types of deep neural network learning for speech recognition and related applications: An overview: 2013 IEEE international conference on acoustics, speech and signal processing, IEEE, 8599–8603.
  * (8) Di, H., M. Alfarraj, and G. AlRegib, 2018a, Three-dimensional curvature analysis of seismic waveforms and its interpretational implications: Geophysical Prospecting.
  * Di and AlRegib, (2019) Di, H., and G. AlRegib, 2019, Semi-automatic fault/fracture interpretation based on seismic geometry analysis: Geophysical Prospecting, 67, 1379–1391.
  * (10) Di, H., M. Shafiq, and G. AlRegib, 2018b, Multi-attribute k-means clustering for salt-boundary delineation from three-dimensional seismic data: Geophysical Journal International, 215, 1999–2007.
  * (11) Di, H., M. A. Shafiq, Z. Wang, and G. AlRegib, 2019a, Improving seismic fault detection by super-attribute-based classification: Interpretation, 7, SE251–SE267.
  * (12) ——–, 2019b, Improving seismic fault detection by super-attribute-based classification: Interpretation, 7, no. 3, SE251–SE267.
  * Krizhevsky et al., (2017) Krizhevsky, A., I. Sutskever, and G. E. Hinton, 2017, Imagenet classification with deep convolutional neural networks: Commun. ACM, 60, 84–90.
  * Long et al., (2018) Long, Z., Y. Alaudah, M. A. Qureshi, Y. Hu, Z. Wang, M. Alfarraj, G. AlRegib, A. Amin, M. Deriche, S. Al-Dharrab, and H. Di, 2018, A comparative study of texture attributes for characterizing subsurface structures in seismic volumes: Interpretation, 6, no. 4, T1055–T1066.
  * Mustafa et al., (2019) Mustafa, A., M. Alfarraj, and G. AlRegib, 2019, Estimation of acoustic impedance from seismic data using temporal convolutional network: SEG Technical Program Expanded Abstracts 2019, Society of Exploration Geophysicists, 2554–2558.
  * Mustafa and AlRegib, (2020) Mustafa, A., and G. AlRegib, 2020, Joint learning for seismic inversion: An acoustic impedance estimation case study: SEG Technical Program Expanded Abstracts 2020, Society of Exploration Geophysicists, 1686–1690.
  * Rahnemoonfar et al., (2021) Rahnemoonfar, M., T. Chowdhury, A. Sarkar, D. Varshney, M. Yari, and R. R. Murphy, 2021, Floodnet: A high resolution aerial imagery dataset for post flood scene understanding: IEEE Access, 9, 89644–89654.
  * Ribeiro et al., (2016) Ribeiro, M. T., S. Singh, and C. Guestrin, 2016, ”why should i trust you?”: Explaining the predictions of any classifier: Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, Association for Computing Machinery, 1135–1144.
  * Roden et al., (2012) Roden, R., M. Forrest, and R. Holeywell, 2012, Relating seismic interpretation to reserve/resource calculations: Insights from a dhi consortium: The Leading Edge, 31, 1066–1074.
  * Shafiq et al., (2017) Shafiq, M. A., Z. Wang, G. AlRegib, A. Amin, and M. Deriche, 2017, A texture-based interpretation workflow with application to delineating salt domes: Interpretation, 5, SJ1–SJ19.
  * Shen et al., (2017) Shen, D., G. Wu, and H.-I. Suk, 2017, Deep learning in medical image analysis: Annual review of biomedical engineering, 19, 221–248.
  * Singh et al., (2017) Singh, S. P., A. Kumar, H. Darbari, L. Singh, A. Rastogi, and S. Jain, 2017, Machine translation using deep learning: An overview: 2017 international conference on computer, communications and electronics (comptelix), IEEE, 162–167.
  * Szegedy et al., (2013) Szegedy, C., A. Toshev, and D. Erhan, 2013, Deep neural networks for object detection: Presented at the Advances in Neural Information Processing Systems, Curran Associates, Inc.
  * Wang et al., (2015) Wang, Z., T. Hegazy, Z. Long, and G. AlRegib, 2015, Noise-robust detection and tracking of salt domes in postmigrated volumes using texture, tensors, and subspace learning: Geophysics, 80, WD101–WD116.