# Dynamics of a Nonequilibrium Discontinuous Quantum Phase Transition
in a Spinor Bose-Einstein Condensate

Matthew T. Wheeler  Physics, Faculty of Science, University of East Anglia,
Norwich Research Park, Norwich, NR4 7TJ, UK  Centre for Photonics and Quantum
Science, University of East Anglia, Norwich Research Park, Norwich, NR4 7TJ,
UK  Hayder Salman  School of Mathematics, University of East Anglia, Norwich
Research Park, Norwich, NR4 7TJ, UK  Centre for Photonics and Quantum Science,
University of East Anglia, Norwich Research Park, Norwich, NR4 7TJ, UK  Magnus
O. Borgh  Physics, Faculty of Science, University of East Anglia, Norwich
Research Park, Norwich, NR4 7TJ, UK  Centre for Photonics and Quantum Science,
University of East Anglia, Norwich Research Park, Norwich, NR4 7TJ, UK

###### Abstract

Symmetry-breaking quantum phase transitions lead to the production of
topological defects or domain walls in a wide range of physical systems. In
second-order transitions, these exhibit universal scaling laws described by
the Kibble-Zurek mechanism, but for first-order transitions a similarly
universal approach is still lacking. Here we propose a spinor Bose-Einstein
condensate as a testbed system where critical scaling behavior in a first-
order quantum phase transition can be understood from generic properties. We
generalize the Kibble-Zurek mechanism to determine the critical exponents for:
(1) the onset of the decay of the metastable state on short times scales, and
(2) the number of resulting phase-separated ferromagnetic domains at longer
times, as a one-dimensional spin-1 condensate is ramped across a first-order
quantum phase transition. The predictions are in excellent agreement with
mean-field numerical simulations and provide a paradigm for studying the decay
of metastable states in experimentally accessible systems.

##  I Introduction

Classical and quantum nonequilibrium phase transitions arise in many areas of
physics, ranging from cosmology [1, 2], to condensed matter [3, 4, 5, 6, 7,
8], and to ultracold atomic gases [9, 10, 11, 12]. For a second-order
(continuous) phase transition, a correlation length and time scale can be
identified that characterise the coherence and dynamical response of the
system. As the critical point is approached, both of these exhibit power-law
divergences described by critical exponents [13, 14]. In non-equilibrium phase
transitions, this implies that close to the critical point, the system is no
longer able to adiabatically follow the ground state [15, 16]. Causally
disconnected regions then choose the new broken symmetry state independently.
This results in the formation of topological defects or domain boundaries at a
density related to the quench rate of the control parameter.

The Kibble-Zurek mechanism (KZM) provides a theoretical framework that can
predict the density of these defects or domain walls for a finite quench rate
from universal properties of the continuous phase transition. First introduced
by Kibble in the context of cosmology as a mechanism for the formation of
cosmic strings in the early universe [17, 1], it was subsequently extended by
Zurek to condensed matter systems [18, 19, 20]. It has since been successfully
verified in many settings, including thermally driven transitions [21, 22, 23,
24, 25] and quantum phase transitions (QPTs) [26, 27, 28], and has been
demonstrated to apply to quantum-annealing implementations of quantum
computation [29, 30].

Recently, there has been increasing interest in studying first order QPTs [31,
32, 33, 34, 35], where metastability plays a crucial role, including in cold-
atom systems [36, 37, 38]. A classical example of such metastability is the
transition of supercooled water which remains liquid below the freezing point.
For first-order quantum phase transitions, such ‘supercooling’ like behaviour
can lead to a zero temperature false vacuum. This state plays an important
role in particle physics and cosmology [39, 40, 41], but an understanding of
how the metastable state decays is hampered by the lack of a general
theoretical approach dealing with first-order quantum phase transitions.

Here we propose an experimentally accessible nonequilibrium first-order QPT
where the decay of such metastable states can be understood through the KZM.
In particular, we study the persistence of the metastable state following a
finite quench of the quadratic Zeeman shift in a spin-1 Bose-Einstein
condensate (BEC) with ferromagnetic (FM) interactions. We demonstrate that the
onset of decay of the metastable state representing the false vacuum agrees
with the critical scaling law predicted by a generalisation of the KZM to our
first order QPT.

A key feature of this phase transition in a ferromagnetic spinor BEC is the
formation of phase-separated domains at long times past the transition point.
Therefore, aside from the short time behaviour characterising the decay of the
metastable state, we also apply the KZM to determine the scaling of the number
of phase separated domains at latter times. In contrast to the short time
behaviour, the KZM accurately predicts the scaling of the number of domains
for intermediate quench rates whereas at very slow quenches, a scaling in
agreement with Landau-Zener [42, 27] predictions emerges.

Focusing on an ultracold atomic system has the advantage that the QPT is
easier to control for isolated systems. Atomic BECs in particular are pristine
systems and offer a highly controllable platform where the strength of inter-
atomic interactions and the confining trapping potentials can be tuned.
Consequently, they are already popular example systems for phase-transition
experiments [9, 43, 44, 45], as well as nonequilibrium physics, even in low
dimensions, ranging from relaxation dynamics [46, 47] to quantum quenches [48,
49, 50, 51, 52, 53, 54, 55]. Unlike scalar BECs, the spin degrees of freedom
are not frozen out in spinor BECs. These additional degrees of freedom give
rise to a non-trivial phase diagram even at zero temperature [56, 57, 58, 59,
60, 61] and a correspondingly rich array of topological defects and textures
[56, 57, 62, 48, 63, 64, 65, 66]. For these reasons, studying nonequilibrium
QPTs with spinor BECs has attracted much attention [48, 67, 68, 69, 70, 71,
72, 73, 74, 75, 76, 52].

In contrast to previous studies on QPTs, the first-order QPT between the
broken-axisymmetry (BA) phase and a ferromagnetic (FM) phase in a spin-1 BEC
with FM interactions corresponds to a discontinuous quantum critical point
(DQCP) [77], the quantum analogue of the classical discontinuous critical
point [78, 79]. As it does not meet the general criteria of applicability of
the standard KZ theory, the KZM has not hitherto been studied in this context.
We consider, in particular, a one-dimensional (1D) spin-1 BEC with FM
interactions in a ring-trap geometry. By quenching the quadratic Zeeman shift,
the system can transition from the BA phase into a phase-separated FM phase
where domains of atoms with opposite condensate-spin projection form [56].
Unlike the single-mode scenario in an antiferromagnetic condensate [80], where
there is no domain formation, here phase separation is a consequence of the FM
interactions under conservation of longitudinal magnetisation in a spatially
extended BEC. Its 1D nature, however, has the advantage that once the domains
form, they are frozen in and cannot undergo any coarsening dynamics [81, 82].

##  II Results and Discussion

As our example system, we consider a spin-1 BEC described by the mean-field
condensate spinor wave function $\Psi=(\psi_{1},\psi_{0},\psi_{-1})^{T}$. The
Hamiltonian density then reads [56]

$\displaystyle
H=H_{0}+\frac{c_{0}}{2}n^{2}+\frac{c_{1}}{2}n^{2}|\langle\hat{\mathbf{F}}\rangle|^{2}-pn\langle\hat{F_{z}}\rangle+qn\langle\hat{F}_{z}^{2}\rangle,$
(1)

