#  Stirring the false vacuum via interacting quantized bubbles
on a 5564-qubit quantum annealer

Jaka Vodeb<EMAIL_ADDRESS>Jülich Supercomputing Centre, Institute for
Advanced Simulation, Forschungszentrum Jülich, 52425 Jülich, Germany  Jean-
Yves Desaules<EMAIL_ADDRESS>Institute of Science and
Technology Austria (ISTA), Am Campus 1, 3400 Klosterneuburg, Austria  School
of Physics and Astronomy, University of Leeds, Leeds LS2 9JT, UK  Andrew
Hallam  School of Physics and Astronomy, University of Leeds, Leeds LS2 9JT,
UK  Andrea Rava  Jülich Supercomputing Centre, Institute for Advanced
Simulation, Forschungszentrum Jülich, 52425 Jülich, Germany  RWTH Aachen
University, 52056 Aachen, Germany  Gregor Humar  Department of Complex Matter,
Jožef Stefan Institute, Jamova 39, 1000 Ljubljana, Slovenia  Department of
Physics, Faculty for Mathematics and Physics, Jadranska 19, University of
Ljubljana, 1000 Ljubljana, Slovenia  Dennis Willsch  Jülich Supercomputing
Centre, Institute for Advanced Simulation, Forschungszentrum Jülich, 52425
Jülich, Germany  Fengping Jin  Jülich Supercomputing Centre, Institute for
Advanced Simulation, Forschungszentrum Jülich, 52425 Jülich, Germany  Madita
Willsch  Jülich Supercomputing Centre, Institute for Advanced Simulation,
Forschungszentrum Jülich, 52425 Jülich, Germany  AIDAS, 52425 Jülich, Germany
Kristel Michielsen  Jülich Supercomputing Centre, Institute for Advanced
Simulation, Forschungszentrum Jülich, 52425 Jülich, Germany  RWTH Aachen
University, 52056 Aachen, Germany  AIDAS, 52425 Jülich, Germany  Zlatko Papić
<EMAIL_ADDRESS>School of Physics and Astronomy, University of Leeds,
Leeds LS2 9JT, UK

###### Abstract

False vacuum decay is a potential mechanism governing the evolution of the
early Universe, with profound connections to non-equilibrium quantum physics,
including quenched dynamics, the Kibble-Zurek mechanism, and dynamical
metastability. The non-perturbative character of the false vacuum decay and
the scarcity of its experimental probes make the effect notoriously difficult
to study, with many basic open questions, such as how the bubbles of true
vacuum form, move and interact with each other. Here we utilize a quantum
annealer with 5564 superconducting flux qubits to directly observe quantized
bubble formation in real time – the hallmark of false vacuum decay dynamics.
Moreover, we develop an effective model that describes the initial bubble
creation and subsequent interaction effects. We demonstrate that the effective
model remains accurate in the presence of dissipation, showing that our
annealer can access coherent scaling laws in driven many-body dynamics of 5564
qubits for over $1\mu$s, i.e., more than 1000 intrinsic qubit time units. This
work sets the stage for exploring late-time dynamics of the false vacuum at
computationally intractable system sizes, dimensionality, and topology in
quantum annealer platforms.

Nearly half a century ago, Coleman proposed the idea that our Universe may
have cooled down into a metastable “false vacuum” state after the Big Bang and
the time of tunneling to the ground state or “true vacuum” was estimated to be
comparable to the lifetime of the Universe [1]. The idea was then further
developed and applied to various cosmological observations and theories [2, 3,
4, 5, 6, 7, 8], with ongoing attempts to observe the signatures of false
vacuum decay in gravitational waves [9].

The dynamics of false vacuum decay are believed to consist of “bubbles” of
true vacuum forming in the background of false vacuum, where the size of a
bubble is determined by balancing the energy gain proportional to the bubble
volume and energy loss proportional to the bubble surface. Bubbles are
typically assumed to undergo isolated quantum tunneling events, growing
classically at a model-dependent speed [9]. The quantum process is difficult
to study due to the non-perturbative nature of the dynamics. To circumvent
this issue, early theoretical works have explored the possibility of directly
creating new Universes in a laboratory setting [10] and in engineered
platforms based on condensed matter systems [11]. With the advances in
ultracold atomic gases, certain aspects of the false-vacuum decay dynamics can
now be studied in table-top experiments [12].

Recently, there has been a flurry of interest in simulating quantum field
theories using synthetic platforms of ultracold atoms in optical lattices,
superconducting circuits, trapped ions and Rydberg atoms [13, 14, 15], with
different proposals addressing specifically the decay of the false vacuum [16,
17, 18, 19, 20, 21, 22]. Two main approaches to quantum simulation involve
either using quantum gates on a digital quantum computer to directly emulate
the quantum field theory in question, or setting up an analogous system that
exhibits a controllable first-order quantum phase transition, where it is
possible to initialize in the false vacuum. In this paper, we take the latter
approach and set up a quantum annealer with $5564$ superconducting flux
qubits, which had previously been used to study the spin glass transition [23]
and the Kibble-Zurek mechanism [24, 25, 26]. We arrange the qubits in a ring
by coupling them via ferromagnetic interactions in the presence of a
transverse magnetic field, thus realizing the quantum Ising model. By then
tuning the uniform longitudinal field, we initialize the system in the
metastable false vacuum state and observe the decay into the true vacuum. The
discrete nature of the qubit lattice gives us a direct window into the
quantized bubble creation, whereby a cascade of bubble sizes is seen to emerge
by tuning the longitudinal field. Moreover, the longitudinal field in the
quantum annealer exhibits intrinsic modulation throughout the decay, driving
the dynamics and extending the regime where we observe the same scaling laws
as in coherent quantum dynamics up to 1000 qubit time units.

Quench dynamics of the Ising chain in transverse and longitudinal fields have
recently attracted much interest due to the confinement effect imposed by the
longitudinal field [27, 28, 29, 30, 31]. The latter has direct implications
for false vacuum decay enabling analytic predictions of the decay rate [32,
33]. Our simulation targets a different regime where quantized bubbles
dominate the out-of-equilibrium dynamics, originally proposed in the context
of the generalized Kibble-Zurek effect [34]. This enables us to access false
vacuum decay dynamics beyond the initial bubble creation and into the
previously unexplored regime of interacting bubbles. In contrast to the
typical false vacuum decay mechanism [1, 9, 33], we find that a large
quantized bubble cannot spread in isolation. It is only through the
interaction of two neighboring bubbles that one bubble can enlarge itself by
reducing the size of the other. Once reduced to the smallest size of one
lattice site, the bubble can then move freely along the system. These results
suggest a new physical picture of the false vacuum dynamics as a heterogeneous
gas of bubbles, where the smallest “light” bubbles bounce around in the
background of larger “heavy” bubbles that directly interact with each other.

##  I Quantum simulation of false vacuum decay

We study the ferromagnetic quantum Ising model in transverse and longitudinal
fields on a ring with $N$ sites:

$\hat{H}=-J\sum_{j=1}^{N}\hat{\sigma}^{z}_{j}\hat{\sigma}^{z}_{j+1}-h_{x}\sum_{j=1}^{N}\hat{\sigma}^{x}_{j}-h_{z}\sum_{j=1}^{N}\hat{\sigma}^{z}_{j},$
(1)

where $\hat{\sigma}^{\alpha}$ are the Pauli matrices, $J>0$ is the
ferromagnetic interaction strength between nearest-neighbor spins, $h_{x}$ and
$h_{z}$ are the transverse and longitudinal fields, respectively. We apply
periodic boundary conditions by identifying spin
$\hat{\sigma}_{N{+}1}^{z}\equiv\hat{\sigma}_{1}^{z}$. The field $h_{x}$ drives
the quantum dynamics of the system, while $h_{z}$ imposes an energy bias
between the states $\ket{\uparrow}$ and $\ket{\downarrow}$.

In the regime $0\leq h_{x}\ll J$ and $h_{z}=0$, there are two degenerate
ground states $\ket{\uparrow...\uparrow}$ and $\ket{\downarrow...\downarrow}$.
When $h_{z}>0$, the $\ket{\uparrow...\uparrow}$ state becomes the ground or
true vacuum state and $\ket{\downarrow...\downarrow}$ a metastable or false
vacuum state, see Fig. 1a. By first setting $h_{z}>0$ and adiabatically
turning on $h_{x}$ to a small value $h_{x}\ll J$, we initialize the system in
the $\ket{\uparrow...\uparrow}$ product state. Then we induce a first-order
quantum phase transition by flipping the sign of $h_{z}$, swapping the true
and false vacuum states, and observe the dynamics for a time duration $t$.
Finally, we turn $h_{x}$ back to zero as fast as possible and measure the spin
configuration in the $\hat{\sigma}^{z}$ basis. Fig. 1b illustrates the
described protocol, while Fig. 1c shows the embedding of the spin chain in a
qubit array used in our quantum simulations. We note here that $h_{z}(t)$ was
determined experimentally through single-qubit measurements and exhibits large
modulation around the final target value after the flip. This modulation
extends up to $t\sim 0.75\mu s$ in the evolution time and it will play an
important role in the interpretation of our data.

Figure 1:  Realizing false vacuum decay on a quantum annealer. a The semi-
classical energy landscape $V$ as a function of magnetization $M$ of a
ferromagnetic Ising chain in a transverse $h_{x}$ and longitudinal $h_{z}$
field. The landscape exhibits a local metastable minimum dubbed as the false
vacuum, represented by the polarized $\ket{\uparrow\uparrow\cdots\uparrow}$
state. The global minimum or true vacuum is the other polarized
$\ket{\downarrow\downarrow\cdots\downarrow}$ state. The false vacuum decay
unfolds via creation of quantized true-vacuum bubbles of size $n$, determined
by the energy balance between the surface ($4J$) and volume energy
contributions ($2h_{z}n$). b False vacuum decay observation protocol. We
initialize all qubits in the $\ket{\uparrow\uparrow\cdots\uparrow}$ state by
setting $h_{z}>0$ and adiabatically switch $h_{x}$ from $0$ to a small value
$h_{x}\ll J$ over $10\mu s$ ($0.27\mu s$ is used in the plot for clarity
purposes). Then we flip the sign of $h_{z}$, swapping the true and false
vacuum states, and observing the dynamics for a time duration $t$. Finally, we
turn $h_{x}$ back to $0$ as fast as possible ($\geq 0.18$ $\mu s$) and measure
the spin configuration in the $\hat{\sigma}^{z}$ basis. This protocol is
repeated $1000$ times for each value of $t$. c Embedding of a $5564$-qubit
ring on the Pegasus graph of the $5614$-qubit device D-Wave
$Advantage\\_system5.4$. The Pegasus graph contains $15\times 15\times 3$
$8$-qubit Chimera cells with complete bipartite connectivity (colored crosses)
that are coupled by additional external and odd couplers (gray lines) [35],
such that each qubit is connected to 15 other qubits on average. Qubits within
the $8$-qubit cells are connected along randomly sampled one-dimensional
chains (inset). d Spin configurations measured in our quantum simulation. The
inner ring shows the initial false vacuum state comprised of $5564$ spins (for
clarity, only $1000$ out of $5564$ spins in a single configuration are shown).
The outer 3 rings show configurations measured at $h_{z}=-0.1,-0.5,-2$ with
$h_{z}$ decreasing radially. An example of a large $n=306$ quantized bubble
shown in purple highlights the extent of the observed bubble sizes. e
Magnetization $M$ heat profile versus time $t$ and longitudinal field
magnitude $h_{z}$ at transverse field strength $h_{x}=0.002$. The color scheme
is split into two separate linear scales, a larger scale from $-1$ to $0.999$
(bottom half) and a smaller scale from $0.999$ to $1$ (upper half). The
adiabatic dynamics and the $n=1$-bubble resonance are easily observed on the
larger scale, while the $n=2$-bubble resonance can only be resolved in the 4th
decimal of $M$, due to the decrease of the rate of dynamics by an order of
magnitude. The apparent resonance at $h_{z}=-4$ is identified with adiabatic
dynamics rather than bubble creation, where the system follows an
instantaneous ground state during the evolution.

Our quantum simulations are performed in the small $h_{x}\ll J$ regime, where
we can apply semiclassical intuition based on the diagonal part of the
Hamiltonian in the $z$-basis. In this case, it is useful to gather possible
configurations of the system into sectors with the same value of
magnetization, $M=\langle\sum_{i}\hat{\sigma}_{i}^{z}/N\rangle$, separated by
energy gaps determined by the value of $h_{z}$. For general values of $h_{z}$,
the initial $\ket{\uparrow...\uparrow}$ state stays an eigenstate in its own
$M$-sector after the $h_{z}$ sign flip and no dynamics of $M$ are observed.
This is due to the large energy separation between different $M$ sectors that
cannot be hybridized by a small $h_{x}$. However, for specific values of
$h_{z}=-2J/n$, where $n>0$ is an integer, the surface energy cost for flipping
a domain of $n$ spins, $4J$, is exactly balanced out by the volume energy
gain, $2h_{z}n$ [34]. Hence, an arbitrarily small $h_{x}$ is sufficient to
hybridize the classical computational basis states into eigenstates consisting
of a superposition of the $\ket{\uparrow...\uparrow}$ state and so-called
$n$-bubbles, i.e., domain walls in the background of
$\ket{\uparrow...\uparrow}$. For example,
$\ket{\uparrow\uparrow\uparrow\downarrow\downarrow\downarrow\uparrow\uparrow\uparrow}$
is a state with a single 3-bubble. Fig. 1d shows the spin configurations
measured in our quantum simulations with bubble sizes up to $306$ spins, which
is consistent with the theoretical prediction in which we can form
increasingly larger bubbles by decreasing $h_{z}$ according to $h_{z}=-2J/n$.
For these discrete values of $h_{z}$, the initial state is no longer an
eigenstate and undergoes nontrivial quantum dynamics, resulting in large
changes in $M$. Fig. 1e shows the observed $n=1$ and $n=2$ resonances, where
large changes in $M$ can be seen at $h_{z}=-2J$ and $h_{z}=-J$, respectively,
in contrast to other values of $h_{z}$.

We note that significant changes in $M$ can also be observed in Fig. 1e at
values of $h_{z}\approx-4J$, where no dynamical resonances are expected. Such
a large magnitude of $h_{z}$ leads to thermally assisted adiabatic dynamics
[36], where the system can follow the instantaneous ground state during time
evolution. The adiabatic theorem is applicable if the time scale of
Hamiltonian changes is slower than or comparable to
$t_{a}\propto\Delta_{min}^{-2}$, where $\Delta_{min}$ is the minimum gap
between the instantaneous ground and first excited state. In the case
$|h_{z}|\gg J,h_{x}$, the gap becomes large enough for the time scale of
$h_{z}(t)$ to match $t_{a}$. Therefore, no bubble creation takes place and the
spins turn simultaneously and in accordance with $h_{z}(t)$ initially,
changing the state from fully polarized and triggering more complex resonant
processes, see Supplementary Section 2.

##  II Observation of quantized bubbles and dynamical scaling laws

To ascertain which bubbles are involved in magnetization changes, we measured
the $n$-bubble density
$\lambda_{n}=(1/N)\sum_{i=1}^{N}\braket{\hat{P}^{\uparrow}_{i}[\prod_{j=1}^{n}\hat{P}^{\downarrow}_{i+j}]\hat{P}^{\uparrow}_{i+n+1}}$,
where $\hat{P}^{\sigma}=\ket{\sigma}\bra{\sigma}$ is a projector on the local
$\sigma=\uparrow,\downarrow$ spin state. Figures 2a-d show the detected
$n=1,2,3,4,5,6$-bubble resonances. We observe a strong suppression of all
other bubble sizes except the expected ones. According to the theoretical
analysis presented in Methods, the leading-order effective Hamiltonian
describing an $n$-bubble resonance is proportional to $h_{x}^{n}$. If we
assume $h_{x}<J$, 1-bubbles are the fastest, then 2-bubbles, etc., arbitrarily
slowing down the dynamics as $n$ increases. Figures 2a-d show that we need to
increase $h_{x}$ by at least two orders of magnitude to begin to observe hints
of higher resonances through low-density bubble formation, which is consistent
with the theoretical prediction.