where $H_{0}=(\hbar^{2}/2M)|\nabla\Psi|^{2}+nV(z)$ for atomic mass $M$ and
external trapping potential $V(z)$. Here, $n=\sum_{m}\psi_{m}^{*}\psi_{m}$ is
the total atomic density. The condensate spin operator
$\hat{\mathbf{F}}\equiv(\hat{F}_{x},\hat{F}_{y},\hat{F}_{z})$ is the vector of
spin-1 Pauli-type matrices such that
$\langle\hat{F}_{\mu}\rangle=\tfrac{1}{n}\sum_{mm^{\prime}}\psi_{m}^{*}(\hat{F}_{\mu})_{mm^{\prime}}\psi_{m^{\prime}}$
for $\mu=x,y,z$. Conservation of angular momentum in $s$-wave scattering means
that the longitudinal magnetisation $M_{z}=\int\langle\hat{F}_{z}\rangle\,dz$
is conserved on experimental time scales. The spin-independent and -dependent
interaction strengths arise from the $s$-wave scattering lengths
$a_{\mathcal{F}}$ in the spin-$\mathcal{F}$ channels of colliding spin-1 atoms
as $c_{0}=4\pi\hbar^{2}(a_{0}+2a_{2})/3M$ and
$c_{1}=4\pi\hbar^{2}(a_{2}-a_{0})/3M$, respectively. Linear and quadratic
Zeeman shifts of strengths $p$ and $q$, respectively, may arise from an
applied magnetic field along the $z$-direction, or in the latter case be
induced by an AC Stark shift [83, 84], which gives precise control over both
strength and sign. Due to conservation of $M_{z}$, a uniform linear Zeeman
shift only causes precession of the spin, under which the Hamiltonian is
invariant. We therefore only consider effects from the quadratic Zeeman shift.

We consider atoms with $c_{1}<0$, such as 87Rb [85], which provides an
interesting phase diagram arising from the competition between the third and
last term in Eq. (1). A three-component BA phase with zero longitudinal
magnetisation occurs for $0<Q=q/(|c_{1}|n_{0})<2$ [56] where $n_{0}$ is the
background density in a uniform system:

$\displaystyle\psi_{\pm
1}=\frac{\sqrt{2n_{0}}}{4}\sqrt{2-Q},\qquad\psi_{0}=\frac{\sqrt{n_{0}}}{2}\sqrt{2+Q}.$
(2)

In addition, a FM state occurs for $Q<0$: $\Psi=(\sqrt{n_{0}},0,0)^{T}$ or
$\Psi=(0,0,\sqrt{n_{0}})^{T}$. Since $M_{z}$ is conserved, a BA initial
condition with $M_{z}=0$ implies that the FM phase results in the formation of
phase-separated domains with opposite spin projection as $Q$ is ramped across
the phase transition. The associated instability that leads to the emergence
of phase separated domains when $c_{1}<0$ is not captured in the single-mode
approximation [60, 61].

Moreover, the transition point between the BA and FM phases at
$(q_{c},p_{c})=(0,0)$ is a DQCP. Such a discontinuous phase transition
satisfies five conditions [77], which permit scaling arguments to be applied.
A critical point at $q=q_{c}$ and $p=p_{c}$, where $p$ functions as a
symmetry-breaking field, is a DQCP of this kind if (1) the energy density
$\epsilon(q,p)$ across the transition is continuous,
$\epsilon(q_{c}^{+},p_{c})-\epsilon(q_{c}^{-},p_{c})=0$, but (2) its
derivative is discontinuous, $\partial\epsilon(q_{c}^{+},p_{c})/\partial
q-\partial\epsilon(q_{c}^{-},p_{c})/\partial q\neq 0$. Here, $q_{c}^{+}$ and
$q_{c}^{-}$ correspond to approaching $q_{c}$ from above or below,
respectively. Further, (3) the order parameter
$m=-\partial\epsilon(q,p)/\partial p$ must exhibit a discontinuous jump with
respect to $q$ as the critical point is crossed,
$|m(q_{c}^{-},p_{c})|>|m(q_{c}^{+},p_{c})|=0$, and additionally, (4) the order
parameter is also discontinuous with respect to $p$:
$|m(q_{c},p_{c}^{\pm})|>0$. Finally, we require (5) that the derivative of the
energy be bounded as the critical point is approached:
$|\partial\epsilon(q_{c}^{\pm},p_{c})/\partial q|<\infty$.

The energy densities per particle for the BA and FM states are, respectively
[56],

$\displaystyle\epsilon_{\mathrm{BA}}$
$\displaystyle=\frac{(-p^{2}+q^{2}+2qc_{1}n_{0})^{2}}{8c_{1}n_{0}q^{2}}+\frac{1}{2}c_{0}n_{0},$
(3) $\displaystyle\epsilon_{\mathrm{FM}_{\pm}}$ $\displaystyle=\mp
p+q+\frac{1}{2}n_{0}(c_{0}+c_{1}),$ (4)

where the subscript $+(-)$ denotes the FM phase with spin pointing up (down).
These energies are continuous at the critical point $(q_{c},p_{c})=(0,0)$. The
derivatives, however, are discontinuous, but remain bounded. Meanwhile, the
relevant order parameter for the BA and FM states is
$m_{\mathrm{BA}}=p(p^{2}-q^{2}-2qc_{1}n_{0})/8c_{1}n_{0}q$ and
$m_{\mathrm{FM}_{\pm}}=\pm 1$, respectively, which is precisely the local
magnetisation $F_{z}=|\psi_{1}|^{2}-|\psi_{-1}|^{2}$ in both phases [56]. For
$p_{c}=0$ this order parameter is zero in the BA phase and becomes non-zero in
the FM phase. A non-zero $p$, however, causes the order parameter to be non-
zero in both phases. This means that the BA to FM transition satisfies all the
conditions for a DQCP.

Next, we recall the key arguments in the KZM. A continuous second-order phase
transition can be characterised by the divergence of an instantaneous
correlation length $\xi\sim|q(t)-q_{c}|^{-\nu}$ and an instantaneous
relaxation time $\tau\sim[\xi(t)]^{z}$, where $\nu$ and $z$ are the
correlation-length and dynamical critical exponents, respectively. In the
standard QPT scenario, the relaxation time is set by the inverse of the energy
gap $\Delta$ between the ground state and the first excited state of a gapped
mode [42, 67]: $\tau\simeq\Delta^{-1}$. A system initially prepared in the
ground state follows this state adiabatically as long as the relaxation time
remains small. However, as the critical point is approached, the energy gap
vanishes and the relaxation time diverges, which leads to the breakdown of the
adiabatic regime. This divides the dynamics into three stages: adiabatic,
frozen, and adiabatic again as the system crosses the critical point. Assuming
a quench of the form $|q(t)-q_{c}|\sim|t/\tau_{Q}|$, where $\tau_{Q}^{-1}$ is
the quench rate, the freezing time is estimated to be
$|\,\tilde{t}\,|=\tau(\tilde{t}\,)$. This leads to a scaling for the freezing
time given by

$\displaystyle|\,\tilde{t}\,|\sim\tau_{Q}^{z\nu/(z\nu+1)}.$ (5)

It follows from this and the scaling of the correlation length that at the
freezing time, $\tilde{\xi}\sim\xi(\,\tilde{t}\,)\sim\tau_{Q}^{\nu/(z\nu+1)}$.
This provides an estimate for the density of defects or domains:
$N_{d}\sim\tilde{\xi}^{-d}\sim\tau_{Q}^{-d\nu/(z\nu+1)}$.

Of the three Bogoliubov modes, $E_{\mathbf{k},{f_{z}}}$ and
$E_{\mathbf{k},{\pm}}$, for the BA phase [86] (see also Appendix), only
$E_{\mathbf{k},{+}}$ is gapped in the long-wavelength limit. This mode
determines the scaling associated with the KZM for the transition from the
polar to the BA phase [68]. By contrast, the relevant mode for the BA to FM
transition is

$\displaystyle E_{k,f_{z}}(k)=\sqrt{\epsilon_{k}(\epsilon_{k}+q)},$ (6)

where $\epsilon_{k}=\hbar^{2}k^{2}/2M$ is the kinetic energy. This spectrum is
gapless in the long-wavelength limit. The imaginary part of $E_{k,f_{z}}$
together with that of $E_{\mathbf{k},{+}}$ is shown in Fig. 1 and clearly
illustrates that an instability can occur at $Q=0$ only for modes with $k\neq
0$. These unstable modes are responsible for the formation of phase separated
domains in the FM phase.

Figure 1: Imaginary parts of $E_{k,f_{z}}$ and $E_{k,+}$. $E_{k,f_{z}}$ is
unstable (positive imaginary part) for $Q<0$, corresponding to the BA to FM
QDCP, while $E_{k,+}$ instead becomes unstable at the BA to polar second-order
transition at $Q=2$.