Figure 2:  Observation of quantized bubbles. a-d Bubble density measurements
at $J=1$ and different $h_{z}$ magnitudes, with $h_{x}=0.002$, $t=2\mu s$ in
a,b, $h_{x}=0.05$, $t=1\mu s$ in c, and $h_{x}=0.1$, $t=1\mu s$ in d. The
bubble sizes $n=1,2,\ldots,6$ are seen to be dominant around their respective
resonances $h_{z}=-2J/n$, indicated by vertical dotted lines. Error bars
across the entire figure are smaller than the size of the symbols.

In a two-level approximation [34], tunneling events to different $n$-bubbles
can be thought of as Landau-Zener transitions, where the metastable state
$\ket{\uparrow...\uparrow}$ and an $n$-bubble state at the appropriate
resonant conditions are the two states involved in the anticrossing. According
to Landau-Zener theory, it follows that $n$-bubble density
$\lambda_{n}\propto\tau_{Q}h_{x}^{n}$ should be proportional to the product of
the time it takes for the Hamiltonian to traverse the anticrossing $\tau_{Q}$,
determined by $h_{z}(t)$ in our case, and the $n$-th power of $h_{x}$. Using
our single-qubit measurements we show that the time it takes for $h_{z}(t)$ to
reach zero during its sign flip is proportional to the square of its magnitude
($\tau_{Q}\propto|h_{z}|^{2}$), see Supplementary Section 3 for details. This
means that $\lambda_{n}(t)$ curves measured at different pause times between
the initialization and measurement ramp $t$ should collapse onto a single
curve if we multiply $t$ by $h_{z}^{2}$. Figure 3a shows that the
$\lambda_{2}$ curves indeed exhibit a collapse according to this law.

Figure 3:  Scaling laws for bubble dynamics. a 2-bubble density at
$h_{x}=0.002$ as a function of time at various $h_{z}$ magnitudes (color bar).
Inset shows the collapse of different curves when time is rescaled by
$h_{z}^{2}$ in accordance with Landau-Zener theory [34]. b-c Magnetization at
$h_{z}=-J$ resonance as a function of rescaled time $h_{x}^{2}t$, for
different values of $h_{x}$ indicated on the color bar. Both the measured
magnetization curves in b and the 3-spin Bloch-Redfield numerical simulations
in c follow the same $h_{x}^{2}$ scaling law, suggesting that the effective
Hamiltonian governing the dynamics is proportional to $h_{x}^{2}$. The inset
in b shows the raw data obtained on the quantum annealer without rescaling.
The unscaled results of the numerical simulations are shown in the inset of c,
where the black curves show the magnetization in the effective model
describing the $h_{z}=-J$ resonance (see Methods). Error bars across the
entire figure are smaller than the size of the symbols.

Nevertheless, in order to fully understand the dynamics in our quantum
annealer, it is necessary to account for all dominant processes and not only
the creation of bubbles. We will focus on the dynamics from the initial state
$\ket{\uparrow\cdots\uparrow}$, for which the creation of $n$-bubbles happens
at $h_{z}=-2J/n$. For each resonance, we have derived the corresponding
effective model using the Schrieffer-Wolff transformation [37] and we present
the effective Hamiltonians at leading orders in the Methods section. The
effective Hamiltonian describing the dynamics at the 2-bubble resonance at
$h_{z}=-J$ is proportional to $h_{x}^{2}$. Figure 3b shows magnetization
measurements taken at this resonance using the quantum annealer and how the
$M(t)$ curves collapse when scaling the time axis with $h_{x}^{2}$. Our
numerical emulation of the quantum annealer in Fig. 3c suggests that the
$h_{x}^{2}$ scaling law is the same as in the case of coherent quantum
dynamics. Figure 3b shows only the initial behavior of $M(t)$, which follows
the $h_{z}(t)$ modulation at later times; however, after the $h_{z}(t)$
modulation stops, an $h_{x}^{3}$ scaling law emerges as a consequence of
thermalization combined with a relatively slow quantum simulation measurement
ramp, see Supplementary Sections 6-9 for more details.

##  III Bubble interactions

The measured dynamics of different bubble sizes at $h_{z}=-2J$ resonance in
Fig. 4a is consistent with the picture that 1-bubbles remain approximately
quantized and do not grow with time. On the quantum annealer, this persists
until thermalization kicks in and 1-bubbles start to transform into 3- and
5-bubbles, with 2- and 4-bubbles remaining suppressed throughout the time
evolution. The exploration of this peculiar thermalization effect is beyond
the scope of this work, as thermalization and bubble interaction effects
cannot be easily separated from each other in our quantum annealer due to
decoherence effects. Nevertheless, we now argue that bubble interactions play
a crucial role in the dynamics at higher $n>1$ resonances if a system is
perfectly isolated from the environment. We will demonstrate this in the
framework of effective models, presented in Methods and previously mentioned
in the context of Fig. 3.

Figure 4:  a Time series measurements of the bubble density at the 1-bubble
resonance $h_{z}=-2J$ and $h_{x}=0.002$. During the initial $h_{z}(t)$
modulation, the profile of which is shown by the black curve on the right
axis, the 1-bubble density (color bar) governs the dynamics. After about $\sim
0.75\mu s$ (dashed line), thermalization effects take over by transforming
1-bubbles into 3- and 5-bubbles. b Measurement of the emergent blockade
$Q_{B}$ and magnetization $M$ (right axis) at $h_{x}=0.002$ and $t=0.38\mu s$,
plotted as a function of $h_{z}$ magnitude. The blockade condition is violated
(deviates from 0) only at $h_{z}$ values significantly off the 1-bubble
resonance ($h_{z}\leq-3.5$), accompanied by large changes in $M$. Near the
resonance ($h_{z}\approx-2$), even though large changes in $M$ occur, the
blockade condition is respected. c-d Dynamics at the resonance $h_{z}=-J$ with
fixed $h_{x}=0.0203$ and $h_{x}=0.002$, respectively. Matrix-product state
simulation with 100 qubits in c captures some of the key aspects of the data
obtained on the quantum annealer in d. The black dashed lines in c show the
sudden change in the number of 2-bubbles when $h_{z}(t)$ is exactly at the
resonance point. The inset magnifies the low-density regime, where only 1- and
3-bubbles can be seen. The increase of 1- and 3-bubbles at later times is
likely due to 2-bubbles interacting. The quantum simulation using
$h_{x}=0.002$ in d shows good agreement with the theoretical prediction in c.
For all the measured data, error bars are smaller than the symbol size. e
Matrix-product-state simulation of the dynamics after an instantaneous quench
from a product state shown at the bottom, containing two large bubbles
($n_{1}=23$ and $n_{2}=24$ spins) next to each other in a system with 50 spins
in total. The system undergoes coherent evolution with fixed $h_{z}=-1$ and
$h_{x}=0.02$, and the color bar shows the “bubble interface density”,
$\langle\hat{P}_{j-1}^{\downarrow}\hat{P}_{j}^{\uparrow}\hat{P}_{j+1}^{\downarrow}\rangle$,
on a log scale and for all sites $j$. The moving front corresponds to the two
bubbles exchanging $\downarrow$-spins and changing their sizes. The final
state at the end of the evolution is a quantum superposition, with one of the
classical configurations shown at the top.