To derive a KZ scaling when the relevant mode is gapless, consider the more
general spectrum,
$E_{k}^{2}\sim|q(t)-q_{c}|^{\alpha}\epsilon_{k}^{\eta}+\epsilon_{k}^{2z}$, of
which Eq. (6) is a special case. To find scaling solutions consistent with the
KZM where $E_{k}\sim|q(t)-q_{c}|^{z}$, we make the ansatz
$k\sim|q(t)-q_{c}|^{\nu}$ (corresponding to $k\sim\xi^{-1}$) and equate the
two terms in the dispersion relation to derive the condition
$\alpha=\nu(2z-\eta)$. For our gapless mode, the adiabatic-impulse
approximation states that the impulse regime begins when
$E_{k}^{2}=\dot{E}_{k}$. This yields

$\displaystyle|\,\tilde{t}\,|\sim\tau_{Q}^{\alpha/(2+\alpha)}\tilde{k}^{-\eta/(2+\alpha)}\sim\tau_{Q}^{\nu
z/(1+\nu z)},$ (7)

for the freezing time upon using the above scaling assumed for $k$. This
immediately implies the characteristic momentum scale
$\tilde{k}\sim\tau_{Q}^{-\nu/(z\nu+1)}$ from which we extract the defect
density where $d$ is the dimensionality of the system (see also Appendix).
Hence, the scaling relations of the KZM for a gapped mode is recovered for a
gapless spectrum.

Specifically for our 1D system ($d=1$ and $q_{c}=0$), Eq. (6) implies
$\alpha=1$, $\eta=2$ and $z=2$. This is equivalent to setting $z=2$ and
$\nu=1/2$ corresponding to a defect-density scaling

$\displaystyle N_{d}\sim\tau_{Q}^{-1/4}.$ (8)

Therefore, despite originating in the same model Hamiltonian, this scaling is
clearly different from that reported for the KZM in continuous phase
transitions through a QCP in spinor BECs [68, 70, 71, 87, 73]. Our results
indicate a new scaling regime that is associated with the DQCP of this system.

To check our prediction, we numerically evolve the time-dependent spin-1
Gross-Pitaevskii equations (GPEs) obtained from Eq. (1) using a symplectic
algorithm [88] (see Appendix).

Figure 2: Component densities in a $64\xi_{s}$ subregion of the $\psi_{1}$
(top), $\psi_{0}$ (middle) and $\psi_{-1}$ (bottom) components for a typical
simulation with $\tau_{Q}=1000$.

Typical results for $\tau_{Q}=1000$ are shown in Fig. 2. We see the clear
formation of FM domains after crossing the critical point at $t=0$.

Figure 3: Normalised atom number for the $\psi_{0}$ component using the
analytical prediction (black line) in Eq. (2). Also plotted are numerically
calculated values for $\tau_{Q}=500,1000,5000$. Inset: Normalised atom number
for the $\psi_{1}$ component using the same analytical and numerical data as
the main figure. We note that $N_{1}/L$ approaches 0.5 since the longitudinal
magnetisation $M_{z}=0$ is conserved. The critical point is marked with a
black dashed line.

Figure 3 shows the normalised atom number for the $\psi_{0}$ component as the
system is quenched for various $\tau_{Q}$. Initially, the system tracks the
BA-phase ground state [80], Eq. (2). After passing the critical point, the
system is no longer able to adiabatically track the true ground state. Rather,
it evolves in a metastable BA state, even for $t/\tau_{Q}>0$, until it emerges
from the impulse regime at a time clearly dependent on the quench rate. At
this point, the metastable state decays with an associated abrupt drop in the
density of the $\psi_{0}$ component, signalling a discontinuous phase
transition to the FM phase. This coincides with an increase in the $\psi_{\pm
1}$ components, where the FM domains start to form (Fig. 3 inset).

To determine the freezing time as well as the short time scaling behaviour, we
introduce $\hat{a}_{k,\pm 1}$, defined as the Fourier transforms of the
$\psi_{\pm 1}$ components, respectively. Since the transition to the FM phase
results in phase-separated domains, driven by an instability associated with a
Bogoliubov mode related to
$\hat{a}_{k,{f_{z}}}=(\hat{a}_{k,1}-\hat{a}_{k,-1})/\sqrt{2}$ (see Appendix),
we extract the critical time $\tilde{t}$ such that
$|\hat{a}_{k,{f_{z}}}(\tilde{t}\,)|=0.01\times\mathrm{max}\\{|\hat{a}_{k,{f_{z}}}(t)|:t\\}$.
The critical Zeeman value is defined as $Q_{a}=|Q(\tilde{t}\,)|$. The inset in
Fig. 4 shows the typical growth of $|\hat{a}_{k,{f_{z}}}|$, which demonstrates
that it remains zero until the system passes the critical point. Thereafter,
it undergoes growth with strong oscillations. The choice of $0.01$ when
extracting $\tilde{t}$ is arbitrary, but we find qualitatively similar results
in tests with values up to $0.1$.

Figure 4 reveals a clear $Q_{a}\sim\tau_{Q}^{-1/2}$ power law for a large
range of $\tau_{Q}$. Despite the discontinuous nature of the phase transition,
scaling behaviour is still observed [89, 90, 91]. However, this scaling is
different from that observed in numerical and experimental results concerning
the continuous phase transitions in spin-1 BECs [68, 75]. The observed scaling
is consistent with the Kibble-Zurek scaling presented above for our system.
Taking $\nu=1/2$ and $z=2$, and using Eq. (5), we obtain
$\tilde{t}\sim\tau_{Q}^{1/2}$. Combining with the relation
$Q_{a}=|\tilde{t}/\tau_{Q}|$, we also recover the $Q_{a}\sim\tau_{Q}^{-1/2}$
scaling seen in our simulations.

Figure 4: a) Scaling of the critical value $Q_{a}$ versus $\tau_{Q}$. Overlaid
is the power-law scaling $\tau_{Q}^{-1/2}$ (dashed line). Inset: the quantity
$|a_{k,{f_{z}}}|$ for $\tau_{Q}=5000$ and $k=1$. b) Same but with $Q_{a}$
calculated using the deviation from the analytically predicted value of
$\psi_{0}$ (see Appendix) reproducing the same scaling. For both methods, each
point is averaged over 50 runs, and the error bars represent $\pm 1$ standard
deviation.

To reinforce our conclusions, we recover the scaling laws by directly
linearising the GP equations around the critical point. In this case, the
temporal dependence of the quadratic Zeeman terms is treated directly.
Following [68], we begin with a wave function close to the BA ground state,
$\Psi^{T}=(\psi_{1}+\delta\psi_{1}(t),\psi_{0}+\delta\psi_{0}(t),\psi_{-1}+\delta\psi_{-1}(t))\exp(-i\mu
t)$ where $\psi_{\pm 1},\psi_{0}$ are defined in Eq. (2) with $Q=Q_{0}$. Here,
$\mu=c_{0}+c_{1}(2-Q_{0})/2$ is the chemical potential, $0\leq Q_{0}\leq 2$ is
a constant, and $|\delta\psi_{m}(t)|\ll 1$. The noise terms have to satisfy
$\int\sum_{m}\delta\psi_{m}+\delta\psi_{m}^{*}\,\differential z=0$ to ensure
the proper normalisation of the wave function and
$\int(\delta\psi_{1}+\delta\psi_{1}^{*}+\delta\psi_{-1}+\delta\psi_{-1}^{*})\differential
z=0$ to enforce conservation of magnetisation.

Linearizing the spin-1 GPEs about the state corresponding to $Q=Q_{0}$ (see
Appendix), we obtain

$\displaystyle\\!\\!\\!i\hbar\derivative{t}G_{y}=\left[-\frac{\hbar^{2}}{2M}\derivative[2]{z}-c_{1}n_{0}\left(Q-\frac{Q_{0}}{2}\right)\right]G_{y}-\frac{c_{1}n_{0}Q}{2}G_{y}^{*},$
(9)