At the $n=1$ resonance, 1-bubbles can be created with a rate proportional to
$h_{x}$. These 1-bubbles can then hop along the chain with a rate proportional
to $h_{x}^{2}/J$, but they cannot merge with each other to create large
bubbles. In fact, even when accounting for higher-order processes, there is
still no path to create $(n>1)$-bubbles when starting from the
$\ket{\uparrow\uparrow\uparrow\cdots}$ state. This dictates that, at the $n=1$
resonance, there can never be two $\downarrow$-spins next to each other. The
system therefore experiences an emergent kinetic constraint, reminiscent of
the Rydberg blockade phenomenon [38]. We quantify the blockade by measuring
the operator
$\hat{Q}_{B}=(1/N)\sum_{j}\hat{P}^{\downarrow}_{j}\hat{P}^{\downarrow}_{j+1}$,
which counts the density of neighboring $\downarrow$-spins. We expect
$\langle\hat{Q}_{B}\rangle$ to be strongly suppressed around $h_{z}=-2J$,
rising towards 0.5 in other dynamical settings. Figure 4b shows a good match
between these predictions and quantum simulation data. Meanwhile, total
magnetization strongly deviates from the initial value of $1$, showing that
the lack of neighboring excitations is not trivially due to frozen dynamics.

Now we consider higher-order resonances $h_{z}=-2J/n$, with $n>1$, where
bubbles contain $n$ spins and are created at a rate proportional to
$h_{x}^{n}/J^{n-1}$. Once these large bubbles are created, they cannot hop
around. However, bubbles can now take or give $\downarrow$-spins to
neighboring bubbles, allowing them to change size. This occurs with a rate
$\propto h_{x}^{2}/J$ which, for large $n$, is much faster than bubble
creation. This can also lead to a bubble shrinking all the way down to a
1-bubble, which can then hop along the chain. This means that, despite large
bubbles being stuck in place, information can still flow through inter-bubble
interactions and movement of 1-bubbles.

Our theoretical predictions imply that, in a quantum simulation tuned to a
$n>1$ resonance, the size of bubbles is not limited to $n$, even if the system
is perfectly isolated from the environment. This can be seen in a fully
coherent matrix-product state simulation of a system with $N=100$ spins at
$n=2$ resonance in Fig. 4c. While 2-bubbles dominate, 1- and 3-bubbles are
also visible. This is expected as they are produced by interactions of
2-bubbles. Qualitatively similar behavior is also seen in the quantum
simulation data in Fig. 4d.

The data in Figs. 4c-d also highlights another important property for $n>1$:
the number of 2-bubbles changes abruptly at some times, while staying
approximately constant during the rest of the simulation. The timings of
abrupt changes coincide exactly with $h_{z}(t)$ hitting the appropriate
resonant value, while the rest of the time the system is slightly away from
the resonance. This highlights the sensitivity to the detuning,
$\delta=h_{z}+2J/n$, which competes with $h_{x}^{n}/J^{n-1}$. As $h_{x}/J\ll
1$, even a small $\delta$ is enough to overpower the bubble creation terms for
$n>1$. As the detuning is a diagonal contribution, it leads to the suppression
of all dynamical processes, including bubble creation. This pattern of sudden
changes due to the fluctuation of $h_{z}$ is clearly captured in the numerical
simulation in Fig. 4c but it is also visible in the annealer data in Fig. 4d.
We note that this sensitivity to detuning is expected to be less strong for
$n=1$, as in that case $\delta$ only competes with $h_{x}$.

To further highlight the importance of bubble interactions, we have studied a
closed system with two large bubbles next to each other and essentially
occupying the entire system, see Fig. 4e. This setup leaves no room for new
bubbles to appear, and the only resonant process left is the exchange between
the two bubbles. We can then track the interface between them by measuring
$\hat{P}_{j-1}^{\downarrow}\hat{P}_{j}^{\uparrow}\hat{P}_{j+1}^{\downarrow}$,
which is plotted on a log scale in Fig. 4e at the $h_{z}=-J$ resonance. While
the interface density is 1 at a single location at time $t=0$ and zero
everywhere else, as time goes on the interface steadily delocalizes due to the
bubbles exchanging $\downarrow$-spins and thereby changing their sizes. We
expect similar behavior to hold at other $n>1$ resonances.

##  IV Discussion and outlook

We have performed a quantum simulation of the false vacuum decay and
identified its underlying mechanism – the formation of quantized bubbles of
true vacuum. The large size of our 5564-qubit quantum annealer allowed us to
observe considerable bubble sizes comprising up to 300 spin flips, confirming
the standard cosmological scenario where the size of the formed bubble is
determined by the competition between the volume energy gain and surface
energy loss. Our central finding is that interactions between bubbles are the
key next-order effect after bubble creation, and we have argued that their
understanding is crucial for a comprehensive description of of false vacuum
decay. Previous studies of false vacuum decay in quantum spin chains [32, 33]
have explored a different parameter regime where $h_{x}$ is not small and the
energy spectrum forms a continuum. While the possibility of resonances was
pointed out as a subleading effect [32], these analytical considerations still
assume a dilute bubble picture, neglecting interactions between bubbles. Our
results therefore call for a deeper understanding of the interaction effects
between bubbles, not only in microscopic models such as the one studied here,
but also in quantum field theory approaches, including cosmological models of
the Big Bang.

More broadly, our work showcases that current quantum annealing devices can be
useful in probing complex many-body dynamics, e.g., as demonstrated here
through external field modulation on a scale of $1000$ individual qubit time
units. Extensions of our model to two or three spatial dimensions with various
lattice topologies are, in principle, straightforward on the same type of
device, potentially reaching intractable computational complexity with a
multitude of implications. Let us mention a few examples of other interesting
non-equilibrium phenomena that can be accessed in the platform established
here. False vacuum decay, as a specific instance of a first-order quantum
phase transition, allows to probe generalizations of the Kibble-Zurek scaling
laws [39, 40, 11, 41] in such transitions, in particular comparing the
predictions for the rate of defect formation after crossing the transition
with quantum simulations. Quantum metastability – the cornerstone of the false
vacuum decay phenomenon – also underlies reaction rate theory [42, 43, 44, 45,
46, 47], allowing the use of quantum simulation for estimating the transition
rate of decay processes from a metastable minimum to a lower energy state in
the presence of temperature, which is challenging to compute otherwise. In the
regime of stronger longitudinal fields, confinement effects are expected to
become important, possibly localizing bubbles in space and giving rise to an
emergent prethermalization regime [48]. Finally, at the 1-bubble resonance,
our model displays an emergent kinetic constraint that maps exactly to the so-
called PXP model [49, 50] (see Methods), which hosts quantum many-body scars
[38, 51, 52], and possibly other types of ergodicity breaking, such as Hilbert
space fragmentation and many-body localization, in higher dimensions [53, 54].
This opens the way to probing non-ergodic phenomena in large systems in the
presence of dissipation and potentially new types of scars in constrained
models at other $n>1$ resonances.

##  V Methods

###  V.1 Quantum simulation on D-Wave’s quantum annealer

Our quantum simulations utilized D-Wave’s quantum annealing device
$Advantage\\_system5.4$, which features $N_{q}=5614$ qubits and is kept at a
cryostat temperature of $16.4\pm 0.1\mathrm{mK}$. The annealer implements the
Hamiltonian

$\displaystyle\hat{H}_{\mathrm{DW}}$
$\displaystyle=-\frac{A(s)}{2}\left(\sum_{i=1}^{N_{q}}\hat{\sigma}_{i}^{x}\right)+\frac{B(s)}{2}\left(g(t)\sum_{i=1}^{N_{q}}h_{i}\hat{\sigma}_{i}^{z}+\sum_{i<j}^{N_{q}}J_{ij}\hat{\sigma}_{i}^{z}\hat{\sigma}_{j}^{z}\right),$
(2)

where $\hat{\sigma}_{i}^{x,z}$ are the Pauli matrices for the $i$th qubit,
$h_{i}$ is the longitudinal external field at qubit $i$, $J_{ij}$ are the
couplings between qubits $i$ and $j$, which are non-zero and user-tunable only
if they are physically connected in the quantum processing unit (Fig. 1c).
$A(s)$ and $B(s)$ represent the energy scales of their respective terms and
are driven in time by the annealing schedule $s(t)$, which is linearly
interpolated from a series of user-specified points $[(t_{i},s_{i})]$.