where $G_{y}=\delta\psi_{1}-\delta\psi_{-1}$. Next, we transform to momentum
space and split $G_{y}$ into real and imaginary parts, where
$a_{y}=\int\mathrm{Re}(G_{y})e^{ikz}dz$ and
$b_{y}=\int\mathrm{Im}(G_{y})e^{ikz}dz$. We then solve for $Q=-t/\tau_{Q}$, by
deriving the equation for $d^{2}a_{y}/dt^{2}$ across a DQCP. Rescaling time as
$t\rightarrow t\lambda$ with $\lambda=\sqrt{\tau_{s}\tau_{Q}}$, we arrive at

$\displaystyle\derivative[2]{a_{y}}{t}=\frac{-1}{\left(2\kappa^{2}-t\right)}\derivative{a_{y}}{t}-\frac{1}{4}\left(\kappa^{4}-2\kappa^{2}t+\frac{3t^{2}}{4}\right)a_{y},$
(10)

where $\kappa^{2}=\xi_{s}^{2}k^{2}\sqrt{\tau_{Q}/\tau_{s}}$. This scaling
ensures that the last term is independent of $\tau_{Q}$. The remaining
dependence on $\tau_{Q}$ is eliminated if we require that $\kappa$ is
constant, which implies $k\sim\tau_{Q}^{-1/4}$. Only under these conditions
can we expect scaling solutions, and they are also consistent with the scaling
of the correlation length and the dynamical exponents derived earlier based on
the KZM.

The number of domains formed in the transition is a measurable quantity, which
has been investigated in other works on the KZM [81, 73]. We numerically count
the number of density peaks at the end of the simulation (Fig. 2).

Figure 5: The number of ferromagnetic domains as a function of the quench
time, $\tau_{Q}$. Overlaid is the power-law scaling $\tau_{Q}^{-1/4}$ (dashed
line). Each point is averaged over 50 runs, and the error bars give $\pm 1$
standard deviation. Inset: number of ferromagnetic domains in a quench that
spans two phase transitions for the same quench times as the main figure.

Fig. 5 shows the total number of FM domains for a broad range of $\tau_{Q}$.
For sufficiently fast quenches ($\tau_{Q}<1000$), a clear power-law scaling
$N_{d}\sim\tau_{Q}^{-1/4}$ emerges, which agrees well with Eq. (8) as well as
the scaling obtained with Eq. (10). As with $\tilde{t}$, the scaling for the
DQCP is again different from that in analogous transitions across a continuous
critical point [81].

Unlike $Q_{a}$, for slow quenches ($\tau_{Q}>1000$), there is a clear
deviation from the predicted KZ scaling. Such differences in the scaling of
observables measured at much later times from the critical point have also
been observed in [21, 73]. In our context, this deviation can be attributed to
a cross-over into a different regime that follows predictions of a Landau-
Zener model [42, 27, 92, 93].

Finally, we tested the robustness of the scaling by considering a quench that
crosses two phase transitions. Starting from the polar state with $Q=2.5$, we
simulate a quench through the BA and then the FM states. Results are
qualitatively similar (Fig. 5 inset), with a $\tau_{Q}^{-1/4}$ scaling again
emerging for the same range of $\tau_{Q}$, and with similar deviation for slow
quenches.

In conclusion, we have shown that the KZM can be generalised to this
discontinuous phase transition, leading to scaling laws that differ from those
observed for phase transitions across continuous quantum critical points for
the same spin-1 BEC model. We find excellent agreement with numerical
simulation for both the short-time growth of the unstable excitations and the
subsequent number of domains formed on longer time scales. Our results hold
for experimentally accessible parameter regimes allowing these extensions of
the KZM to be realized in current experiments on spinor BECs, which therefore
emerge as prime candidates for testbed systems for investigating critical
scaling in first-order quantum phase transitions, including as laboratory
emulators for understanding false-vacuum decay [36, 37, 38, 94].

###### Acknowledgements.

The numerical simulations were carried out on the High Performance Computing
Cluster supported by the Research and Specialist Computing Support service at
the University of East Anglia. MOB acknowledges support from EPSRC under grant
number EP/V03832X/1.

*

## Appendix A

###  A.1 Numerical Simulation.

We measure length and time in units of the spin healing length
$\xi_{s}=\hbar/\sqrt{2M|c_{1}|n_{0}}$ and the spin time
$\tau_{s}=\hbar/2|c_{1}|n_{0}$, respectively. Our simulations are performed on
a 1D grid of $N_{x}=16384$ points with a spacing of $\Delta_{x}=0.125\xi_{s}$,
considering a ring-shaped geometry by assuming periodic boundary conditions
and $V(z)=0$. We start from Eq. (2), adding small noise terms,
$\delta\psi_{m}$, to each component, where $|\delta\psi_{m}|\ll 1$. The real
and imaginary parts of $\delta\psi_{m}$ are drawn from the probability
distribution $p(z)=\exp(-z^{2}/2\sigma^{2})/(\sqrt{2\pi}\sigma)$, with
$\sigma=10^{-4}$ to remain close to the BA ground state. We vary the quadratic
Zeeman shift as $Q(t)=-t/\tau_{Q}$ for a range of quench times $\tau_{Q}$,
starting at $Q=1$ and ending the simulation at $Q=-2.5$.

###  A.2 Bogoliubov modes for the broken-axisymmetry phase.

The broken-axisymmetry (BA) phase of a spin-1 BEC exhibits three Bogoliubov
modes [86]. Here we rederive each mode explicitly from the relevant Bogoliubov
transformations and explain why $E_{\mathbf{k},f_{z}}$ is the relevant mode
for the BA-to-ferromagnetic (FM) transition.

The broken-axisymmetry phase can be parameterised as

$\displaystyle\zeta^{\text{BA}}=\left(\frac{\sin\theta}{\sqrt{2}},\cos\theta,\frac{\sin\theta}{\sqrt{2}}\right),$
(11)

where $\sin\theta=\sqrt{1/2+q/(4nc_{1})}$. The fluctuation operators for this
state are then defined as [86]:

$\displaystyle\hat{a}_{\mathbf{k},d}$
$\displaystyle=\frac{\sin\theta}{\sqrt{2}}(\hat{a}_{\mathbf{k},1}+\hat{a}_{\mathbf{k},-1})+\cos\theta\hat{a}_{\mathbf{k},0},$
(12) $\displaystyle\hat{a}_{\mathbf{k},f_{z}}$
$\displaystyle=\frac{1}{\sqrt{2}}(\hat{a}_{\mathbf{k},1}-\hat{a}_{\mathbf{k},-1}),$
(13) $\displaystyle\hat{a}_{\mathbf{k},\theta}$
$\displaystyle=\frac{\cos\theta}{\sqrt{2}}(\hat{a}_{\mathbf{k},1}+\hat{a}_{\mathbf{k},-1})-\sin\theta\hat{a}_{\mathbf{k},0},$
(14)

where on the right-hand side $\hat{a}_{\mathbf{k},m}$ is the annihilation
operator for a spin-$1$ boson in magnetic level $m$ (for $m=-1,0,+1$),
determined by expanding the wave-function field operator as

$\displaystyle\hat{\psi}_{m}(\mathbf{x})=\frac{1}{\sqrt{V}}\sum_{\mathbf{k}}\hat{a}_{\mathbf{k},m}e^{i\mathbf{k}\cdot\mathbf{x}},$
(15)

where $V$ is the volume of the system.

The sub-Hamiltonian for the spin fluctuation mode $\hat{a}_{\mathbf{k},f_{z}}$
can be diagonalized using the transformation

$\displaystyle\hat{b}_{\mathbf{k},f_{z}}=\sqrt{\frac{\epsilon_{\mathbf{k}}+q/2+E_{\mathbf{k},f_{z}}}{2E_{\mathbf{k},f_{z}}}}\hat{a}_{\mathbf{k},f_{z}}+\sqrt{\frac{\epsilon_{\mathbf{k}}+q/2-E_{\mathbf{k},f_{z}}}{2E_{\mathbf{k},f_{z}}}}\hat{a}_{-\mathbf{k},f_{z}}^{\dagger},$
(16)

where $\epsilon_{\mathbf{k}}=\hbar^{2}|\mathbf{k}|^{2}/2M$ is the kinetic
energy and the Bogoliubov spectrum is given by

$\displaystyle
E_{\mathbf{k},{f_{z}}}=\sqrt{\epsilon_{\mathbf{k}}(\epsilon_{\mathbf{k}}+q)}.$
(17)

The sub-Hamiltonians for the density fluctuation mode $\hat{a}_{\mathbf{k},d}$
and the $\theta$ mode $\hat{a}_{\mathbf{k},\theta}$ can be similarly
diagonalized using operators $\hat{b}_{\mathbf{k},+}$ and
$\hat{b}_{\mathbf{k},+}$, which yields the remaining two Bogoliubov modes
[86]:

$\displaystyle
E_{\mathbf{k},\pm}=\sqrt{\epsilon_{\mathbf{k}}^{2}+n(c_{0}-c_{1})\epsilon_{\mathbf{k}}+2{(nc_{1})}^{2}(1-\tilde{q}^{2})\pm
E_{1}(\mathbf{k})},$ (18)

where $\tilde{q}=-q/2c_{1}n$ and

$\displaystyle
E_{1}(\mathbf{k})=\sqrt{{\left[n{(c_{0}+3c_{1})}\epsilon_{\mathbf{k}}+2{(c_{1}n)}^{2}(1-\tilde{q}^{2})\right]}^{2}-4c_{1}(c_{0}+2c_{1}){(n\tilde{q}\epsilon_{\mathbf{k}})}^{2}}.$
(19)

The final, diagonalized Hamiltonian then reads

$\displaystyle\begin{split}\hat{H}^{\text{BA}}=\,&E_{0}^{\text{BA}}+\sum_{\mathbf{k}\neq
0}[E_{\mathbf{k},f_{z}}\hat{b}^{\dagger}_{\mathbf{k},f_{z}}\hat{b}_{\mathbf{k},f_{z}}\\\
&+E_{\mathbf{k},-}\hat{b}^{\dagger}_{\mathbf{k},-}\hat{b}_{\mathbf{k},-}+E_{\mathbf{k},+}\hat{b}^{\dagger}_{\mathbf{k},+}\hat{b}_{\mathbf{k},+}],\end{split}$
(20)

where $E_{0}^{\text{BA}}$ is the ground state energy for the BA phase, which
is explicitly derived in Ref. [86].

We now consider our 1D system. In the long-wavelength limit, $k\rightarrow 0$,
the only non-zero (gapped) mode is
$E_{k,+}=\sqrt{4{(c_{1}n)}^{2}(1-\tilde{q}^{2})}$ which has the form
$E_{k,+}\sim\sqrt{{q_{c}^{\prime}}^{2}-q^{2}}$ with $q_{c}^{\prime}=2c_{1}n$.
The relevant mode of the BA to FM transition is found from the imaginary parts
of the Bogoliubov energies (Fig 1). For $|Q|>2$, where $Q\equiv q/(|c_{1}|n)$,
$E_{k,+}$ has a positive imaginary part, indicating instability. The critical
point $Q=2$ ($q=q_{c}^{\prime}$) corresponds to the second-order transition
between the polar and BA phases, for which $E_{k,+}$ therefore is the
corresponding Bogoliubov energy.

However, for $Q<0$ the imaginary part of $E_{k,f_{z}}$ mode becomes non-zero
and positive, and thus unstable. This corresponds precisely to the transition
from the BA to the FM phase at $Q=0$ that we are interested in here, and
therefore $E_{k,f_{z}}$ is the relevant mode to study. Note that the
$E_{k,f_{z}}$ mode does not give rise to instability at $k=0$. Therefore,
studies focusing on this mode at $k=0$ only do not capture the phase
transition that occurs at $Q=0$ [60, 80, 61]. In contrast, the $k=0$ mode
corresponds to the most unstable mode for $E_{k,+}$, and thus it suffices to
choose this Bogoliubov energy to capture the phase transition that occurs at
$Q=2$. In practice, the $Q=-2$ transition is not realized since the
instability of $E_{k,f_{z}}$ at any $k\neq 0$ will typically arise at $Q=0$
when $Q$ is quenched from positive to negative values.

###  A.3 Scaling of Density of Defects from the Kibble-Zurek Mechanism.

Given the form of our dispersion relation together with the scaling relation
for the relaxation time provided by the KZM, which assumes $\tau\sim\xi^{z}$,
we will consider a dispersion relation of the general form [77]

$\displaystyle\omega^{2}\sim|q-q_{c}|^{\alpha}k^{\eta}+k^{2z}$ (21)

for a critical point at $q=q_{c}$. For a scaling solution to arise that is
consistent with the scaling of the relaxation time, we require $\omega\sim
k^{z}$ with $k\sim\xi^{-1}\sim|q-q_{c}|^{\nu}$. Combining these relations, Eq.
(21) implies that

$\displaystyle\alpha=\nu(2z-\eta).$ (22)

In our system $q_{c}=0$ and $q=-t/\tau_{Q}$. Therefore, Eq. (21) simplifies to

$\displaystyle\omega^{2}\sim\frac{|t|^{\alpha}}{\tau_{Q}^{\alpha}}k^{\eta}.$
(23)

Now let us consider the adiabatic-impulse approximation for a gapped spectrum,
with an energy gap given by $\Delta(t)$. In this case, $\tau\sim 1/\Delta$ and
so far from the critical point where $\Delta$ is large, the relaxation time is
small and the evolution is adiabatic, meaning that the system is able to
adjust quickly enough and can therefore track the true ground state of the
system. As the critical point is approached, $\Delta$ vanishes and the
relaxation time diverges. At some instant $\tilde{t}$, the reaction time
becomes comparable to the transition time, $\Delta/\dot{\Delta}$ and the
system is no longer able to adiabatically track the evolution. This point is
the onset of the impulse regime, where the state becomes frozen. The freezing
time is therefore often evaluated from the condition

$\displaystyle\frac{1}{\Delta(\tilde{t})}\sim\frac{\Delta(\tilde{t})}{\dot{\Delta}(\tilde{t})}.$
(24)

For a gapless spectrum, which is the case for the BA-to-FM transition, we work
with the form of the dispersion relation as given by Eq. (23) and determine
the onset of the impulse regime with the condition $\omega^{2}=\dot{\omega}$.
This yields

$\displaystyle\frac{\alpha}{2\tau_{Q}^{\alpha/2}}|\tilde{t}|^{(\alpha/2-1)}\tilde{k}^{\eta/2}\sim\frac{|\tilde{t}|^{\alpha}}{\tau_{Q}^{\alpha}}\tilde{k}^{\eta},$
(25)

which leads us directly to the scaling relation for the freezing time

$\displaystyle|\tilde{t}|\sim\tau_{Q}^{\alpha/(2+\alpha)}\tilde{k}^{-\eta/(2+\alpha)}.$
(26)

To obtain the characteristic momentum scale, we substitute Eq. (26) back into
Eq. (25), which yields

$\displaystyle\tilde{k}\sim\tau_{Q}^{\alpha/(\eta-(2+\alpha)z)}.$ (27)

###  A.4 Extracting the freezing time

In order to extract the freezing time $\tilde{t}$ of the system, an
appropriate quantity must be chosen. Since the transition to the FM phase
causes the formation of phase-separated domains, a natural choice is to
measure fluctuations in the difference of the populations of $\psi_{\pm 1}$
[86]. To do this, we first construct the Fourier transforms of the $\psi_{\pm
1}$ components as $\hat{a}_{\pm 1}(k)=\int\psi_{\pm 1}e^{-2\pi
ik\cdot\mathbf{x}}\,d\mathbf{x}$. After passing through the critical point
into the FM phase, the difference
$\hat{a}_{k,f_{z}}(k)=\hat{a}_{1}(k)-\hat{a}_{-1}(k)/\sqrt{2}$ generates
measurable fluctuations as FM domains with opposite spin start to form [see
inset of Fig. 4(a)], while before the transition it remains zero (due to the
absence of domains). To measure the freezing time, we extract the time at
which $|\hat{a}_{k,f_{z}}(k)|$ exceeds some appropriately chosen value. In our
numerical simulations, we take this value to be 1% of the maximum value of
$|\hat{a}_{k,f_{z}}(k)|$ over the entire simulation.