Finding a ring embedding in a given graph is an instance of an NP-complete
Hamiltonian circuit problem [55]. We generate our ring embedding on $5564$
qubits of the $Advantage\\_system5.4$ graph by first connecting all $8$-qubit
Chimera cells in the Pegasus topology, see Fig. 1c. We start in the upper-left
corner and proceed horizontally, changing the horizontal direction at the end
of every row, until we reach the bottom-right corner. The chain of qubits
within each $8$-qubit Chimera cell is chosen along a random suitable path, see
Fig. 1c inset. The ring is closed by proceeding along the outer qubits at the
right and top edge of the graph (black part of the chain in Fig. 1c). This
procedure yields a ring of $5446$ qubits. We then iteratively add qubits to
the chain from the set of omitted remaining qubits by adding detours into the
ring until we obtain the final $5564$-qubit closed chain. We note that a few
of the qubits and couplers in the full Pegasus graph are not present on the
device due to fabrication defects; these are accounted for individually.

We are interested in probing the dynamics of $\hat{H}$ in Eq. (1) at a certain
value of $h_{x}$. Therefore, we specify the annealing schedule as
$[(0,1),(irt,s_{h_{x}}),(irt+pt,s_{h_{x}}),(irt+pt+mt,1)]$, where $s_{h_{x}}$
is obtained from the relation $h_{x}=A(s_{h_{x}})/B(s_{h_{x}})$, obtained by
rewriting $\hat{H}_{DW}$ as $(B(s)/2)\hat{H}$. We choose uniform $h_{i}=h$,
$J_{ij}=-1$ and define $h_{x}=A(s)/B(s)$, $h_{z}=-g(t)h$. At time $t=0$, we
specify the initial state for all qubits as the product state
$\ket{\uparrow...\uparrow}$. Then, within the initial ramp time $irt$, we
bring the system to the desired $h_{x}$ value, which drives the dynamics we
are interested in, and keep it constant for time $pt$. We replace $pt$ with
$t$ in all plots. Finally, we bring $h_{x}$ to $0$ within time $mt$, which
constitutes a measurement. Only after $h_{x}$ is brought back to $0$ is it
possible to read out the state of the qubits in the computational or
$\hat{\sigma}^{z}$ basis.

Typical time scales that we used on the D-Wave device are $irt=10\mu s$,
$mt=272ns$, and $pt$ ranging from $0$ to $2\mu s$. After the initial state
preparation, the system remains in the $\ket{\uparrow...\uparrow}$ state due
to the small values of $h_{x}$ compared to $h_{z}$. During the entire time
evolution, which lasts for a time $irt+pt+mt$, the system is subject to open
system dynamics, governed by two main effects; measurement by the environment
and thermalization. Our single spin measurements show that measurement by the
environment is dominant whenever the system is being driven by the
longitudinal external field $h_{z}$. Whenever $h_{z}$ becomes constant,
thermalization effects become more evident and are heavily dependent on the
value of $h_{x}$, which drives the quantum dynamics of the system – see
Supplementary Sections 7 and 9.

###  V.2 Simulations of thermalization dynamics

To capture thermalization effects on the system’s dynamics, we employed the
Bloch-Redfield master equation [56]

$\frac{d}{dt}\rho_{ab}(t)=-i\omega_{ab}\rho_{ab}(t)+\sum_{cd}^{sec}R_{abcd}\rho_{cd},$
(3)

where $\hat{\rho}$ is the density matrix of the system,
$\omega_{ab}\equiv\omega_{a}-\omega_{b}$ with $\omega_{a}=E_{a}/\hbar$ and
$E_{a}$ being the eigenenergies of the system. $sec$ denotes the secular
approximation, which states that we can neglect all fast-rotating terms in the
sum, and $R_{abcd}$ is the Bloch-Redfield tensor [56]

$\displaystyle R_{abcd}=$
$\displaystyle-\frac{1}{2\hbar^{2}}\sum_{\alpha,\beta}\Biggl{\\{}\delta_{bd}\sum_{n}A_{an}^{\alpha}A_{nc}^{\beta}S_{\alpha\beta}(\omega_{cn})-A_{ac}^{\alpha}A_{db}^{\beta}S_{\alpha\beta}(\omega_{ca})+\delta_{ac}\sum_{n}A_{dn}^{\alpha}A_{nb}^{\beta}S_{\alpha\beta}(\omega_{dn})-A_{ac}^{\alpha}A_{db}^{\beta}S_{\alpha\beta}(\omega_{db})\Biggr{\\}},$
(4)

where $A_{ab}^{\alpha}$ are the matrix elements in the system’s eigenbasis of
the operator $A^{\alpha}$ that couples bilinearly to the bath. Here we choose
$A^{\alpha}=\sigma_{\alpha}^{z}$, where $\alpha$ runs through all the spins of
the system.
$S_{\alpha\beta}(\omega)=\eta\omega\theta(\omega)\exp{(\omega/\omega_{c})}$ is
the noise power spectrum of the bath, chosen to be Ohmic in our case, where
$\theta(\omega)$ is the Heaviside step function, $\eta$ the coupling strength
of the system-bath coupling that ranges from $0.1$ to $0.2$ in our case, and
$\omega_{c}$ a cutoff frequency higher than any other relevant energy scale.

The numerical simulations in Figs. 4c,e were performed under the assumption of
a closed system using matrix-product state (MPS) formalism [57]. For
efficiency, the simulated system has open boundary conditions, but we discard
the boundary sites when computing observable expectation values in order to
minimize boundary effects. To reach the long times required for the
simulation, a 4th-order time-evolving block decimation (TEBD4) was used [58,
59]. For Fig. 4c, the timestep is $\delta t=0.01$ while the maximum MPS bond
dimension is $\chi=128$, which was never saturated during the simulation. For
Fig. 4e, the timestep is t=0.05 while the maximum bond dimension is
$\chi=200$.

###  V.3 Effective models at different resonances

To fully understand the dynamics beyond bubble creation in the vicinity of
resonances, we have derived the corresponding effective Hamiltonians using the
Schrieffer-Wolff transformation [37]. We quote the main results here, while
the derivation and detailed analysis of the models are provided in
Supplementary Section 5. For $n=1$, in the sector containing the state
$\ket{\uparrow\cdots\uparrow}$, the combined effective Hamiltonian at first
and second order reads:

$\hat{H}^{(1,2)}_{\mathrm{eff},n{=}1}=-h_{x}\sum_{j=1}^{N}\hat{P}^{\uparrow}_{j-1}\hat{\sigma}^{x}_{j}\hat{P}^{\uparrow}_{j+1}-\delta\sum_{j=1}^{N}\hat{\sigma}^{z}_{j}+\frac{h_{x}^{2}}{4J}\Big{[}\sum_{j=1}^{N}\hat{P}^{\uparrow}_{j-1}(\hat{\sigma}^{+}_{j}\hat{\sigma}^{-}_{j+1}+\hat{\sigma}^{-}_{j}\hat{\sigma}^{+}_{j+1})\hat{P}^{\uparrow}_{j+2}+2\sum_{j=1}^{N}\hat{P}^{\downarrow}_{j}-\frac{3}{2}\sum_{j=1}^{N}\hat{P}^{\downarrow}_{j-1}\hat{P}^{\downarrow}_{j+1}\Big{]},$
(5)

where $\delta=h_{z}+2J$ is the (weak) detuning away from the $n=1$ resonance,
$\hat{\sigma}^{\pm}=(\hat{\sigma}^{x}\pm i\hat{\sigma}^{y})/2$ are the
standard spin raising and lowering operators, and
$\hat{P}^{\downarrow}=\ket{\downarrow}\bra{\downarrow}$,
$\hat{P}^{\uparrow}=\ket{\uparrow}\bra{\uparrow}$ are local spin projectors.