An alternative choice is the population of the $\psi_{0}$ component. Here,
instead of measuring the growth of a quantity, we now extract the freezing
time as the time required for the $\psi_{0}$ component to deviate from its
analytically calculated value in the (metastable) BA phase
$\psi_{0}=\sqrt{2+Q}/2$ (Fig. 3). In particular, we choose the freezing time
to be the time at which the deviation reaches 1% of the analytically predicted
value, yielding Fig. 4(b). We see that, despite using an entirely different
quantity to measure the freezing time, the resulting scaling of $Q_{a}$ is the
same as when calculated from $|\hat{a}_{k,f_{z}}(k)|$.

###  A.5 Deriving scaling near the critical point.

The spin-1 Gross-Pitaevskii equations (GPEs) are given as [56]:

$\displaystyle
i\hbar\partialderivative{\Psi}{t}=\left[-\frac{\hbar^{2}\nabla^{2}}{2M}-p\hat{F_{z}}+q\hat{F}_{z}^{2}+c_{0}n+c_{1}n\langle\hat{\mathbf{F}}\rangle\cdot\hat{\mathbf{F}}\right]\Psi.$
(28)

Recall that we start from a BA phase of the form
$\Psi^{T}=(\psi_{1}+\delta\psi_{1}(t),\psi_{0}+\delta\psi_{0}(t),\psi_{-1}+\delta\psi_{-1}(t))\exp(-i\mu
t)$. Substituting this expression into the GPEs and keeping leading order
terms in $\delta\psi_{m}$ yields the following equations for $\delta\psi_{\pm
1}$ ($p=0$)

$\begin{split}i\hbar\partialderivative{\delta\psi_{1}}{t}=&\left[-\frac{\hbar^{2}}{2M}\derivative[2]{z}+q-\mu+\frac{n_{0}(10c_{0}+6c_{1}-(c_{0}-c_{2})Q)}{8}\right]\delta\psi_{1}+\frac{n_{0}\sqrt{2(4-Q^{2})}}{8}\left[(c_{0}+3c_{1})\delta\psi_{0}+(c_{0}+c_{1})\delta\psi_{0}^{*}\right]\\\
&+\frac{n_{0}(2-Q)}{8}\left[(c_{0}-c_{1})\delta\psi_{-1}+(c_{0}+c_{1})\delta\psi_{1}^{*}\right]+\frac{n_{0}}{8}\left[(2-Q)c_{0}+(2+3Q)c_{1}\right]\delta\psi_{-1}^{*},\end{split}$
(29)
$\begin{split}i\hbar\partialderivative{\delta\psi_{-1}}{t}=&\left[-\frac{\hbar^{2}}{2M}\derivative[2]{z}+q-\mu+\frac{n_{0}(10c_{0}+6c_{1}-(c_{0}-c_{2})Q)}{8}\right]\delta\psi_{-1}+\frac{n_{0}\sqrt{2(4-Q^{2})}}{8}\left[(c_{0}+3c_{1})\delta\psi_{0}+(c_{0}+c_{1})\delta\psi_{0}^{*}\right]\\\
&+\frac{n_{0}(2-Q)}{8}\left[(c_{0}-c_{1})\delta\psi_{1}+(c_{0}+c_{1})\delta\psi_{-1}^{*}\right]+\frac{n_{0}}{8}\left[(2-Q)c_{0}+(2+3Q)c_{1}\right]\delta\psi_{1}^{*}.\end{split}$
(30)

Subtracting Eq. (30) from Eq. (29) results in the differential equation for
$G_{y}=\delta\psi_{1}-\delta\psi_{-1}$:

$\displaystyle
i\hbar\derivative{G_{y}}{t}=\left[-\frac{\hbar^{2}}{2M}\derivative[2]{z}+q-\mu+n_{0}(c_{0}+c_{1})\right]G_{y}-\frac{c_{1}n_{0}Q}{2}G_{y}^{*}.$
(31)

Additionally, to calculate the chemical potential, we take the $\psi_{0}$
component of Eq. (28) keeping lead order terms and assuming a stationary
state, which leads to $\mu=c_{0}n_{0}+c_{1}n_{0}(2-Q_{0})/2$ where $Q_{0}$ is
a constant. Substituting this expression and $q(t)=-c_{1}n_{0}Q(t)$ into Eq.
(31) yields

$\displaystyle
i\hbar\derivative{G_{y}}{t}=\left[-\frac{\hbar^{2}}{2M}\derivative[2]{z}-c_{1}n_{0}\left(Q-\frac{Q_{0}}{2}\right)\right]G_{y}-\frac{c_{1}n_{0}Q}{2}G_{y}^{*}.$
(32)

To progress, we split $G_{y}$ into real and imaginary parts and take the
Fourier transform: $a_{y}=\int\mathrm{Re}(G_{y})e^{ikz}dz$ and
$b_{y}=\int\mathrm{Im}(G_{y})e^{ikz}dz$. Substituting into Eq. (32) yields the
matrix equation

To solve the above equation, we construct the equation for
$\derivative[2]{a_{y}}{t}$ and take $Q_{0}=0$, which yields

$\displaystyle\derivative[2]{a_{y}}{t}=\frac{c_{1}n_{0}}{2\hbar\tau_{Q}}b_{y}+\left(\frac{\hbar^{2}k^{2}}{2M}-\frac{c_{1}n_{0}Q}{2\hbar}\right)\derivative{b_{y}}{t}.$
(33)

Expressions for $b_{y}$ and $db_{y}/dt$ are found from Eq. (A.5). Substituting
these in yields the following equation for $d^{2}a_{y}/dt^{2}$:

$\derivative[2]{a_{y}}{t}=\frac{1}{\tau_{Q}\left(\frac{\hbar^{2}k^{2}}{Mc_{1}n_{0}}-Q\right)}\derivative{a_{y}}{t}-\left(\frac{\hbar^{2}k^{4}}{4M^{2}}-\frac{k^{2}c_{1}n_{0}Q}{M}+\frac{3c_{1}^{2}n_{0}^{2}Q^{2}}{4\hbar^{2}}\right)a_{y}.$
(34)

To simplify the above expression, we use the spin healing length
$\xi_{s}=\hbar/\sqrt{2|c_{1}|n_{0}}$ and the spin time
$\tau_{s}=\hbar/|c_{1}|n_{0}$:

$\displaystyle\\!\\!\\!\\!\\!\\!\derivative[2]{a_{y}}{t}=\frac{-1}{\left(2\xi_{s}^{2}k^{2}\tau_{Q}-t\right)}\derivative{a_{y}}{t}-\left(\frac{\xi_{s}^{4}k^{4}}{4\tau_{s}^{2}}-\frac{\xi_{s}^{2}k^{2}t}{2\tau_{s}^{2}\tau_{Q}}+\frac{3t^{2}}{16\tau_{s}^{2}\tau_{Q}^{2}}\right)a_{y}.$
(35)

Rescaling time as $t\rightarrow t\lambda$ with
$\lambda=\sqrt{\tau_{s}\tau_{Q}}$, leads to the differential equation

$\displaystyle\derivative[2]{a_{y}}{t}=\frac{-1}{\left(2\kappa^{2}-t\right)}\derivative{a_{y}}{t}-\frac{1}{4}\left(\kappa^{4}-2\kappa^{2}t+\frac{3t^{2}}{4}\right)a_{y},$
(36)

where $\kappa^{2}=\xi_{s}^{2}k^{2}\sqrt{\tau_{Q}/\tau_{s}}$.