The dynamics generated by the Hamiltonian in Eq. (5) can be understood as
follows. The $\hat{P}^{\uparrow}\hat{\sigma}^{x}\hat{P}^{\uparrow}$ term
allows the creation of single-site bubbles (i.e., single $\downarrow$-spins in
a background of $\uparrow$-spins), while the
$\hat{P}^{\uparrow}\left(\hat{\sigma}^{+}\hat{\sigma}^{-}+\hat{\sigma}^{-}\hat{\sigma}^{+}\right)\hat{P}^{\uparrow}$
allows these bubbles to hop around. A sequence of allowed processes is
illustrated in Fig. 5. Importantly, due to the projectors, the bubbles cannot
merge to form larger ones. This is also impossible to do using higher-order
processes. A simple argument is that there are no states with larger bubbles
at the same classical energy (i.e., the energy contribution of the
$\hat{\sigma}^{z}$ terms) as the $\ket{\uparrow\uparrow\uparrow\cdots}$ state,
so it is impossible to reach these states resonantly.

In the main text, we have demonstrated that one measurable consequence of the
effective Hamiltonian in Eq. (5) is a robust emergent kinetic constraint
reminiscent of the Rydberg blockade [38]. The quality of this emergent
blockade can be assessed using the operator $\hat{Q}_{B}$ introduced in the
main text, which measures the density of neighboring $\downarrow$-spins and
can be equivalently expressed in the spin language as
$\hat{Q}_{B}=1/4+(1/(4N))\sum_{j}\hat{\sigma}^{z}_{j}\hat{\sigma}^{z}_{j+1}-(1/(2N))\sum_{j}\hat{\sigma}^{z}_{j}$.

Figure 5: Schematic of the allowed dynamical processes at $n=1$ and $n=2$
resonances. At $n=1$, the 1-bubbles can be directly created and move along the
chain. However, they cannot merge and their size is fixed. By contrast, for
$n>1$ (here illustrated for $n=2$) only $n$-bubbles can form. While these
larger bubbles cannot move, neighboring bubbles can exchange
$\downarrow$-spins to change their size. If 1-bubbles are formed that way,
they can hop along the chain unlike bigger bubbles.

For $n>1$ resonances, the bubble creation term is no longer dominant as it
happens at order $n$ according to

$\hat{H}_{\mathrm{eff},n}=c_{n}\frac{h_{x}^{n}}{J^{n-1}}\sum_{j=1}^{N}\hat{P}^{\uparrow}_{j}\left(\prod_{k=1}^{n}\hat{\sigma}^{-}_{j+k}\right)\hat{P}^{\uparrow}_{j+n+1}+\mathrm{h.c.},$
(6)

where $c_{n}$ is a coefficient that depends on the multiple subprocesses
involved, e.g., we have $c_{2}=-1$ and $c_{3}=-81/64$. Instead, regardless of
$n$ there are always other terms at order one and two that read

$\displaystyle\hat{H}^{(1,2)}_{\mathrm{eff},n{=}2}$
$\displaystyle=-\delta\sum_{j=1}^{N}\hat{\sigma}^{z}_{j}+\frac{h_{x}^{2}n}{4J}\sum_{j=1}^{N}\Big{(}\frac{\hat{P}^{\downarrow}_{j-1}\hat{\sigma}^{z}_{j}\hat{P}^{\downarrow}_{j+1}}{n+1}{+}\hat{P}^{\uparrow}_{j-1}\hat{\sigma}^{z}_{j}\hat{P}^{\downarrow}_{j+1}{+}\hat{P}^{\downarrow}_{j-1}\hat{\sigma}^{z}_{j}\hat{P}^{\uparrow}_{j+1}{-}\frac{\hat{P}^{\uparrow}_{j-1}\hat{\sigma}^{z}_{j}\hat{P}^{\uparrow}_{j+1}}{n-1}\Big{)}$
(7)
$\displaystyle{+}\frac{h_{x}^{2}n^{2}}{4J(n-1)}\sum_{j=1}^{N}\hat{P}^{\uparrow}_{j{-}1}\left(\hat{\sigma}^{+}_{j}\hat{\sigma}^{-}_{j{+}1}{+}\hat{\sigma}^{-}_{j}\hat{\sigma}^{+}_{j{+}1}\right)\hat{P}^{\uparrow}_{j{+}2}{-}\frac{h_{x}^{2}n^{2}}{4J(n+1)}\sum_{j=1}^{N}\hat{P}^{\downarrow}_{j-1}(\hat{\sigma}^{+}_{j}\hat{\sigma}^{-}_{j+1}+\hat{\sigma}^{-}_{j}\hat{\sigma}^{+}_{j+1})\hat{P}^{\downarrow}_{j+2}.$

The terms on the second line create different dynamics. The first one leads to
1-bubbles hopping (as for $n=1$), while the second one allows larger bubbles
to exchange down-spins in order to grow or shrink. This allows bubbles of size
other than $n$ to develop. This includes the case of larger bubbles shrinking
all the way down to 1-bubble and then moving on their own. A sequence of these
allowed processes is illustrated in Fig. 5 for $n=2$, with higher $n$ values
displaying qualitatively similar behaviors. The dynamics at $n>1$ resonances
are clearly much richer than at $n=1$. Indeed, while the bubble interaction
term should also be present for $n=1$, it cannot act between two 1-bubbles.
This would require one of them to shrink to 0, which is not resonant. Thus, in
the sector of the $\ket{\uparrow\cdots\uparrow}$ state where only 1-bubbles
appear, the bubble interaction term vanishes.

Finally, it is worth noting that the first term of the effective Hamiltonian
at the $n=1$ resonance [Eq. (5)], up to a global spin flip, is identical to
the PXP model used to describe chains of Rydberg atoms [49, 50]. The second
term can then be recast as $-2\delta\sum_{j}\hat{P}^{\downarrow}_{j}$ up to an
irrelevant constant, and then becomes the chemical potential for the effective
PXP model. On the other hand, to the best of our knowledge, the effective
Hamiltonians for $n>1$ resonances, Eq. (7), do not map to the models
previously studied in the literature.

##  VI Acknowledgements