## References

  * Kibble [1980] T. W. Kibble, Phys. Rep. 67, 183 (1980).
  * Mazumdar and White [2019] A. Mazumdar and G. White, Rep. Prog. Phys. 82, 076901 (2019).
  * Chuang  _et al._ [1991] I. Chuang, R. Durrer, N. Turok, and B. Yurke, Science 251, 1336 (1991).
  * Hendry  _et al._ [1994] P. C. Hendry, N. S. Lawson, R. A. Lee, P. V. McClintock, and C. D. Williams, Nature 368, 315 (1994).
  * Bäuerle  _et al._ [1996] C. Bäuerle, Y. M. Bunkov, S. N. Fisher, H. Godfrin, and G. R. Pickett, Nature 382, 332 (1996).
  * Ruutu  _et al._ [1996] V. M. H. Ruutu, V. B. Eltsov, A. J. Gill, T. W. B. Kibble, M. Krusius, Y. G. Makhlin, B. Plaçais, G. E. Volovik, and W. Xu, Nature 382, 334 (1996).
  * Sondhi  _et al._ [1997] S. L. Sondhi, S. M. Girvin, J. P. Carini, and D. Shahar, Rev. Mod. Phys. 69, 315 (1997).
  * Polkovnikov  _et al._ [2011] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Rev. Mod. Phys. 83, 863 (2011).
  * Hadzibabic  _et al._ [2006] Z. Hadzibabic, P. Krüger, M. Cheneau, B. Battelier, and J. Dalibard, Nature 441, 1118 (2006).
  * Langen  _et al._ [2015] T. Langen, R. Geiger, and J. Schmiedmayer, Annu. Rev. Condens. Matter Phys. 6, 201 (2015).
  * Fletcher  _et al._ [2015] R. J. Fletcher, M. Robert-de-Saint-Vincent, J. Man, N. Navon, R. P. Smith, K. G. H. Viebahn, and Z. Hadzibabic, Phys. Rev. Lett. 114, 255302 (2015).
  * Liu  _et al._ [2018] I.-K. Liu, S. Donadello, G. Lamporesi, G. Ferrari, S.-C. Gou, F. Dalfovo, and N. P. Proukakis, Commun. Phys. 1, 24 (2018).
  * Hohenberg and Halperin [1977] P. C. Hohenberg and B. I. Halperin, Rev. Mod. Phys. 49, 435 (1977).
  * Goldenfeld [1992] N. Goldenfeld,  _Lectures on Phase Transitions and the Renormalization Group_ (CRC Press, 1992).
  * Dziarmaga [2010] J. Dziarmaga, Adv. Phys. 59, 1063 (2010).
  * Del Campo and Zurek [2014] A. Del Campo and W. H. Zurek, Int. J. Mod. Phys. A 29, 1430018 (2014).
  * Kibble [1976] T. W. Kibble, J. Phys. Math. Gen. 9, 1387 (1976).
  * Zurek [1985] W. H. Zurek, Nature 317, 505 (1985).
  * Zurek [1993] W. Zurek, Acta Phys. Pol. B 24, 1301 (1993).
  * Zurek [1996] W. H. Zurek, Phys. Rep. 276, 177 (1996).
  * Su  _et al._ [2013] S.-W. Su, S.-C. Gou, A. Bradley, O. Fialko, and J. Brand, Phys. Rev. Lett. 110, 215302 (2013).
  * Lamporesi  _et al._ [2013] G. Lamporesi, S. Donadello, S. Serafini, F. Dalfovo, and G. Ferrari, Nat. Phys. 9, 656 (2013).
  * Liu  _et al._ [2014] C.-W. Liu, A. Polkovnikov, and A. W. Sandvik, Phys. Rev. B 89, 054307 (2014).
  * Donadello  _et al._ [2016] S. Donadello, S. Serafini, T. Bienaimé, F. Dalfovo, G. Lamporesi, and G. Ferrari, Phys. Rev. A 94, 023628 (2016).
  * Beugnon and Navon [2017] J. Beugnon and N. Navon, J. Phys. B: At. Mol. Opt. Phys. 50, 022002 (2017).
  * Dziarmaga [2005] J. Dziarmaga, Phys. Rev. Lett. 95, 245701 (2005).
  * Damski [2005] B. Damski, Phys. Rev. Lett. 95, 035701 (2005).
  * Yi  _et al._ [2020] C.-R. Yi, S. Liu, R.-H. Jiao, J.-Y. Zhang, Y.-S. Zhang, and S. Chen, Phys. Rev. Lett. 125, 260603 (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, Nat. Phys. 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, Nature 617, 61 (2023).
  * Coulamy  _et al._ [2017] I. B. Coulamy, A. Saguia, and M. S. Sarandy, Phys. Rev. E 95, 022127 (2017).
  * Shimizu  _et al._ [2018] K. Shimizu, T. Hirano, J. Park, Y. Kuno, and I. Ichinose, New J. Phys. 20, 083006 (2018).
  * Pelissetto  _et al._ [2018] A. Pelissetto, D. Rossini, and E. Vicari, Phys. Rev. E 97, 052148 (2018).
  * Pelissetto  _et al._ [2020] A. Pelissetto, D. Rossini, and E. Vicari, Phys. Rev. E 102, 012143 (2020).
  * Sinha  _et al._ [2021] A. Sinha, T. Chanda, and J. Dziarmaga, Phys. Rev. B 103, L220302 (2021).
  * Billam  _et al._ [2022] T. P. Billam, K. Brown, and I. G. Moss, Phys. Rev. A 105, L041301 (2022).
  * Song  _et al._ [2022] B. Song, S. Dutta, S. Bhave, J.-C. Yu, E. Carter, N. Cooper, and U. Schneider, Nat. Phys. 18, 259 (2022).
  * [38] A. Zenesini, A. Berti, R. Cominotti, C. Rogora, I. G. Moss, T. P. Billam, I. Carusotto, G. Lamporesi, A. Recati, and G. Ferrari, Observation of false vacuum decay via bubble formation in ferromagnetic superfluids, arXiv:2305.05225 [hep-ph] .
  * Coleman [1977] S. Coleman, Phys. Rev. D 15, 2929 (1977).
  * Fialko  _et al._ [2015] O. Fialko, B. Opanchuk, A. I. Sidorov, P. D. Drummond, and J. Brand, Europhys. Lett. 110, 56001 (2015).
  * Devoto  _et al._ [2022] F. Devoto, S. Devoto, L. D. Luzio, and G. Ridolfi, J. Phys. G: Nucl. Part. Phys. 49, 103001 (2022).
  * Zurek  _et al._ [2005] W. H. Zurek, U. Dorner, and P. Zoller, Phys. Rev. Lett. 95, 105701 (2005).
  * Krüger  _et al._ [2007] P. Krüger, Z. Hadzibabic, and J. Dalibard, Phys. Rev. Lett. 99, 040402 (2007).
  * Weiler  _et al._ [2008] C. N. Weiler, T. W. Neely, D. R. Scherer, A. S. Bradley, M. J. Davis, and B. P. Anderson, Nature 455, 948 (2008).
  * Chomaz  _et al._ [2015] L. Chomaz, L. Corman, T. Bienaimé, R. Desbuquois, C. Weitenberg, S. Nascimbène, J. Beugnon, and J. Dalibard, Nat. Commun. 6, 6162 (2015).
  * Gring  _et al._ [2012] M. Gring, M. Kuhnert, T. Langen, T. Kitagawa, B. Rauer, M. Schreitl, I. Mazets, D. A. Smith, E. Demler, and J. Schmiedmayer, Science 337, 1318 (2012).
  * Reeves  _et al._ [2022] M. T. Reeves, K. Goddard-Lee, G. Gauthier, O. R. Stockdale, H. Salman, T. Edmonds, X. Yu, A. S. Bradley, M. Baker, H. Rubinsztein-Dunlop, M. J. Davis, and T. W. Neely, Phys. Rev. X 12, 011031 (2022).
  * Sadler  _et al._ [2006] L. E. Sadler, J. M. Higbie, S. R. Leslie, M. Vengalattore, and D. M. Stamper-Kurn, Nature 443, 312 (2006).
  * Barnett  _et al._ [2011] R. Barnett, A. Polkovnikov, and M. Vengalattore, Phys. Rev. A 84, 023606 (2011).
  * Navon  _et al._ [2015] N. Navon, A. L. Gaunt, R. P. Smith, and Z. Hadzibabic, Science 347, 167 (2015).
  * Symes and Blakie [2017] L. M. Symes and P. B. Blakie, Phys. Rev. E 95, 013311 (2017).
  * Kang  _et al._ [2017] S. Kang, S. W. Seo, J. H. Kim, and Y. Shin, Phys. Rev. A 95, 053638 (2017).
  * Prüfer  _et al._ [2018] M. Prüfer, P. Kunkel, H. Strobel, S. Lannig, D. Linnemann, C. M. Schmied, J. Berges, T. Gasenzer, and M. K. Oberthaler, Nature 563, 217 (2018).
  * Schmied  _et al._ [2019] C.-M. Schmied, T. Gasenzer, and P. B. Blakie, Phys. Rev. A 100, 033603 (2019).
  * Liu  _et al._ [2020] I.-K. Liu, J. Dziarmaga, S.-C. Gou, F. Dalfovo, and N. P. Proukakis, Phys. Rev. Res. 2, 033183 (2020).
  * Kawaguchi and Ueda [2012] Y. Kawaguchi and M. Ueda, Phys. Rep. 520, 253 (2012).
  * Stamper-Kurn and Ueda [2013] D. M. Stamper-Kurn and M. Ueda, Rev. Mod. Phys. 85, 1191 (2013).
  * Murata  _et al._ [2007] K. Murata, H. Saito, and M. Ueda, Phys. Rev. A 75, 013607 (2007).
  * Borgh  _et al._ [2014] M. O. Borgh, J. Lovegrove, and J. Ruostekoski, New. J. Phys. 16, 053046 (2014).
  * Matuszewski  _et al._ [2009] M. Matuszewski, T. J. Alexander, and Y. S. Kivshar, Phys. Rev. A 80, 023602 (2009).
  * Mirkhalaf  _et al._ [2021] S. S. Mirkhalaf, D. Benedicto Orenes, M. W. Mitchell, and E. Witkowska, Phys. Rev. A 103, 023317 (2021).
  * Leanhardt  _et al._ [2003] A. E. Leanhardt, Y. Shin, D. Kielpinski, D. E. Pritchard, and W. Ketterle, Phys. Rev. Lett. 90, 140403 (2003).
  * Seo  _et al._ [2015] S. W. Seo, S. Kang, W. J. Kwon, and Y. I. Shin, Phys. Rev. Lett. 115, 015301 (2015).
  * Kang  _et al._ [2019] S. Kang, S. W. Seo, H. Takeuchi, and Y. Shin, Phys. Rev. Lett. 122, 095301 (2019).
  * Weiss  _et al._ [2019] L. S. Weiss, M. O. Borgh, A. Blinova, T. Ollikainen, M. Möttönen, J. Ruostekoski, and D. S. Hall, Nat. Commun. 10, 4772 (2019).
  * Xiao  _et al._ [2022] Y. Xiao, M. O. Borgh, A. Blinova, T. Ollikainen, J. Ruostekoski, and D. S. Hall, Nat. Commun. 13, 4635 (2022).
  * Damski and Zurek [2006] B. Damski and W. H. Zurek, Phys. Rev. A 73, 063405 (2006).
  * Damski and Zurek [2007] B. Damski and W. H. Zurek, Phys. Rev. Lett. 99, 130402 (2007).
  * Lamacraft [2007] A. Lamacraft, Phys. Rev. Lett. 98, 160404 (2007).
  * Saito  _et al._ [2007a] H. Saito, Y. Kawaguchi, and M. Ueda, Phys. Rev. A 75, 013621 (2007a).
  * Saito  _et al._ [2007b] H. Saito, Y. Kawaguchi, and M. Ueda, Phys. Rev. A 76, 043613 (2007b).
  * Vengalattore  _et al._ [2008] M. Vengalattore, S. R. Leslie, J. Guzman, and D. M. Stamper-Kurn, Phys. Rev. Lett. 100, 170403 (2008).
  * Świsłocki  _et al._ [2013] T. Świsłocki, E. Witkowska, J. Dziarmaga, and M. Matuszewski, Phys. Rev. Lett. 110, 045303 (2013).
  * Witkowska  _et al._ [2013] E. Witkowska, J. Dziarmaga, T. Świsłocki, and M. Matuszewski, Phys. Rev. B 88, 054508 (2013).
  * Anquez  _et al._ [2016] M. Anquez, B. A. Robbins, H. M. Bharath, M. Boguslawski, T. M. Hoang, and M. S. Chapman, Phys. Rev. Lett. 116, 155301 (2016).
  * Williamson and Blakie [2016] L. A. Williamson and P. B. Blakie, Phys. Rev. A 94, 063615 (2016).
  * Suzuki and Dutta [2015] S. Suzuki and A. Dutta, Phys. Rev. B 92, 064419 (2015).
  * Nauenberg [1975] M. Nauenberg, J. Phys. A: Math. Gen. 8, 925 (1975).
  * Fisher and Berker [1982] M. E. Fisher and A. N. Berker, Phys. Rev. B 26, 2507 (1982).
  * Qiu  _et al._ [2020] L.-Y. Qiu, H.-Y. Liang, Y.-B. Yang, H.-X. Yang, T. Tian, Y. Xu, and L.-M. Duan, Sci. Adv. 6, eaba7292 (2020).
  * Sabbatini  _et al._ [2011] J. Sabbatini, W. H. Zurek, and M. J. Davis, Phys. Rev. Lett. 107, 230402 (2011).
  * Sabbatini  _et al._ [2012] J. Sabbatini, W. H. Zurek, and M. J. Davis, New J. Phys. 14, 095030 (2012).
  * Gerbier  _et al._ [2006] F. Gerbier, A. Widera, S. Fölling, O. Mandel, and I. Bloch, Phys. Rev. A 73, 041602(R) (2006).
  * Santos  _et al._ [2007] L. Santos, M. Fattori, J. Stuhler, and T. Pfau, Phys. Rev. A 75, 053606 (2007).
  * Klausen  _et al._ [2001] N. N. Klausen, J. L. Bohn, and C. H. Greene, Phys. Rev. A 64, 053602 (2001).
  * Uchino  _et al._ [2010] S. Uchino, M. Kobayashi, and M. Ueda, Phys. Rev. A 81, 063632 (2010).
  * Saito  _et al._ [2013] H. Saito, Y. Kawaguchi, and M. Ueda, J. Phys. Condens. Matter 25, 404212 (2013).
  * Symes  _et al._ [2016] L. M. Symes, R. I. McLachlan, and P. B. Blakie, Phys. Rev. E 93, 053309 (2016).
  * Turban and Iglói [2002] L. Turban and F. Iglói, Phys. Rev. B 66, 014440 (2002).
  * Continentino and Ferreira [2004] M. A. Continentino and A. S. Ferreira, Phys. Stat. Mech. Its Appl. 339, 461 (2004).
  * Continentino [2017] M. Continentino,  _Quantum Scaling in Many-Body Systems: An Approach to Quantum Phase Transitions_, 2nd ed. (Cambridge University Press, Cambridge, 2017).
  * Pellegrini  _et al._ [2008] F. Pellegrini, S. Montangero, G. E. Santoro, and R. Fazio, Phys. Rev. B 77, 140404(R) (2008).
  * Divakaran  _et al._ [2008] U. Divakaran, A. Dutta, and D. Sen, Phys. Rev. B 78, 144301 (2008).
  * [94] G. Lagnese, F. M. Surace, S. Morampudi, and F. Wilczek, Detecting a long lived false vacuum with quantum quenches, arXiv:2308.08340 [cond-mat.stat-mech] .