J.V., D.W. and M.W. acknowledge support from the project Jülich UNified
Infrastructure for Quantum computing (JUNIQ) that has received funding from
the German Federal Ministry of Education and Research (BMBF) and the Ministry
of Culture and Science of the State of North Rhine-Westphalia. A.R.
acknowledges support from the project HPCQS (101018180) of the European High-
Performance Computing Joint Undertaking (EuroHPC JU). J.-Y.D., A.H. and Z.P.
acknowledge support by the Leverhulme Trust Research Leadership Award
RL-2019-015 and EPSRC grants EP/R513258/1, EP/W026848/1. J.-Y.D. acknowledges
support from the European Union’s Horizon 2020 research and innovation
programme under the Marie Skłodowska-Curie Grant Agreement No.101034413. This
research was supported in part by grants NSF PHY-1748958 and PHY-2309135 to
the Kavli Institute for Theoretical Physics (KITP). Computational portions of
this research work were carried out on ARC3 and ARC4, part of the High-
Performance Computing facilities at the University of Leeds, UK. G.H. would
like to acknowledge the financial support from ARIS, P1-0040 Nonequilibrium
Quantum System Dynamics. The authors gratefully acknowledge the Jülich
Supercomputing Centre (https://www.fzjuelich.de/ias/jsc) for funding this
project by providing computing time on the D-Wave Advantage™ System JUPSI
through the Jülich UNified Infrastructure for Quantum computing (JUNIQ). We
would like to acknowledge the helpful theoretical discussions with Gianluca
Lagnese and the quantum simulation related discussions with D-Wave’s
experimental team, in particular Allison MacDonald, Gabriel Poulin-Lamarre,
Axel Daian and Andrew Berkley. We also thank Victoria Goliber and Andy Mason
for patiently organizing and mediating the corresponding meetings that enabled
the discussions with D-Wave’s team.

##  VII Author Contributions Statement

J.V. conceptualized and performed the quantum simulations. J.V., D.W., F.J.
and G.H. designed and analyzed the quantum simulations, whose coherent
emulation was performed by A.R. and M.W. J-Y.D., A.H. and Z.P. conducted the
theoretical analysis. J.V., J-Y.D., D.W. and Z.P. co-wrote the manuscript with
input from other authors. All authors participated in the discussions of the
results and development of the manuscript.

##  VIII Competing Interests Statement

The authors declare no competing interests.

## References

  * Coleman [1977] S. Coleman, Fate of the false vacuum: Semiclassical theory, Phys. Rev. D 15, 2929 (1977).
  * Kobsarev  _et al._ [1974] I. Kobsarev, L. B. Okun, and M. V. Voloshin,  _Bubbles in metastable vacuum_, Tech. Rep. (Moscow Institute for Theoretical and Experimental Physics, USSR, 1974) iTEP–81.
  * Linde [1981] A. Linde, Fate of the false vacuum at finite temperature: Theory and applications, Physics Letters B 100, 37 (1981).
  * Guth [1981] A. H. Guth, Inflationary universe: A possible solution to the horizon and flatness problems, Phys. Rev. D 23, 347 (1981).
  * Hawking and Moss [1982] S. Hawking and I. Moss, Supercooled phase transitions in the very early universe, Physics Letters B 110, 35 (1982).
  * Abdalla  _et al._ [2022] E. Abdalla, G. F. Abellán, A. Aboubrahim, A. Agnello, Ö. Akarsu, Y. Akrami, G. Alestas, D. Aloni, L. Amendola, L. A. Anchordoqui, _et al._ , Cosmology intertwined: A review of the particle physics, astrophysics, and cosmology associated with the cosmological tensions and anomalies, Journal of High Energy Astrophysics 34, 49 (2022).
  * Isidori  _et al._ [2001] G. Isidori, G. Ridolfi, and A. Strumia, On the metastability of the Standard Model vacuum, Nuclear Physics B 609, 387 (2001).
  * Degrassi  _et al._ [2012] G. Degrassi, S. Di Vita, J. Elias-Miro, J. R. Espinosa, G. F. Giudice, G. Isidori, and A. Strumia, Higgs mass and vacuum stability in the standard model at NNLO, Journal of High Energy Physics 2012, 1 (2012).
  * Caprini  _et al._ [2020] C. Caprini, M. Chala, G. C. Dorsch, M. Hindmarsh, S. J. Huber, T. Konstandin, J. Kozaczuk, G. Nardini, J. M. No, K. Rummukainen, _et al._ , Detecting gravitational waves from cosmological phase transitions with LISA: an update, Journal of Cosmology and Astroparticle Physics 2020 (03), 024.
  * Farhi  _et al._ [1990] E. Farhi, A. H. Guth, and J. Guven, Is it possible to create a universe in the laboratory by quantum tunneling?, Nuclear Physics B 339, 417 (1990).
  * Zurek [1996] W. H. Zurek, Cosmological experiments in condensed matter systems, Physics Reports 276, 177 (1996).
  * Zenesini  _et al._ [2024] A. Zenesini, A. Berti, R. Cominotti, C. Rogora, I. G. Moss, T. P. Billam, I. Carusotto, G. Lamporesi, A. Recati, and G. Ferrari, False vacuum decay via bubble formation in ferromagnetic superfluids, Nature Physics 20, 558 (2024).
  * Bañuls  _et al._ [2020] M. C. Bañuls, R. Blatt, J. Catani, A. Celi, J. I. Cirac, M. Dalmonte, L. Fallani, K. Jansen, M. Lewenstein, S. Montangero, C. A. Muschik, B. Reznik, E. Rico, L. Tagliacozzo, K. Van Acoleyen, F. Verstraete, U.-J. Wiese, M. Wingate, J. Zakrzewski, and P. Zoller, Simulating lattice gauge theories within quantum technologies, The European Physical Journal D 74, 165 (2020).
  * Bauer  _et al._ [2023] C. W. Bauer, Z. Davoudi, A. B. Balantekin, T. Bhattacharya, M. Carena, W. A. de Jong, P. Draper, A. El-Khadra, N. Gemelke, M. Hanada, D. Kharzeev, H. Lamm, Y.-Y. Li, J. Liu, M. Lukin, Y. Meurice, C. Monroe, B. Nachman, G. Pagano, J. Preskill, E. Rinaldi, A. Roggero, D. I. Santiago, M. J. Savage, I. Siddiqi, G. Siopsis, D. Van Zanten, N. Wiebe, Y. Yamauchi, K. Yeter-Aydeniz, and S. Zorzetti, Quantum simulation for high-energy physics, PRX Quantum 4, 027001 (2023).
  * Halimeh  _et al._ [2023] J. C. Halimeh, M. Aidelsburger, F. Grusdt, P. Hauke, and B. Yang, Cold-atom quantum simulators of gauge theories (2023), arXiv:2310.12201 [cond-mat.quant-gas] .
  * Billam  _et al._ [2019] T. P. Billam, R. Gregory, F. Michel, and I. G. Moss, Simulating seeded vacuum decay in a cold atom system, Phys. Rev. D 100, 065016 (2019).
  * Billam  _et al._ [2020] T. P. Billam, K. Brown, and I. G. Moss, Simulating cosmological supercooling with a cold-atom system, Phys. Rev. A 102, 043324 (2020).
  * Abel and Spannowsky [2021] S. Abel and M. Spannowsky, Quantum-field-theoretic simulation platform for observing the fate of the false vacuum, PRX Quantum 2, 010349 (2021).
  * Ng  _et al._ [2021] K. L. Ng, B. Opanchuk, M. Thenabadu, M. Reid, and P. D. Drummond, Fate of the false vacuum: Finite temperature, entropy, and topological phase in quantum simulations of the early universe, PRX Quantum 2, 010350 (2021).
  * Milsted  _et al._ [2022] A. Milsted, J. Liu, J. Preskill, and G. Vidal, Collisions of false-vacuum bubble walls in a quantum spin chain, PRX Quantum 3, 020316 (2022).
  * Darbha  _et al._ [2024a] S. Darbha, M. Kornjača, F. Liu, J. Balewski, M. R. Hirsbrunner, P. Lopes, S.-T. Wang, R. V. Beeumen, D. Camps, and K. Klymko, False vacuum decay and nucleation dynamics in neutral atom systems (2024a), arXiv:2404.12360 [quant-ph] .
  * Darbha  _et al._ [2024b] S. Darbha, M. Kornjača, F. Liu, J. Balewski, M. R. Hirsbrunner, P. Lopes, S.-T. Wang, R. V. Beeumen, K. Klymko, and D. Camps, Long-lived oscillations of false and true vacuum states in neutral atom systems (2024b), arXiv:2404.12371 [quant-ph] .
  * Harris  _et al._ [2018] R. Harris, Y. Sato, A. J. Berkley, M. Reis, F. Altomare, M. H. Amin, K. Boothby, P. Bunyk, C. Deng, C. Enderud, S. Huang, E. Hoskinson, M. W. Johnson, E. Ladizinsky, N. Ladizinsky, T. Lanting, R. Li, T. Medina, R. Molavi, R. Neufeld, T. Oh, I. Pavlov, I. Perminov, G. Poulin-Lamarre, C. Rich, A. Smirnov, L. Swenson, N. Tsai, M. Volkmann, J. Whittaker, and J. Yao, Phase transitions in a programmable quantum spin glass simulator, Science 361, 162 (2018).
  * Bando  _et al._ [2020] Y. Bando, Y. Susa, H. Oshiyama, N. Shibata, M. Ohzeki, F. J. Gómez-Ruiz, D. A. Lidar, S. Suzuki, A. del Campo, and H. Nishimori, Probing the universality of topological defect formation in a quantum annealer: Kibble-zurek mechanism and beyond, Phys. Rev. Res. 2, 033369 (2020).
  * King  _et al._ [2022] A. D. King, S. Suzuki, J. Raymond, A. Zucca, T. Lanting, F. Altomare, A. J. Berkley, S. Ejtemaee, E. Hoskinson, S. Huang, E. Ladizinsky, A. J. R. MacDonald, G. Marsden, T. Oh, G. Poulin-Lamarre, M. Reis, C. Rich, Y. Sato, J. D. Whittaker, J. Yao, R. Harris, D. A. Lidar, H. Nishimori, and M. H. Amin, Coherent quantum annealing in a programmable 2,000 qubit ising chain, Nature Physics 18, 1324 (2022).
  * King  _et al._ [2023] A. D. King, J. Raymond, T. Lanting, R. Harris, A. Zucca, F. Altomare, A. J. Berkley, K. Boothby, S. Ejtemaee, C. Enderud, E. Hoskinson, S. Huang, E. Ladizinsky, A. J. R. MacDonald, G. Marsden, R. Molavi, T. Oh, G. Poulin-Lamarre, M. Reis, C. Rich, Y. Sato, N. Tsai, M. Volkmann, J. D. Whittaker, J. Yao, A. W. Sandvik, and M. H. Amin, Quantum critical dynamics in a 5,000-qubit programmable spin glass, Nature 617, 61 (2023).
  * Kormos  _et al._ [2017] M. Kormos, M. Collura, G. Takács, and P. Calabrese, Real-time confinement following a quantum quench to a non-integrable model, Nat. Phys. 13, 246 (2017).
  * Liu  _et al._ [2019] F. Liu, R. Lundgren, P. Titum, G. Pagano, J. Zhang, C. Monroe, and A. V. Gorshkov, Confined quasiparticle dynamics in long-range interacting quantum spin chains, Phys. Rev. Lett. 122, 150601 (2019).
  * Tan  _et al._ [2021] W. L. Tan, P. Becker, F. Liu, G. Pagano, K. Collins, A. De, L. Feng, H. Kaplan, A. Kyprianidis, R. Lundgren, _et al._ , Domain-wall confinement and dynamics in a quantum simulator, Nat. Phys. 17, 742 (2021).
  * Vovrosh and Knolle [2021] J. Vovrosh and J. Knolle, Confinement and entanglement dynamics on a digital quantum computer, Scientific reports 11, 11577 (2021).
  * Lagnese  _et al._ [2023] G. Lagnese, F. M. Surace, S. Morampudi, and F. Wilczek, Detecting a long lived false vacuum with quantum quenches (2023), arXiv:2308.08340 [cond-mat.stat-mech] .
  * Rutkevich [1999] S. B. Rutkevich, Decay of the metastable phase in $d=1$ and $d=2$ ising models, Phys. Rev. B 60, 14525 (1999).
  * Lagnese  _et al._ [2021] G. Lagnese, F. M. Surace, M. Kormos, and P. Calabrese, False vacuum decay in quantum spin chains, Phys. Rev. B 104, L201106 (2021).
  * Sinha  _et al._ [2021] A. Sinha, T. Chanda, and J. Dziarmaga, Nonadiabatic dynamics across a first-order quantum phase transition: Quantized bubble nucleation, Phys. Rev. B 103, L220302 (2021).
  * D-Wave Systems [2020] D-Wave Systems,  _Technical Description of the D-Wave Quantum Processing Unit_, Tech. Rep. (D-Wave Systems Inc., Burnaby, BC, Canada, 2020) D-Wave User Manual 09-1109A-V.
  * Amin  _et al._ [2008] M. H. S. Amin, P. J. Love, and C. J. S. Truncik, Thermally assisted adiabatic quantum computation, Phys. Rev. Lett. 100, 060503 (2008).
  * Bravyi  _et al._ [2011] S. Bravyi, D. P. DiVincenzo, and D. Loss, Schrieffer-wolff transformation for quantum many-body systems, Annals of Physics 326, 2793 (2011).
  * Bernien  _et al._ [2017] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, _et al._ , Probing many-body dynamics on a 51-atom quantum simulator, Nature 551, 579 (2017).
  * Polkovnikov  _et al._ [2011] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Colloquium: Nonequilibrium dynamics of closed interacting quantum systems, Rev. Mod. Phys. 83, 863 (2011).
  * Kibble [1976] T. W. Kibble, Topology of cosmic domains and strings, Journal of Physics A: Mathematical and General 9, 1387 (1976).
  * Dziarmaga [2005] J. Dziarmaga, Dynamics of a quantum phase transition: Exact solution of the quantum Ising model, Phys. Rev. Lett. 95, 245701 (2005).
  * Langer [1969] J. S. Langer, Statistical theory of the decay of metastable states, Annals of Physics 54, 258 (1969).
  * Affleck [1981] I. Affleck, Quantum-statistical metastability, Phys. Rev. Lett. 46, 388 (1981).
  * Caldeira and Leggett [1981] A. O. Caldeira and A. J. Leggett, Influence of dissipation on quantum tunneling in macroscopic systems, Phys. Rev. Lett. 46, 211 (1981).
  * Leggett [1984] A. J. Leggett, Quantum tunneling in the presence of an arbitrary linear dissipation mechanism, Phys. Rev. B 30, 1208 (1984).
  * Leggett  _et al._ [1987] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, Dynamics of the dissipative two-state system, Rev. Mod. Phys. 59, 1 (1987).
  * Hänggi  _et al._ [1990] P. Hänggi, P. Talkner, and M. Borkovec, Reaction-rate theory: fifty years after Kramers, Rev. Mod. Phys. 62, 251 (1990).
  * Birnkammer  _et al._ [2022] S. Birnkammer, A. Bastianello, and M. Knap, Prethermalization in one-dimensional quantum many-body systems with confinement, Nat. Commun. 13, 7663 (2022).
  * Fendley  _et al._ [2004] P. Fendley, K. Sengupta, and S. Sachdev, Competing density-wave orders in a one-dimensional hard-boson model, Phys. Rev. B 69, 075106 (2004).
  * Lesanovsky and Katsura [2012] I. Lesanovsky and H. Katsura, Interacting Fibonacci anyons in a Rydberg gas, Phys. Rev. A 86, 041601(R) (2012).
  * Turner  _et al._ [2018] C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn, and Z. Papić, Weak ergodicity breaking from quantum many-body scars, Nat. Phys. 14, 745 (2018).
  * Su  _et al._ [2023] G.-X. Su, H. Sun, A. Hudomal, J.-Y. Desaules, Z.-Y. Zhou, B. Yang, J. C. Halimeh, Z.-S. Yuan, Z. Papić, and J.-W. Pan, Observation of many-body scarring in a Bose-Hubbard quantum simulator, Phys. Rev. Res. 5, 023010 (2023).
  * Balducci  _et al._ [2022] F. Balducci, A. Gambassi, A. Lerose, A. Scardicchio, and C. Vanoni, Localization and melting of interfaces in the two-dimensional quantum ising model, Phys. Rev. Lett. 129, 120601 (2022).
  * Hart and Nandkishore [2022] O. Hart and R. Nandkishore, Hilbert space shattering and dynamical freezing in the quantum ising model, Phys. Rev. B 106, 214426 (2022).
  * Karp [1972] R. M. Karp, Reducibility among combinatorial problems, in  _Complexity of Computer Computations_, The IBM Research Symposia Series, edited by R. E. Miller, J. W. Thatcher, and J. D. Bohlinger (Springer US, Boston, MA, 1972) pp. 85–103.
  * Breuer and Petruccione [2007] H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2007).
  * Schollwöck [2011] U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Annals of Physics 326, 96 (2011).
  * Vidal [2003] G. Vidal, Efficient classical simulation of slightly entangled quantum computations, Phys. Rev. Lett. 91, 147902 (2003).
  * Paeckel  _et al._ [2019] S. Paeckel, T. Köhler, A. Swoboda, S. R. Manmana, U. Schollwöck, and C. Hubig, Time-evolution methods for matrix-product states, Annals of Physics 411, 167998 (2019).