# Caustics in the sine-Gordon model from quenches in coupled 1D Bose gases

Aman Agarwal<EMAIL_ADDRESS>Department of Physics and Astronomy,
McMaster University, 1280 Main St. W., Hamilton, Ontario, Canada L8S 4M1
BITS-Pilani, K. K. Birla Goa Campus, NH17B, Bypass Road, Zuarinagar, Goa
403726, India  International Centre for Theoretical Sciences, Tata Institute
of Fundamental Research, Bengaluru – 560089, India  Perimeter Institute for
Theoretical Physics, Waterloo, Ontario, Canada, N2L 2Y5  Department of
Physics, University of Guelph, Guelph, Ontario, Canada, N1G 2W1  Institute of
Physics, University of Greifswald, 17489 Greifswald, Germany  Manas Kulkarni
<EMAIL_ADDRESS>International Centre for Theoretical Sciences,
Tata Institute of Fundamental Research, Bengaluru – 560089, India  D. H. J.
O’Dell<EMAIL_ADDRESS>Department of Physics and Astronomy, McMaster
University, 1280 Main St. W., Hamilton, Ontario, Canada L8S 4M1

(August 27, 2024)

###### Abstract

Caustics are singularities that occur naturally in optical, hydrodynamic and
quantum waves, giving rise to high amplitude patterns that can be described
using catastrophe theory. In this paper we study caustics in a statistical
field theory setting in the form of the sine-Gordon model that describes a
variety of physical systems including coupled 1D superfluids. Specifically, we
use classical field simulations to study the dynamics of two ultracold 1D Bose
gases (quasi-condensates) that are suddenly coupled to each other and find
that the resulting non-equilibrium dynamics are dominated by caustics. Thermal
noise is included by sampling the initial states from a Boltzmann distribution
for phononic excitations. We find that caustics pile up over time in both the
number and phase difference observables leading to a characteristic non-
thermal ‘circus tent’ shaped probability distribution at long times.

coupled condensates, sine-Gordon hamiltonian, cigar shaped 1-D BEC, caustics

##  I Introduction

Wave focusing is ubiquitous in nature and leads to localized regions of high
amplitude called caustics that dominate wavefields. Everyday examples are
provided by rainbows and also the bright lines on the bottom of water pools
which are caused by the focusing of sunlight by raindrops and surface water
waves, respectively [1]. Caustics also occur in water waves themselves as ship
wakes [2] and more dramatically as tsunamis (focused by the topography of the
seabed [3, 4, 5]) and tidal bores (focused by v-shaped bays [6]).
Astrophysical examples include gravitational lensing by matter and the
twinkling of starlight due to time-dependent fluctuations in the density of
Earth’s atmosphere. Natural focusing also leads to the phenomenon of branched
flow [7] and is speculated to have given rise to the filamented nature of the
large scale structure of the universe [8, 9, 10, 11]. In all these systems
caustics give rise to extreme amplitude fluctuations that occur more
frequently than those predicted by gaussian statistics [12].

A remarkable property of caustics is that they commonly take on particular
characteristic shapes. This is because caustics are singularities of the ray
description, i.e. they are places where two or more rays coalesce leading to a
diverging intensity in the short wavelength limit [13]. Such singularities are
described by Thom’s catastrophe theory which rigorously shows that only
certain shapes of singularity are structurally stable against perturbations
and hence occur under ‘natural’ or generic conditions [14, 15, 16]. These
special shapes or catastrophes form a hierarchy organized by dimension where
the higher ones contain the lower ones. Each member of the hierarchy
represents a class of equivalent shapes that can be smoothly transformed into
each other, but each class is distinct and cannot be smoothly transformed into
any of the others. In two dimensions the only structurally stable shape is the
cusp and we shall see it appear frequently when we plot quantities such as
number fluctuations versus time. It is worth noting in this context that the
humble point focus that we associate with lensing is structurally unstable and
unfolds into an extended caustic in the presence of perturbations
(aberrations). Natural lenses are of course never perfect and so typically
produce the shapes predicted by catastrophe theory. The upshot of all this is
that caustics represent a form of universality in nonequilibrium wave
dynamics: they fall into equivalence classes each with their own shapes and
scaling properties analogous to, but a generalization of, equilibrium phase
transitions [13, 17].

Caustics should equally be present in quantum waves where, due to the
probabilistic interpretation, they correspond to regions of high probability
density. Quantum matter wave caustics have been seen in experiments with cold
neutrons [18, 19], electron microscopes [20], atom optics [21, 22, 23], and
most recently in atom lasers [24]. Theoretical works on such matter wave
caustics have also considered their ‘fine structure’ [13] which features a
lattice of vortices [25, 26, 27]. Quantum fields are another area where
caustics are expected to form naturally during dynamics. Early work centred on
the electromagnetic field [28, 29], including an interpretation of Hawking
radiation as a ‘quantum catastrophe’ [30], and more recently this idea has
been extended to quantum many-particle systems including bosonic Josephson
junctions [26, 31, 32], the XY model with long-range interactions (Hamiltonian
mean field model) [33], quantum spin chains [27] and the Bose-Hubbard model
[34]. One point to appreciate is that the caustics in many-body systems can
occur in the wavefunction associated with an entire $N$-body configuration.
Quantum many-particle caustics therefore live in Fock space which can have a
large number of dimensions and hence lead to very complicated catastrophes
[34]. However, catastrophes obey projection identities which means that when
projected down to lower dimensions one obtains either the same catastrophe or
one lower down the hierarchy [35]. Thus, low order correlation functions
obtained by integrating out most of the degrees of freedom will also
generically contain caustics [27].

Figure 1: Schematic of the setup we consider. The top figure shows two quasi
one-dimensional gases that are prepared independently and then suddenly
coupled together. We call this process of sudden coupling a “$J$-quench”.
$\rho_{1}(z)$ and $\rho_{2}(z)$ represent the density (red) in the first and
second condensates, respectively. Similarly, $\phi_{1}(z)$ and $\phi_{2}(z)$
represent the phases (black) of the two condensates. Prior to the $J$-quench,
these fields in the two condensates are independent and contain thermal
fluctuations. The bottom figure shows how a $J$-quench could be implemented by
suddenly reducing the tunneling barrier height in a double well potential from
a higher to a lower value.

In this paper we study caustics in the sine-Gordon (SG) model. The (classical)
SG model obeys the nonlinear wave equation

$\frac{\partial^{2}\phi}{\partial
t^{2}}-c_{0}^{2}\frac{\partial^{2}\phi}{\partial
z^{2}}+\omega_{0}^{2}\sin\phi=0$ (1)

where $\phi=\phi(z,t)$ is a one dimensional field, and $c_{0}$ and
$\omega_{0}$ represent a characteristic speed and frequency, respectively. If
$c_{0}$ is taken to be the speed of light then Eq. (1) is relativistically
covariant, being a nonlinear version of the Klein-Gordon equation and reducing
to it when $\phi\ll 1$ such that $\sin\phi\approx\phi$. The SG model received
attention from the high energy physics community in the 1970s due its soliton
solutions [36, 37, 38, 39], but also describes the low energy physics of a
considerable range of condensed matter systems including crystal dislocations
[40], domain walls in magnetic [41] and binary superfluid [42] systems, the
Heisenberg spin chain with a field induced gap [43, 44, 45], one-dimensional
Bose gases in periodic potentials (that can capture the Mott-insulator to
superfluid transition in one dimension) [46, 47], two-dimensional Bose gases
realizing the XY model [48], trapped ions [49], and two tunnel-coupled one-
dimensional Bose gases [50, 51, 52, 53, 54, 55, 56, 57]. The fact that the SG
model is both nonlinear and integrable means that attention is often focused
on its soliton solutions, but part of our mission in this paper is to point
out that these same properties also imply that caustics (which are associated
with the existence of tori in phase space [58]) are expected to occur
generically, and we are aware of only one previous study of caustics in this
model [59].

The particular physical realization we have in mind for this paper is a system
composed of two elongated quasi-one dimensional Bose gases coupled by
tunneling along their length; the field $\phi(z,t)$ in Eq. (1) gives the
relative phase between the two quantum gases. Quasi one-dimensional Bose gases
have been created in a number of experiments over the last two decades using
tightly trapped ultracold atoms, and the remarkable tunability of these
systems allows the strongly interacting Tonks-Girardeau regime [60, 61], the
weakly interacting quasi-condensate regime [62, 63, 64, 65], and also the
crossover between the two [66, 67], to be reached. It is important to note
that, in accordance with the Mermin-Wagner theorem [68], one-dimensional Bose
gases do not undergo true Bose-Einstein condensation at low temperature,
unlike three dimensional gases. Instead, they can form quasi-condensates where
density fluctuations are suppressed but phase fluctuations remain [69, 70]. In
this paper we shall work in the weakly interacting regime and assume a state
of the system consisting of a quasi-condensate plus small thermal
fluctuations.

A system comprised of two coupled quasi-one dimensional gases can be made by
taking a single gas and splitting it in two along its long axis by switching
on an elongated double well potential. This is the experimental protocol
typically adopted in a series of experiments conducted by the Vienna group
[63, 71, 72, 73, 74, 75, 76, 77]. The combination of almost complete isolation
from the environment, long relaxation times and spatially resolved
measurements of phase and number difference make these experiments ideal for
investigating many-particle quantum dynamics, including fundamental questions
such as whether and how closed quantum systems reach equilibrium. The gas can
be split slowly so that it always remains close to equilibrium leading to
number squeezed states [78, 79] or it can be split rapidly, leading to a so-
called quantum quench which launches the system into a nonequilibrium state.

In this paper we shall consider the opposite quench where two one-dimensional
gases are suddenly connected together (see schematic representation in Figure
1). This touches on rather fundamental considerations in quantum mechanics
since it describes the build-up of coherence between two initially independent
systems, and is therefore related to the double-slit experiment for many-
particle systems [80, 81, 82, 83, 53]. We shall refer to this as a
“$J$-quench” because $J$ is often used to denote the coupling strength between
the two wells. In a simple two-mode description of a bosonic Josephson
junction, i.e. one that assumes a single mode in each well without the quasi-
continuum of low energy longitudinal modes that are present in highly
elongated traps, such a quench is predicted to result in a periodic collapse
and revival of the atom number distribution between the two wells [84, 85,
86]. Essentially the same behavior, but $\pi/2$ out of phase, occurs in the
relative phase which is the conjugate variable to number difference. In Refs.
26, 31, 32 these revivals are shown to be examples of quantum caustics in a
many-particle system. One of our main aims here is to investigate what happens
to these caustics in the presence of the dispersive longitudinal modes present
in the SG model, and is part of a wider program attempting to understand the
role of caustics in quantum many particle dynamics [31, 17, 26, 33, 32, 27,
34].

Due to the difficulty of solving the fully quantum SG model we take a
semiclassical-style approach based on classical field configurations which are
solutions of Eq. (1). Each configuration is analogous to a single geometric
ray in optics and we include fluctuations by summing many configurations. The
initial conditions for each field configuration are randomly sampled from a
Boltzmann distribution. This approach is similar in spirit to the truncated
Wigner approximation (TWA) [87, 88, 89, 90, 91, 92] which includes quantum
fluctuations around the classical field by summing many rays sampled from a
quantum probability distribution. The TWA has previously been applied to one-
dimensional Bose gases by Martin and Ruostekoski [93, 94] who studied dark
solitons, and also to the connection problem of two zero temperature one-
dimensional Bose gases by Dalla Torre, Demler and Polkovnikov [53], who
proposed a universal scaling form for the phase dynamics after the quench.
More recently, the TWA has been used by Horváth et al. [95] to study the
surprisingly sudden relaxation of the phase seen in the Vienna BEC splitting
experiments [77]. In this paper, we include both the quantum fluctuations
arising from coupling two independent systems and thermal fluctuations arising
from thermal phonons in the longitudinal modes and compare the time evolution
of macroscopic variables (the total number difference and phase difference) in
the SG system against the simpler two mode system [31, 17, 26]. We find that
following a quench caustics dominate the dynamics of the macroscopic variables
of both systems, even in the presence of thermal fluctuations. Due to the
singular nature of caustics, and combined with their structural stability, we
therefore propose that strong nongaussian fluctuations are a generic
phenomenon following a quench in the SG model (and indeed, in integrable or
moderately chaotic many-body systems in general).

The caustics we discuss in this paper also have implications for the question
of relaxation towards equilibrium at long times in many particle systems.
While chaotic (nonintegrable) and open quantum systems should thermalize
(although a complete description is still the subject of active research [96,
97, 98, 99, 100, 101, 102, 103]), closed integrable models do not reach a
conventional Gibbs state. We show here that in the SG model there is a pile-up
of caustics leading to a singular shape for the long time probability
distribution for the macroscopic variables that resembles the shape of a
circus tent and is quite distinct from the thermal equilibrium prediction. We
find that an analytic approximation to the singular distribution based on an
ergodic pendulum (assuming a microcanonical or ‘equal-probability’
distribution) provides a good fit to the numerical data.

The plan for the rest of this paper is as follows. We start in Sec. II by
deriving the SG hamiltonian from the many-body description of two coupled 1D
Bose gases. In Sec. III we describe the natural length and time scales and use
them to write the SG hamiltonian and equations of motion in convenient
dimensionless forms. Subsequently, in Sec. IV we develop a method for finding
the initial conditions for the SG equations of motion. We assume that prior to
the quench the two Bose gases are independent and at thermal equilibrium with
a bath at temperature $T$. The initial conditions are obtained by
stochastically sampling the Fourier modes of a 1D quasicondensate obeying the
Tomonaga-Luttinger liquid theory. With the initial conditions in hand, in Sec.
V we give the main results of this paper which are the dynamics of the
macroscopic number and phase difference variables obtained by solving the
equations of motion numerically. In Sec. VI we consider the bigger picture and
examine the universal aspects of our results including the influence of
caustics on the coherence as well as the long time dynamics and the
establishment of (non-thermal / non-Gaussian) equilibrium. We conclude in Sec.
VII. There are also six appendices where we give the details of the
calculations as well as bench marking our numerical method.

##  II From two coupled condensates to the sine-Gordon plus model

We begin by deriving the SG model as an effective low energy description for
two coupled one-dimensional Bose gases. For the sake of clarity, we list the
main simplifications employed in this work:

  * •

the treatment of a quantum many body problem by a semiclassical method (TWA).

  * •

the neglect of a weak harmonic trap along the long axis which would otherwise
lead to a non-uniform longitudinal density (this can be avoided in box traps
which, although rarer, can be realized [76, 104])

  * •

the assumption of a constant value for the tunnel coupling $J$ along the
entire length of the gases

  * •

the neglect of coupling to symmetric and higher transverse modes. Some more
involved theoretical models do include these effects [56, 57].

These simplifications are not expected to qualitatively alter the main results
of this work due to the structural stability of caustics. In other words,
caustics are known to be robust to perturbations in both the Hamiltonian and
initial conditions.

A theoretical description of two ultracold quasi-one dimensional gases made up
of bosonic atoms of mass $m$, and held parallel to each other so that the
atoms can tunnel between them at rate $J$, can be obtained from the following
microscopic Hamiltonian [50, 74, 51]

$\begin{split}\hat{H}=\sum_{j=1,2}\int_{{{-L/2}}}^{{{L/2}}}dz&\left[-\frac{\hbar^{2}}{2m}\hat{\psi}_{j}^{\dagger}(z)\frac{\partial^{2}}{\partial
z^{2}}\hat{\psi}_{j}(z)+U(z)\,\hat{\psi}_{j}^{\dagger}(z)\hat{\psi}_{j}(z)+\frac{g_{1D}}{2}\,\hat{\psi}_{j}^{\dagger}(z)\hat{\psi}_{j}^{\dagger}(z)\hat{\psi}_{j}(z)\hat{\psi}_{j}(z)\right]\\\
-&\int_{-L/2}^{L/2}dz\,\hbar
J\,\left[\hat{\psi}_{1}^{\dagger}(z)\hat{\psi}_{2}(z)+\hat{\psi}_{2}^{\dagger}(z)\hat{\psi}_{1}(z)\right]\
.\end{split}$ (2)

The indices $j=1,2$ label the two gases and each is assumed to be tightly
trapped in the $x$ and $y$ directions so that those degrees of freedom are
frozen into their ground states. Only the longitudinal degree of freedom $z$
in each gas is taken to be active. In experiments there will usually be a weak
longitudinal trapping potential $U(z)$, although as mentioned above for
simplicity we set it to zero and hence consider a uniform system of length $L$
with periodic boundary conditions. The quantum field operator
$\hat{\psi}_{j}(z)$ annihilates a particle at point $z$ and together with its
hermitian conjugate obeys bosonic commutation relations
$[\hat{\psi}_{j}(z),\hat{\psi}_{j^{\prime}}^{\dagger}(z^{\prime})]=\delta_{jj^{\prime}}\delta(z-z^{\prime})$.
The interaction constant $g_{1D}$ characterizes the effect of atom-atom
scattering within each gas on the longitudinal degree of freedom and can be
controlled both in magnitude and sign either through Feshbach or confinement-
induced scattering resonances [105]. We note in passing that a possible
alternative physical realization of this problem could be a spinor Bose gas in
a single quasi-one dimensional trap [106]. In fact, bosonic Josephson
junctions where the atoms are held in a single trap and two atomic spin states
are used for the two states have already been realized experimentally [107].

A weakly interacting three-dimensional Bose gas at ultracold temperatures will
undergo Bose-Einstein condensation and can be described to high accuracy by a
classical field approximation (Gross-Pitaevskii theory [108]). In a quasi-one
dimensional geometry quantum fluctuations can still be small if the density is
not too low, and under these circumstances the gas can be treated as a quasi-
condensate where the quantum field operators are replaced by classical fields
[69, 109, 110]

$\hat{\psi}_{j}(z)\rightarrow\psi_{j}(z)=\sqrt{n_{1D}+\rho_{j}(z)}\,\exp[i\phi_{j}(z)]\
.$ (3)

Here $n_{1D}=N/L$ is the background density where $N$ is the number of atoms
in each gas (for simplicity we assume an equal number of atoms $N$ in each
gas; the structural stability of caustics means that they are stable to small
differences in $n_{1D}$ between the two gases). $\rho_{j}(z)$ and
$\phi_{j}(z)$ are the atom number density and phase fluctuations at each point
$z$, respectively. These are canonically conjugate variables and can even be
quantized in a semiclassical regime such that they obey the commutation
relations
$[\hat{\rho}_{j}(z),\hat{\phi}_{j^{{}^{\prime}}}(z^{\prime})]\approx\delta_{jj^{\prime}}\delta(z-z^{\prime})$
in a coarse grained sense [110]. However, in the present paper $\rho_{j}(z)$
and $\phi_{j}(z)$ will be purely classical fields subject only to thermal
fluctuations.

We can further decompose the fields into their symmetric and antisymmetric
components

$\displaystyle\rho_{s}(z)=\frac{\rho_{1}(z)+\rho_{2}(z)}{2},\quad\rho_{a}(z)=\frac{\rho_{1}(z)-\rho_{2}(z)}{2}$
$\displaystyle\phi_{s}(z)=\phi_{1}(z)+\phi_{2}(z),\quad\phi_{a}(z)=\phi_{1}(z)-\phi_{2}(z)\
.$ (4)

If the fluctuations are small $\rho_{a}$ will be small whereas $\rho_{s}$ will
be comparatively large. The particle-particle interaction energy will then
typically cause the dynamics of the symmetric modes to occur at higher energy
than the antisymmetric ones, and consequently we can ignore the symmetric
degrees of freedom as long as we restrict attention to low energies [50, 55,
72, 75]. The Hamiltonian purely describing the antisymmetric variables is (see
Appendix A for details)

$\begin{split}H_{\mathrm{SG}+}=&\int_{-L/2}^{L/2}dz\bigg{[}g_{1D}\,\rho_{a}^{2}(z)+\frac{\hbar^{2}n_{1D}}{4m}\left(\frac{\partial\phi_{a}}{\partial
z}\right)^{2}\\\
&+\frac{\hbar^{2}}{4mn_{1D}}\left(\frac{\partial\rho_{a}}{\partial
z}\right)^{2}-2\hbar Jn_{1D}\cos\phi_{a}(z)\bigg{]}.\end{split}$ (5)

We refer to this as the “sine-Gordon plus” (SG+) Hamiltonian because it
includes an extra term (the third term) in comparison to the standard SG
Hamiltonian. This term involves gradients of density fluctuations and results
in an energy cost which automatically suppresses density fluctuations at small
length scales. It is also worth noting that including this term means that the
density and phase fluctuations [the second term in Eq. (5)] are incorporated
on an equal footing. This is also in accordance with Gross-Pitaevskii theory
which suppresses density fluctuations with wavelengths below the healing
length [95]

$\xi_{h}=\frac{\hbar}{\sqrt{mg_{1D}n_{1D}}}\ .$ (6)

However, when $n_{1D}$ is relatively large the third term is naturally
suppressed in comparison to the others and can be dropped as long as the
density gradients are small leading to the SG Hamiltonian [55, 74]

$\begin{split}H_{\mathrm{SG}}=\int_{-L/2}^{L/2}dz\,\bigg{[}&g_{1D}\,\rho_{a}(z)^{2}+\frac{\hbar^{2}n_{1D}}{4m}\left(\frac{\partial\phi_{a}}{\partial
z}\right)^{2}\\\ &-2\hbar Jn_{1D}\cos\phi_{a}(z)\bigg{]}\ .\end{split}$ (7)

The nonlinear piece in both Hamiltonians is the cosine term which originates
from tunneling between the two wells and occurs in all Josephson junction type
problems. It provides an effective potential well for phase configurations
$\phi(z,t)$ that play the role of rays. In fact, as we shall see in Section V,
it acts as an (imperfect) lens that focuses rays excited by the quench to form
caustics. For the sake of brevity, and when we deem no confusion can arise, we
will omit the ‘$a$’ subscript on antisymmetric variables since we will not be
dealing with symmetric degrees of freedom.

The fact that the two fields $\phi(z)$ and $\rho(z)$ form a conjugate pair
means that their equations of motion are given by Hamilton’s equations

$\begin{split}\dot{\phi}&=\frac{1}{\hbar}\frac{\delta\mathcal{H}}{\delta\rho(z)}\\\
\dot{\rho}&=-\frac{1}{\hbar}\frac{\delta\mathcal{H}}{\delta\phi(z)}\end{split}$
(8)

where $\mathcal{H}$ is the Hamiltonian density defined via

$H=\int_{-L/2}^{L/2}{\mathcal{H}\,dz}.$ (9)

Applying these equations to the SG+ Hamiltonian given in Eq. (5) we find the
following of equations of motion

$\begin{split}\frac{d\phi(z,t)}{dt}&=2\ \frac{g_{1D}}{\hbar}\rho(z,t)+2\
\frac{\hbar}{4mn_{1D}}\frac{\partial^{2}\rho(z,t)}{\partial z^{2}}\\\
\frac{d\rho(z,t)}{dt}&=2\ \frac{\hbar
n_{1D}}{4m}\frac{\partial^{2}\phi(z,t)}{\partial
z^{2}}-2Jn_{1D}\sin[\phi(z,t)]\ .\end{split}$ (10)

These are the key equations we use to solve for the dynamics of the field
configurations. They have the form of Josephson’s equations [111] augmented by
second order spatial derivatives $\partial^{2}\phi/\partial z^{2}$ and
$\partial^{2}\rho/\partial z^{2}$ which account for phase and density
fluctuations along the longitudinal direction. Combined with the sine term,
they will cause wavepackets to disperse along $z$. In the absence of these
terms we have exactly the equations of motion for a pendulum where $\phi$ is
the angular displacement from equilibrium and $\rho$ plays the role of angular
momentum. The dependence on $z$ suggests an interpretation in terms of a
continuous chain of many pendula each coupled to its neighbors by the spatial
derivative terms and is reminiscent of the Fermi-Pasta-Ulam-Tsingou problem
[50, 112].

In this paper the coupled equations of motion given in Eq. (10) will be solved
numerically for a system of length $L$. To perform the numerical computations
we discretize the system on a spatial grid with $N_{L}+1$ points which makes
the grid spacing $a=L/N_{L}$. The positions of the grid points are given by
$z=ra$ where $r$ is an integer

$r=-\frac{N_{L}}{2},\ldots,\frac{N_{L}}{2}$ (11)

and $N_{L}$ is chosen to be an even integer.

There is in fact a physical limitation on the grid size. Eq. (10) is classical
and valid only on length scales greater that healing length $\xi_{h}$ [95,
51]. Therefore, any numerics performed on Eq. (10) are meaningful only when
the lattice grid size $a$ is greater than $\xi_{h}$. In particular, $N_{L}$
should be such that $a>\xi_{h}$ which implies

$N_{L}^{2}<\frac{mg_{1D}n_{1D}L^{2}}{\hbar^{2}}.$ (12)

We fulfil the condition given in Eq. (12) in our numerics.

##  III Natural scales

Let us express the SG/SG+ Hamiltonians and equations of motion in terms of the
natural scales for a one-dimensional quantum fluid. For a length scale we
chose the healing length $\xi_{h}$ given in Eq. (6). The ratio of the healing
length to the mean interparticle spacing $1/n_{1D}$ motivates the definition
of the Luttinger parameter

$K=\sqrt{\frac{n_{1D}(\hbar\pi)^{2}}{4g_{1D}m}}.$ (13)

This dimensionless quantity measures how strongly interacting the system is -
when $K\gg 1$ the healing length is much greater than the interparticle
spacing and the system is in the weakly interacting (quasi-condensate) regime.
Another key physical quantity is the speed of sound

$c=\sqrt{\frac{g_{1D}n_{1D}}{m}}.$ (14)

This can be used to define a characteristic energy, namely that associated
with phonons (quanta of sound)

$E=\hbar\omega=\frac{\hbar c}{\xi_{h}}$ (15)

where we have set the natural frequency $\omega$ to be the ratio of the speed
of sound to the healing length.

We therefore transform to the following dimensionless variables

$\begin{split}z&\xrightarrow{}\tilde{z}=\frac{z}{\xi_{h}}\quad,\quad
t\xrightarrow{}\tilde{t}=\frac{c}{\xi_{h}}t\\\
\rho&\xrightarrow{}\tilde{\rho}=\rho\;\xi_{h}\quad,\quad\phi\xrightarrow{}\tilde{\phi}=\phi\end{split}$
(16)

and defining $\tilde{H}_{\mathrm{SG}}=H_{\mathrm{SG}}/E$ and likewise for
$\tilde{H}_{\mathrm{SG}+}$ we obtain the two Hamiltonians in dimensionless
form

$\tilde{H}_{\mathrm{SG}}=\int_{-L/2}^{L/2}d\tilde{z}\,\bigg{[}\Gamma\,\tilde{\rho}^{2}+\epsilon\left(\frac{\partial\tilde{\phi}}{\partial\tilde{z}}\right)^{2}-2\mathcal{J}\cos\tilde{\phi}\bigg{]}$
(17)

and

$\displaystyle\tilde{H}_{\mathrm{SG}+}$ $\displaystyle=$
$\displaystyle\int_{-L/2}^{L/2}d\tilde{z}\,\bigg{[}\Gamma\,\tilde{\rho}^{2}+\epsilon\left(\frac{\partial\tilde{\phi}}{\partial\tilde{z}}\right)^{2}+\frac{\Gamma}{4}\left(\frac{\partial\tilde{\rho}}{\partial\tilde{z}}\right)^{2}$
(18) $\displaystyle-$ $\displaystyle 2\mathcal{J}\cos\tilde{\phi}\bigg{]}$

where the coefficients are given by

$\Gamma=\frac{\pi}{2K}\ ,\ \epsilon=\frac{K}{2\pi}\ ,\
\mathcal{J}=\frac{K}{2\pi}\frac{\xi_{h}^{2}}{\xi_{s}^{2}}\ .$ (19)

In the last term we have introduced the spin healing length

$\xi_{s}=\sqrt{\frac{\hbar}{4mJ}}$ (20)

which provides a measure for the distance over which coherence between the two
gases is restored due to the tunnel coupling $J$ [55]. At finite temperatures
another useful length scale is the thermal phase coherence length

$\lambda_{T}=\frac{2\hbar^{2}n_{1D}}{mk_{B}T}.$ (21)

The dimensionless form of the equations of motion can now be given. For the SG
model we find

$\begin{split}\frac{d\tilde{\phi}}{d\tilde{t}}&=2\Gamma\tilde{\rho}\\\
\frac{d\tilde{\rho}}{d\tilde{t}}&=2\epsilon\frac{\partial^{2}\tilde{\phi}}{\partial\tilde{z}^{2}}-2\mathcal{J}\sin\tilde{\phi}\end{split}$
(22)

and for the SG+ model we obtain

$\begin{split}\frac{d\tilde{\phi}}{d\tilde{t}}&=2\Gamma\tilde{\rho}-\frac{\Gamma}{2}\frac{\partial^{2}\tilde{\rho}}{\partial\tilde{z}^{2}}\\\
\frac{d\tilde{\rho}}{d\tilde{t}}&=2\epsilon\frac{\partial^{2}\tilde{\phi}}{\partial\tilde{z}^{2}}-2\mathcal{J}\sin\tilde{\phi}\
.\end{split}$ (23)

##  IV Initial Conditions

The dynamics we seek to study in this paper start from a $J$-quench where two
independent one-dimensional gases at thermal equilibrium are suddenly coupled.
In order to obtain the initial density and phase fluctuations of these gases
we use the Tomonaga-Luttinger (TL) model that provides the universal low
energy effective theory for one-dimensional systems (low energy limit of the
Lieb-Lininger theory, for example) [53].

###  IV.1 Tomonaga-Luttinger (TL) liquid

In our notation the TL Hamiltonian reads

$\begin{split}H_{\mathrm{TL}}&=\int_{-L/2}^{L/2}{dz\left[g_{1D}\rho_{j}(z)^{2}+\frac{\hbar^{2}n_{1D}}{4m}\left(\frac{\partial\phi_{j}}{\partial
z}\right)^{2}\right]}\end{split}$ (24)

where $j$ labels either of the two gases. We henceforth, omit this label for
the sake of brevity with the understanding that in this section the density
and phase fields refer to just one of the two gases. Eq. (24) has the same
mathematical structure as the SG model but without the tunnelling term. If we
include density fluctuations we find

$\begin{split}H_{\mathrm{TL}+}=\int_{-L/2}^{L/2}dz\bigg{[}g_{1D}\rho(z)^{2}&+\frac{\hbar^{2}n_{1D}}{4m}\left(\frac{\partial\phi}{\partial
z}\right)^{2}\\\
&+\frac{\hbar^{2}}{4mn_{1D}}\left(\frac{\partial\rho}{\partial
z}\right)^{2}\bigg{]}\ .\end{split}$ (25)

The TL model is quadratic and hence its thermal fluctuations can be treated
exactly. To this end it is useful to work in Fourier space and we apply
discrete Fourier transforms defined on the numerical grid with $N_{L}$ points
as discussed at the end of Section II. The phase field $\phi$ and its Fourier
transform $\varphi$ are related by

$\begin{split}&\phi_{r}=\frac{1}{\sqrt{N_{L}+1}}\sum_{k=-N_{L}/2}^{N_{L}/2}\varphi_{k}\
\exp\left[i\frac{2\pi kr}{N_{L}+1}\right]\\\
&\varphi_{k}=\frac{1}{\sqrt{N_{L}+1}}\sum_{r=-N_{L}/2}^{N_{L}/2}\phi_{r}\
\exp\left[-i\frac{2\pi kr}{N_{L}+1}\right]\ .\end{split}$ (26)

The discrete data
$\\{\phi_{r}\\}=\\{\phi_{-N_{L}/2},\ldots,\phi_{0},\ldots,\phi_{N_{L}/2}\\}$
and its transform are located symmetrically about $r=0$ and $k=0$,
respectively. Since the value $\phi_{r}$ of the field at each coordinate space
grid point is a real number the condition $\varphi_{-k}=\varphi_{k}^{\ast}$
must hold. Similarly the density fluctuation field $\rho$ and its Fourier
transform $\varrho$ are related by

$\begin{split}&\rho_{r}=\frac{1}{\sqrt{N_{L}+1}}\sum_{k=-N_{L}/2}^{N_{L}/2}\varrho_{k}\
\exp\left[i\frac{2\pi kr}{N_{L}+1}\right]\\\
&\varrho_{k}=\frac{1}{\sqrt{N_{L}+1}}\sum_{r=-N_{L}/2}^{N_{L}/2}\rho_{r}\
\exp\left[-i\frac{2\pi kr}{N_{L}+1}\right]\end{split}$ (27)

where again the reality of the field in coordinate space requires that
$\varrho_{-k}=\varrho_{k}^{\ast}$. Inserting these transformations in Eq. (25)
we obtain (see Appendix B for details)

$\begin{split}H_{\mathrm{TL}+}&=a\,g_{1D}\,\sum_{k=-N_{L}/2}^{N_{L}/2}\,|\varrho_{k}|^{2}\\\
&+a\ \hbar\ n_{1D}\sum_{k=-N_{L}/2}^{N_{L}/2}\frac{\hbar\pi^{2}k^{2}}{mL^{2}}\
|\varphi_{k}|^{2}\\\ &+a\
\frac{\hbar^{2}}{4mn_{1D}}\sum_{k=-N_{L}/2}^{N_{L}/2}\frac{4\pi^{2}k^{2}}{L^{2}}\
|\varrho_{k}|^{2}\ .\end{split}$ (28)

Before proceeding with further analysis of Eq. (28), it is worth noting that
it can be recast in a standard Luttinger liquid form

$\begin{split}H_{\mathrm{TL}+}=\frac{ac\hbar}{2}\sum_{k=-N_{L}/2}^{N_{L}/2}\,\bigg{[}\frac{K}{\pi}&\frac{4\pi^{2}k^{2}}{L^{2}}|\varphi_{k}|^{2}+\frac{\pi}{K}|\varrho_{k}|^{2}\\\
&+\frac{K}{\pi}\frac{4\pi^{2}k^{2}}{N^{2}}|\varrho_{k}|^{2}\bigg{]}\end{split}$
(29)

where the strength of the terms depends either on $K$ or $1/K$.

Applying the transformations given in Eq. (16), the Fourier space variables
can be written in dimensionless form as

$\begin{split}\varrho_{k}\xrightarrow{}\tilde{\varrho}_{k}=\xi_{h}\varrho_{k}\quad,\quad\varphi_{k}\xrightarrow{}\tilde{\varphi}_{k}=\varphi_{k}\
\end{split}$ (30)

and the TL+ Hamiltonian given in Eq. (28) scaled by the energy $E=\hbar
c/\xi_{h}$ is given by

$\begin{split}\tilde{H}_{\mathrm{TL}+}=\frac{\tilde{L}}{N_{L}}\sum_{k=-N_{L}/2}^{N_{L}/2}\,\bigg{[}&\frac{\epsilon\,4\pi^{2}k^{2}}{\tilde{L}^{2}}|\tilde{\varphi}_{k}|^{2}+\Gamma|\tilde{\varrho}_{k}|^{2}\\\
&+\frac{\Gamma\,\pi^{2}k^{2}}{\tilde{L}^{2}}|\tilde{\varrho}_{k}|^{2}\bigg{]}\end{split}$
(31)

where $\tilde{L}=L/\xi_{h}$ is the ratio of the system size to the healing
length. Comparison with the spatial version of $H_{\mathrm{TL+}}$ given in Eq.
(25) shows where this factor comes from: as the size is increased the range of
the integration increases linearly and this is accounted for by $\tilde{L}$ in
the Fourier transformed version. Note that all parameters and variables in Eq.
(31) are dimensionless.

###  IV.2 Thermal equilibrium

To find the initial conditions on the fields $\rho_{j}(z)$ and $\phi_{j}(z)$
we assume that each gas is at thermal equilibrium such that the excitation
(phonon) modes of the TL+ Hamiltonian are populated with a probability given
by the Boltzmann distribution. The range of temperatures we simulate is listed
in Table 1 along with the values of all the other key parameters, and is
chosen so as to correspond to realistic experimental conditions (the
temperature must be low enough that the quasi-condensate description is
valid).

In the canonical ensemble of statistical mechanics the probability that a
system at thermal equilibrium has the phase space configuration
$s={q_{1},p_{1},q_{2},p_{2}...q_{N},p_{N}}$ is proportional to the Boltzmann
weight $\exp[-\beta H(s)]$, where $\beta=1/k_{B}T$ and
$H=\sum_{i}{p_{i}^{2}/2m+V(q_{i})}$. The Hamiltonian in Eq. (31) is quadratic
and hence the Boltzmann weight becomes that of a series of independent
harmonic oscillators

$\begin{split}e^{-\tilde{\beta}\tilde{H}_{\mathrm{TL+}}}=\prod_{k}e^{-P_{k}^{2}/2\sigma^{2}_{\rho+}}\
e^{-Q_{k}^{2}/2\sigma_{\phi+}^{2}(k)}\end{split}$ (32)

where $\tilde{\beta}=(\hbar c/\xi_{h})/k_{B}T$ is the appropriately scaled
temperature parameter and we have introduced the real variables $Q_{k}$ and
$P_{k}$ which are related to the old variables by

$\tilde{\varphi}_{k}=Q_{k}e^{i\alpha_{k}},\,\quad\,\tilde{\varrho}_{k}=P_{k}e^{i\beta_{k}}.$
(33)

The phases $\alpha_{k}$ and $\beta_{k}$ allow for the fact that
$\tilde{\varphi}_{k}$ and $\tilde{\varrho}_{k}$ can be complex numbers. The
variances in Eq. (32) are given by

$\displaystyle\sigma^{2}_{\rho+}(k)$ $\displaystyle=$
$\displaystyle\frac{N_{L}}{2\tilde{\beta}}\frac{1}{\Gamma\tilde{L}(1+\pi^{2}k^{2}/\tilde{L}^{2})}$
(34) $\displaystyle\sigma^{2}_{\phi+}(k)$ $\displaystyle=$
$\displaystyle\frac{N_{L}}{2\tilde{\beta}}\frac{\tilde{L}}{4\pi^{2}k^{2}\epsilon}\
.$ (35)

The partition function can now be written down as

$\begin{split}Z=&\prod_{k}\int_{-\infty}^{\infty}e^{-\tilde{\beta}\tilde{H}_{\mathrm{TL}+}}\
dP_{k}dQ_{k}\\\
&=\prod_{k}\left(\sigma_{\rho+}\sqrt{2\pi}\right)\left(\sigma_{\phi+}(k)\sqrt{2\pi}\right)\end{split}$
(36)

and hence the probability $\mathcal{P}$ of a particular configuration
$(Q_{1},Q_{2},....,P_{1},P_{2},....)$ is

$\mathcal{P}=\prod_{k}\left(\frac{e^{-P_{k}^{2}/2\sigma^{2}_{\rho+}}}{\sigma_{\rho+}\sqrt{2\pi}}\right)\left(\frac{e^{-Q_{k}^{2}/2\sigma_{\phi+}^{2}(k)}}{\sigma_{\phi+}(k)\sqrt{2\pi}}\right)\
.$ (37)

This is seen to be the total probability distribution for independent random
variables $P_{k}$ and $Q_{k}$ drawn from normal distributions. Thus, the
absolute values of the Fourier coefficients $\tilde{\varrho}_{k}$ and
$\tilde{\varphi}_{k}$ are normally distributed random variables with zero mean
and variances given by Eqns. (34) and (35). We sample these numerically from
normal distributions to generate the initial system configuration. The phases
$\alpha_{k}$ and $\beta_{k}$ given in Eq. (33) do not appear in the Boltzmann
weight and are chosen randomly from the range $[0,2\pi)$. In fact, for both
the phases and the amplitudes we only need to choose the values for terms with
$k\geq 0$ because the reality conditions imply that we can put

$\displaystyle Q_{k}=Q_{-k}\ ,\quad P_{k}=P_{-k},$
$\displaystyle\alpha_{k}=-\alpha_{-k}\ ,\quad\beta_{k}=-\beta_{-k}\ .$ (38)

So far we have only considered the initial state of a single gas. By
subtracting the results for two gases we can obtain the initial values of the
antisymmetric variables $\rho_{a}(z)$ and $\phi_{a}(z)$ defined in Eq. (4).
Actually, due to the fact that the SG+ Hamiltonian with $J=0$ and expressed in
terms of antisymmetric variables as given in Eq. (5) formally has the same
structure as the TL+ Hamiltonian given in Eq. (25), sampling initial data for
two gases is unnecessary and one can obtain $\rho_{a}(z)$ and $\phi_{a}(z)$
directly by sampling them as though they were from one gas described by the
TL+ Hamiltonian. However, in doing so, consideration needs to be given to the
average value of relative phase $\phi_{a}(z)$ because both the SG+ and TL+
Hamiltonians only contain the spatial derivative of the phase but not the
phase itself. Its average value is therefore not determined by energy
considerations and is left to float freely. This is also apparent in the
Fourier transformed version of the TL Hamiltonian given in Eq. (31) where the
$k=0$ term involving $\tilde{\varphi}_{0}$ is absent due to the vanishing of
its coefficient which is proportional to $k^{2}$. To take into account the
random phase difference between the two gases one can chose
$\tilde{\varphi}_{0}$ to be a random number in the range $[-\pi\ldots\pi)$ but
multiplied by a factor of $\sqrt{N_{L}+1}$ in order to respect the
normalization in Eq. (26). This gives values of the average value of
$\phi_{a}(z)$ in the desired range $-\pi$ and $+\pi$.

The random value of the initial phase difference is actually a key feature of
the $J$-quench. It populates the cosine potential landscape in the Hamiltonian
with uniform probability. As the trajectories roll back and forth in this
potential they form caustics. In effect, the cosine potential acts as an
imperfect lens that focuses an initially flat ‘wavefront’ over time.

Symbol | Parameter | Value
---|---|---
$\omega_{\perp}$ | trapping frequency |  $2\pi\times 3$ kHz
$m$ | mass of Rb atom |  $1.41\times 10^{-25}$ kg
$a_{s}$ | scattering length |  $98\times 0.52$ Å
$N$ | number of atoms | 1200
$L$ | system length | 18 $\mu$m
$n_{1D}$ | average density | $6.7\times 10^{7}\mathrm{m}^{-1}$
$g_{1D}$ | 2 $\hbar a_{\mathrm{scat}}\omega_{\perp}$ |  $2\times 10^{-38}$ Jm
$K$ | Luttinger parameter | 25
$T$ | temperature | 2 - 20 nK
$J$ |  $J$-quench | 0 - 30 Hz
$N_{L}$ | number of grid points | 50
$c$ | speed of sound |  $3\times 10^{-3}$ m s-1
a | grid spacing |  $0.36\ \mu$m.
$\xi_{h}$ | healing length |  $0.24\ \mu$m
$\lambda_{T}$ | phase coherence length |  $38-380\ \mu$m
$\xi_{s}$ | spin healing length |  $2.5\ \mu$m
Table 1: Table containing important parameters and their values. The
parameters are chosen to be experimentally feasible and correspond roughly to
those reported in references [72, 75, 74, 73, 76, 77].

###  IV.3 Choice of parameters

There are three constraints which must be satisfied in order to have a quasi-
one dimensional condensate [55]. To ensure minimal scattering into the
transverse modes we need the interaction to be sufficiently weak which implies
$\mu=g_{1D}n_{1D}\ll\hbar\omega_{\perp}$ where $\mu$ is the chemical potential
and $\omega_{\perp}$ is the transverse trapping frequency. Moreover, the
temperature needs to be low enough such that transverse modes are not
thermally excited leading to the inequality $k_{B}T\ll\hbar\omega_{\perp}$.
Finally, in order to have a quasi-condensate which permits a semiclassical
approach we need weak interactions in comparison to the zero-point kinetic
energy associated with the density of the particles. This implies
$n_{1D}g_{1D}\ll\hbar^{2}n_{1D}^{2}/m$ which means the Luttinger parameter
should obey $K\gg 1$. All the parameter values we use satisfy these three
inequalities.

In quasi-one dimensional gases the interatomic interaction parameter $g_{1D}$
is related to the scattering length $a_{s}$ and transverse trapping frequency
as $g_{1D}=2\hbar a_{s}\omega_{\perp}$. For 87Rb atoms we have $a_{s}\approx
98\times 0.52\,$Å[113] and we will assume $\omega_{\perp}=2\,\pi\times$3 kHz
[77]. The full list of parameters used in our simulations is given in Table 1
and roughly corresponds to those used in the experiments by the Vienna group
[72, 75, 74, 73, 76, 77].

For our numerical simulations we choose a grid size that slightly exceeds the
healing length because, as explained above, this cuts off unphysical density
fluctuations [51, 95]. This condition is given in Eq. (12) but can be
expressed succinctly in terms of $\Gamma$ as $N_{L}^{2}<\Gamma N^{2}$. The
magnitudes of $\tilde{\rho}$ and $\tilde{\phi}$ also need to be considered.
The phase difference can take the full range $+\pi$ to $-\pi$, but the number
difference is limited by the condition that the total number difference
(integrated over the entire system) cannot exceed the total number of
particles. In fact, due to the random nature of sampled thermal fluctuations,
the integral of $\tilde{\rho}$ is always approximately zero. However, the
validity of the SG/SG+ model requires that local density fluctuations be small
in comparison to the background density $n_{1D}$, see Appendix A. Translated
into the scaled variables this means that at any point
$\tilde{\rho}(\tilde{z})\ll n_{1D}\xi_{h}$. In practice we choose
$\tilde{\rho}(\tilde{z})\leq 1.6$ so that the fluctuations are an order of
magnitude smaller than the background density.

###  IV.4 Examples of Initial conditions

In Figure 2 we present typical spatial profiles of the initial number
difference field $\tilde{\rho}$ (upper row) and phase difference field
$\tilde{\phi}$ (lower row). Each profile provides the initial conditions for a
single classical field trajectory and is obtained by summing up thermally
activated phonons (Fourier modes) using the Tomonaga-Luttinger model.

|  |
---|---|---
|  |
Figure 2:  Examples of initial spatial profiles of the number difference
$\tilde{\rho}$ (top row) and phase difference $\tilde{\phi}$ (bottom row).
Each profile is obtained by randomly sampling a thermal distribution using the
method described in Section IV.2, and each panel includes ten different
profiles. The parameter values common to all panels include the number of
computational lattice points $N_{L}=50$, grid spacing $a=0.36\mu m$, and
healing length $\xi_{h}=0.24\,\mu$m (the remaining parameters are listed in
Table 1). The difference between the columns is as follows. The left column
has the Luttinger parameter $K=25$, and temperature $T=2$ nK giving a phase
coherence length of $\lambda_{T}=380\,\mu$m. In the middle column $K=25$, but
the temperature is increased to $20\,$nK, giving $\lambda_{T}=38\,\mu$m. In
the right column, the value of $K$ is artificially increased (without changing
any other parameters) to $K=250$ and $T=2$ nK. Increases in temperature excite
stronger fluctuations in the profiles as expected. Increases in the Luttinger
parameter have opposite effects on $\tilde{\rho}$ and $\tilde{\phi}$. The
maximum value and jaggedness of $\tilde{\rho}$ is increased whereas the
jaggedness of $\tilde{\phi}$ is reduced. An explanation of this behavior is
given in the main text.

The different columns show the effect of changing temperature $T$ or Luttinger
parameter $K$. As expected, when $T$ is increased the fluctuations in both
$\tilde{\rho}$ and $\tilde{\phi}$ increase. By contrast, if $K$ is increased
the maximum magnitude and jaggedness of $\tilde{\rho}$ increases but the
jaggedness of $\tilde{\phi}$ decreases. Referring to Eq. (19) we can see that
this is because the coefficient multiplying the density fluctuation term in
the Hamiltonian is $\Gamma=\pi/2K$ which decreases as $K$ increases leading to
increased variance of $\varrho_{k}$ modes according to Eq. (34). The phase
fluctuation term shows the opposite behavior because its coefficient in the
Hamiltonian (which only appears as the spatial gradient of $\tilde{\phi}$) is
$\epsilon=K/2\pi$ which increases as $K$ increases and this reduces the
variance of the $\varphi_{k}$ modes according to Eq. (35), thereby making the
$\tilde{\phi}$ profiles smoother.

##  V Numerical Simulations of the Dynamics

In this section we explore the dynamics following a $J$-quench. Our approach
is inspired by the TWA where multiple classical field configurations are
propagated in time using the classical equations of motion, although in our
case the initial conditions are sampled from a thermal distribution as
described in Section IV rather than a quantum distribution as in the standard
TWA.

$J$-quench dynamics have previously been explored for the simpler case of a
two-mode zero temperature bosonic Josephson junction where it was found that
caustics dominate the number and phase difference probability distributions
[31, 17, 26]. In the two-mode case it is possible to compute the exact quantum
dynamics for some thousands of particles and compare them against the TWA. The
results (see Figure 1 in [31]) show good qualitative agreement and give us
confidence that the TWA can capture the main features of the quantum dynamics.
Furthermore, the inevitable presence of decoherence due to the environment
will tend to reduce the quantum dynamics to their classical limit (this has
been investigated in the two-mode case for a $J$-quench in [32]) increasing
the relevance of semiclassical calculations. In the present work we are
interested in whether the phonons along the long axis disrupt or sustain these
caustics. We will start by reproducing the caustics presented in Ref. 31 for
the two-mode case and then add in the longitudinal modes after that.

###  V.1 Numerical Methods

The initial conditions are generated via random sampling from Gaussian
distributions. We then evolve the equations of motion (Eq. 23 for the case of
the full SG+ model) using a Runge-Kutta solver with a user-defined time step
[114]. The endpoints of our system are treated by imposing periodic boundary
conditions. In Appendix C we demonstrate the numerical convergence of the
solver by varying the temporal and spatial steps by tracking the time
evolution of the total energy (hamiltonian) which should be a constant of the
motion and obtain the fiducial time and space resolution for all our
calculations.

---
Figure 3:  Dynamics of the number difference $\tilde{\rho}$ (left), phase
difference $\tilde{\phi}$ (middle), and phase space distribution (right)
following a $J$-quench from $J=0$ to $J=30$ Hz in the two mode approximation
governed by the Josephson equations given in Eq. (39). The other parameter
values are given in Table 1. Each panel contains 150 trajectories: each
trajectory starts with $\tilde{\rho}=0$ at time $\tilde{t}=0$ but has an
initial phase randomly sampled from $[-\pi,\pi)$. Both number and phase
difference variables display a series of cusp shaped caustics given by the
envelopes of families of trajectories; to guide the eye we have outlined the
first cusp caustic in the $\tilde{\rho}$ variable with a black curve. In the
right panel three different time slices of the results are plotted in phase
space ($\tilde{\rho}$ versus $\tilde{\phi}$). Each dot corresponds to a
different initial condition (trajectory) and the colors indicate the time:
$\tilde{t}$=0 (red), $\tilde{t}$=50 (green), $\tilde{t}$=100 (blue). During
time evolution the initial horizontal line winds into a whorl and the caustics
in the $\tilde{\rho}$ and $\tilde{\phi}$ plots occur due to horizontal and
vertical segments of a whorl, respectively.

###  V.2 Special case: two-mode approximation

In the two-mode approximation only a single mode in each well is taken into
account. This description is relevant to the SG/SG+ model in the limit where
the entire length of each quasicondensate is perfectly synchronized so that
the fields $\tilde{\rho}(\tilde{z})$ and $\tilde{\phi}(\tilde{z})$ do not
depend on $\tilde{z}$. In this case the spatial derivative terms vanish and
the equations of motion in Eq. (23) reduce to

$\frac{d\tilde{\phi}}{d\tilde{t}}=2\Gamma\tilde{\rho}\quad,\quad\frac{d\tilde{\rho}}{d\tilde{t}}=-2\mathcal{J}\sin\tilde{\phi}\
.$ (39)

These are the standard Josephson equations of motion and also correspond to
those of a classical pendulum [115]. Such synchronization can occur at very
low temperatures or when the coefficients $\epsilon$ and $\Gamma$ are large
enough that they suppress spatial fluctuations in the initial conditions.

In Figure 3 we display the post-quench dynamics in the two-mode approximation.
The left hand and central panels show the time dependence of $150$ independent
solutions of Eq. (39) which give the trajectories for the number difference
and phase difference, respectively. Note that in this paper we use the color
blue for trajectories calculated within the two mode approximation and reserve
red for the trajectories of the full many mode model. In accordance with our
assumption that the two wells start with an equal number of atoms, each
solution starts with $\tilde{\rho}=0$. And as discussed in Section IV.2, the
initial value of $\tilde{\phi}$ is randomly chosen from the range $[-\pi,\pi)$
because the two condensates are independent before the $J$-quench.

The most striking feature of Figure 3 is the series of cusp-shaped caustics
that form in both variables. In order to guide eye, we have have outlined the
first cusp caustic in the number difference variable using a black curve (the
calculation for this curve is given in Appendix D). Like in optics, caustics
are regions of high intensity formed by the envelopes of families of rays
(trajectories). Each caustic is born at the centre of the distribution at the
tip of a cusp before spreading out in two arms that move towards the edges of
the distribution. The fact that they are cusp shaped is in agreement with the
prediction of catastrophe theory that in two dimensions the only structurally
stable and hence generic singularities are cusps.

Each trajectory represents a single experimental run. The idea behind the TWA
is that the number of trajectories reaching a point $\tilde{\rho}$ at time
$\tilde{t}$ is proportional to the probability that a measurement of the true
quantum system would yield that value of $\tilde{\rho}$. An equivalent
interpretation holds for the $\tilde{\phi}$ trajectories. The caustics have
the highest probability density and hence give the values most likely to be
observed. Of course, if we only consider the average values of $\tilde{\rho}$
or $\tilde{\phi}$ we would get zero in both cases due to the symmetry of the
distributions and hence miss the caustics. Many experimental runs must be
performed in order to obtain the probability distribution where these patterns
live.

Figure 4:  Dynamics of the number difference $\tilde{\rho}$ (left), phase
difference $\tilde{\phi}$ (middle), and the phase space distribution (right)
in the linearized version of the two-mode approximation [Eq. (40)] following a
$J$-quench from $J=0$ to $J=30$ Hz. Like in Figure 3, there are 150
trajectories shown in each panel corresponding to different values of the
initial value of $\tilde{\phi}$. However, in this linearized case we obtain a
series of perfect focus points (revivals of the initial state). This is
because linearization gives rise to rigid rotation in phase space without
whorls. Unlike the extended cusp caustics seen in Figure 3 (which will be
qualitatively robust to details of the nonlinearity) perfect focus points are
nongeneric because they are unstable to perturbations such as the effects of
nonlinearity. All parameter values and color labels are the same as Figure 3.

The mechanism underlying caustics can be understood from a phase space
perspective, as shown in the right hand panel of Figure 3. Each dot gives the
number and phase difference at a particular time for a different initial
condition. The red dots are the initial values which lie in a horizontal line
because at $\tilde{t}=0$ all trajectories have $\tilde{\rho}=0$. As time
evolves the dots rotate around the origin: the green and blue dots show two
successively later times. However, the nonlinearity of the Josephson equations
means dots further from the origin rotate more slowly and this leads to the
formation of a spiral or whorl. At places where the whorl has a vertical
segment a range of different solutions all have the same value of
$\tilde{\phi}$ and this stationarity of the distribution with respect to
changes in the initial conditions is what generates a caustic, in this case a
$\tilde{\phi}$-caustic. Conversely, horizontal segments give rise to
$\tilde{\rho}$-caustics.

In the absence of nonlinearity the equations reduce to those of a harmonic
oscillator

$\frac{d\tilde{\phi}}{d\tilde{t}}=2\Gamma\tilde{\rho}\quad,\quad\frac{d\tilde{\rho}}{d\tilde{t}}=-2\mathcal{J}\tilde{\phi}$
(40)

giving rise to rigid rotation in phase space and the formation of perfect
focal points in the number and phase difference variables, as shown in Figure
4. However, these perfect revivals of the initial state are not stable: any
nonlinearity will cause the focal points to evolve into the extended cusp
caustics shown in Figure 3.

The frequency of the linearized motion is known in Josephson junction
terminology as the plasma frequency. In our notation it reads

$\omega_{p}=\sqrt{4\Gamma\mathcal{J}}$ (41)

and the period of the motion is therefore given by $2\pi/\omega_{p}$. For the
case shown in Figure 4 we have $\Gamma=0.063$ and $\mathcal{J}=0.037$ giving a
period $\approx 65$. In fact, the tips of the cusps in the nonlinear case also
occur with this period since they are formed from small amplitude trajectories
that only experience the quadratic bottom of the cosine potential.

---
Figure 5:  Dynamics of the spatially averaged number difference
$\langle\tilde{\rho}\rangle_{z}$ (left), phase difference
$\langle\tilde{\phi}\rangle_{z}$ (middle), and phase space distribution
(right) for the full many-mode SG+ model following a $J$-quench from $J=0$ to
$J=30$ Hz. Each panel contains 150 trajectories which are solutions of Eq.
(23). The initial conditions are randomly sampled thermal phonons with the
same parameter values as those shown in the top row of Figure 2 and described
in Table 1. In particular, the number of numerical lattice points is
$N_{L}=50$ separated by a grid spacing of $a=0.36\ \mu$m, and the temperature
is $T=2$ nK. The healing length is $\xi_{h}=0.24\ \mu$m, the spin healing
length is $\xi_{s}=2.5\ \mu$m and the phase coherence length is
$\lambda_{T}=380\ \mu$m. The different colors on the phase space plot
correspond to the same time slices as in the previous phase space plots.

###  V.3 General case: many-mode SG+ model

Simulations of the full SG+ model are shown in Figure 5, which represents one
of the main results of this paper. The trajectories in the left panel give the
spatially averaged number difference
$\langle\tilde{\rho}(\tilde{t})\rangle_{z}$ as a function of time obtained by
solving the equations of motion given in Eq. (23) for the many-mode system and
then averaging over its length. The trajectories in the middle panel of Figure
5 give the equivalent spatial average of the phase difference
$\langle\tilde{\phi}(\tilde{t})\rangle_{z}$, and the right-hand panel is the
phase space picture. Each trajectory is evolved from a single randomly sampled
field configuration (describing thermally activated phonons) such as those
shown in the top row of Figure 2 and for the parameters given in Table 1. We
observe that despite the inclusion of longitudinal modes and the randomness of
the initial conditions, the caustics survive and are quite similar to those of
the two-mode approximation shown in Figure 3. This suggests that caustics are
a generic feature of many particle dynamics following quenches, at least for
systems whose underlying physics is based on coupled nonlinear oscillators.
Each oscillator starts with a random phase and a noisy momentum but the quench
acts so as to give all the oscillators a momentum kick at the same time
$\tilde{t}=0$ leading to an initial partial synchronization. As the system
evolves in time after the kick the different periods of nonlinear oscillators
leads to cusp catastrophes in the distribution of trajectories. If we had
instead calculated only the expectation values of the number and phase
differences then this underlying structure would not have been visible because
it lives in the probability distribution rather than the mean values.

A slice at fixed time through the probability distribution for the spatially
averaged phase variable $\langle\tilde{\phi}\rangle_{z}$ is shown in Figure 6.
This is obtained by sorting the $\langle\tilde{\phi}\rangle_{z}(\tilde{t})$
trajectories into bins each of which covers a small range of
$\langle\tilde{\phi}\rangle_{z}$ and counting the number of trajectories in
each bin. The result is noisy due to the thermal fluctuations but the caustics
are clearly visible as strong peaks. These peaks display the characteristic
‘square root’ divergence of fold caustics [1]

$\mathcal{P}(\langle\tilde{\phi}\rangle_{z})\propto\frac{1}{\sqrt{\tilde{\phi}_{c}-\langle\tilde{\phi}\rangle_{z}}}$
(42)

where $\mathcal{P}(\langle\tilde{\phi}\rangle_{z})$ is the probability density
and $\tilde{\phi}_{c}$ is the location of the caustic. The blue dashed lines
in Figure 6 are fits of Eq. (42) to the numerical data and we see that the
agreement is good. Although the height of the singularities predicted by Eq.
(42) is infinite at the caustic, this function is integrable so that a
probability distribution with caustics is still normalizable (of course, the
peaks in the numerical data are of finite height because the number of
trajectories is finite). A very similar pattern of square root singularities
at each caustic is obtained for a time slice through the probability density
for the number difference variable so we shall not show it here.

---
Figure 6:  The probability density (red curve) as a function of
$\langle\tilde{\phi}\rangle_{z}$ obtained from the density of trajectories at
time $\tilde{t}=162$ for the SG+ model. This corresponds to a slice at fixed
time through the middle panel of Figure 5, although calculated using 10000
trajectories to improve the statistics and averaged over a short time window
of $\Delta\tilde{t}=1$ to remove rapid time fluctuations. The red curve has
been drawn with a bin width $d\tilde{\phi}=0.04$ and is normalised such that
the area under the graph is 1. The caustics are clearly visible as diverging
peaks and are well fitted (blue dashed curves) by the inverse square root form
given in Eq. (42) that is expected for fold catastrophes [1] (the satellite
caustics also have this shape but the fit is not shown to avoid obscuring the
data). A very similar profile is obtained for the probability density in the
$\langle\tilde{\rho}\rangle_{z}$ variable (not shown).

###  V.4 Effect of dispersion on the caustics

The double derivative terms in the SG+ equations of motion given in Eq. 23 are
responsible for transmitting wave disturbances along the longitudinal axis and
are not present in the simpler two-mode case discussed in Section V.2. Initial
thermal fluctuations in the SG+ model will therefore disperse in $z$ over time
and it is interesting to see what difference this makes to the caustics;
comparison of Figures 3 and 5 suggests it makes little difference to spatially
averaged variables. However, this observation is for only one choice of the
parameters $\epsilon$ and $\Gamma$ that govern the size of the derivative
terms and also for relatively short times. In particular, in Figure 5 the
parameters are $\epsilon\approx 4$ and $\Gamma\approx 0.06$ which were chosen
to match experimental values [72, 75, 74, 73, 76, 77]. In Figure 7 we compare
the long time dynamics of the two-mode approximation and the SG+ model for the
case where $\epsilon$ in the SG+ model has been artificially increased by a
factor of 10 (without changing any other parameters), thereby increasing the
effect of spatial dispersion. Apart from this change, the initial conditions
and $J$-quench are similar to those used in Figure 5. Note that we only use
this increased value of $\epsilon$ for the time propagation and not for the
generation of the thermal initial conditions. This avoids changing the
starting phase fluctuations from those used in Figure 5 which would otherwise
be energetically suppressed and would also lead to significantly different
dynamics but is not the comparison we would like to make here. From Figure 7
we see that the strong coupling of neighboring ‘pendula’ does seem to largely
wash out the caustics at long times in comparison to the dispersionless two-
mode case, although some faint structure is still present which underlines the
structural stability of caustics. The long time behavior will be further
analyzed in Section VI.

---
Figure 7:  Comparison of the long-time behavior of the phase difference in the
two-mode approximation (upper) and many-mode SG+ model (lower). Both panels
contain 150 different runs and the initial conditions and $J$-quench are
similar to those of Figure 5 except that $\epsilon$ has been artificially
multiplied by 10 (without changing any other parameters) in the lower panel.
This enhances the effect of the spatial derivative term in $\phi$ in the SG+
model (this term does not appear in the two mode model). We see that in the
upper panel the caustics are still visible. By contrast, the stronger spatial
interaction causes dispersion and makes the caustics much less visible in the
lower panel.

###  V.5 Effect of $J$ on the caustics

Another parameter that affects the dynamics is the tunnel coupling strength
$J$ [or its dimensionless version $\mathcal{J}$ which is defined in Eq. (19)]
that becomes non-zero after the quench. The quench itself creates a strongly
nonequilibrium phase difference where all values of $\tilde{\phi}$ are equally
probable independently of the value of $J$ by virtue of the fact that before
the quench there is no phase coherence between the two quasicondensates.
However, $J$ does control the post-quench dynamics. One way it does this is
via the frequency of the Josephson oscillations. The cusps occur with a
frequency given by the plasma frequency in Eq. (41) which goes as $\sqrt{J}$.

---
Figure 8:  Effect of quench strength $J$ for $J=0$ Hz, 3 Hz, and 30 Hz (from
left to right). The top row shows the dynamics of
$\langle\tilde{\rho}\rangle_{z}$ with initial conditions sampled in the same
way as in Figure 5. The bottom row plots the corresponding phase space
distributions. Like in previous figures, the different colors give different
time instants: $\tilde{t}$=0 (red), $\tilde{t}$=50 (green), $\tilde{t}$=100
(blue). The dots with intense colors are the spatially averaged values. We
have also included the raw data (without spatial averaging) as faint dots.
This gives an idea of the size of the statistical fluctuations due to the
thermal initial conditions and is the same for all values of $J$. In the left
column there is no coupling between the two quasicondensates and hence no time
evolution of the spatially averaged data (the intense red, green, and blue
dots sit on top of each other) although there can be evolution of unaveraged
data due to intrawell dynamics, i.e. without the $J$ term in Eq. (10). As we
increase the magnitude of $J$ time evolution leads to whorls with a greater
vertical extent because more energy can be extracted from the cosine potential
in Eq. (18) giving larger values of
$\langle\tilde{\rho}\rangle_{z}^{\mathrm{max}}$.

In Figure 8 we examine the effect of quenching to different $J$ values, with
the value of $J$ increasing from left to right. We can see the expected
increase in frequency. The amplitude of the motion also increases with $J$
because immediately after the quench each trajectory finds itself at a random
point on the cosine potential energy surface whose depth between valley top
and valley bottom is $2\mathcal{J}$. The initial potential energy of a field
configuration is therefore
$-2\mathcal{J}\langle\cos\tilde{\phi}_{0}\rangle_{z}$, where
$\tilde{\phi}_{0}$ is the phase field $\tilde{\phi}(\tilde{z},\tilde{t})$ at
the initial time. This configuration evolves under the full Hamiltonian and
upon spatial averaging is seen to execute oscillations about the potential
minimum. The upper row in Figure 8 plots the spatially averaged number
difference and according to Eq. (18) the maximum amplitude this can have is

$\langle\tilde{\rho}\rangle_{z}^{\mathrm{max}}=\sqrt{\frac{2\mathcal{J}(1-\langle\cos\tilde{\phi}_{0}\rangle_{z})}{\Gamma}}$
(43)

where we have ignored the effects of spatial coupling (second order derivative
terms). Thus, $\langle\tilde{\rho}\rangle_{z}^{\mathrm{max}}$ also scales as
$\sqrt{J}$, and this is in correspondence with Figure 8.

The lower row of Figure 8 shows the behavior in phase space. In these figures
we have also included the unaveraged data, i.e. the $\tilde{\rho}$ and
$\tilde{\phi}$ values of each grid point at the three selected times. This
gives a sense of the size of the statistical fluctuations due to the spatial
degrees of freedom. In the left hand column $J$ remains zero for all time and
the only dynamics that can occur is along the long-axis of each
quasicondensate individually. The middle and right hand panels, which have
$J=3$ and $J=30$ Hz, respectively, have the same initial statistical
fluctuations as the left hand one because, as mentioned above, the initial
distribution is set by the pre-quench thermal fluctuations in the two
quasicondensates and is independent of $J$. However, as time evolves the
effects of $J$ described by Eq. (43) become apparent because larger $J$ allows
a greater value of $\langle\tilde{\rho}\rangle_{z}^{\mathrm{max}}$ and this
stretches the distribution along the vertical direction in comparison to a
smaller value of $J$. For a whorl to become apparent
$\langle\tilde{\rho}\rangle_{z}^{\mathrm{max}}$ should at least exceed the
width of the statistical fluctuations and becomes better and better defined as
$J$ is increased.

##  VI Universality and caustics

We have already discussed the relationship between nonlinearity and caustics
in the preceding section. As motivated earlier, and expounded in Refs. 31, 17,
33, 26, 34, caustics also have implications for the universal dynamics of
quantum systems. We explore a few of these effects in this section.

###  VI.1 Long time distribution: the circus tent

The quench generates collective excitations that lead to caustics as shown in
Figures 3 and 5 for the two nonlinear models (two mode and SG+) discussed
above. The caustics are born at the center of the probability distribution (in
either the $\tilde{\rho}$ or the $\tilde{\phi}$ variable) at intervals of the
plasma period and move out to the edges over time. Figure 6 plots the
probability distribution for the SG+ model as a function of
$\langle\tilde{\phi}\rangle_{z}$ at an intermediate time where four pairs of
fold caustics are discernible and shows how they diminish in strength but are
still present as they move to the edges. The question then naturally arises as
to what happens at long times $\tilde{t}\rightarrow\infty$ when the
distribution comprises of a large number of caustics and whether it tends to a
characteristic shape? The answer is yes, and is shown in Figure 9 which is
made in the same way as Figure 6 but this time by calculating the density of
$\langle\tilde{\rho}\rangle_{z}$ trajectories and averaging over a time window
extending between $\tilde{t}=800$ and $\tilde{t}=980$ in order to remove rapid
fluctuations. The probability distribution takes a shape reminiscent of a
‘circus tent’ or ‘big top’ and can be understood as follows. The strongest
singularities present are the cusp tips born at the center of the distribution
which leads to this being the highest point. Each cusp then splits into two
fold arms (which according to catastrophe theory are lower singularities) that
move outwards, reducing in height as they go, before accumulating at the edges
where there is a sharp drop to zero. The position of the outer edge is set by
the maximum energy that can be extracted from the quench and is given by Eq.
43.

An analytic expression for the circus tent distribution is given by the
integral

$\mathcal{P}_{\mathrm{CT}}(\tilde{\rho})=\frac{1}{2\pi
B}\int_{\tilde{\rho}^{2}/B^{2}}^{1}\frac{U(m,\tilde{\rho})}{K(m)}\ dm$ (44)

where

$U(m,\tilde{\rho})=\frac{1}{\sqrt{m(1-m)(m-\tilde{\rho}^{2}/B^{2})(1+\tilde{\rho}^{2}/B^{2}-m)}},$
(45)

$K(m)$ is the complete elliptic integral of the first kind, and
$B=2\sqrt{\mathcal{J}/\Gamma}$. This expression is plotted in Figure 9 as the
dashed line and is derived in Appendix E under the assumption that at long
times we can model the system by an ensemble of independent pendulua where
each pendulum is ergodic. In other words, each pendulum obeys a microcanonical
distribution where there is equal probability for it to be found anywhere on
its energy shell. The nature of the $J$-quench is such that it leads to an
ensemble with an equal probability for any starting angle (this is different
to an equal probability for each energy due to the dependence of the density
of states on angle). As can be seen from Figure 9,
$\mathcal{P}_{\mathrm{CT}}(\tilde{\rho})$ gives a good fit to the numerical
data generated by both the SG+ and two-mode models considered in this paper.

---
Figure 9:  The long time probability distribution for the number difference
$\tilde{\rho}$. The data points are from the different nonlinear models
considered in this paper averaged over the spatial coordinate $z$ and also
over a time window ranging from $\tilde{t}=800$ to $\tilde{t}=980$ to remove
fluctuations. The pink dashed line is the circus tent distribution
$\mathcal{P}_{\mathrm{CT}}$ given in Eq. (44) and derived in Appendix E under
the assumption of ergodicity; the circus tent shape is due to the
proliferation of caustics at long times and gives a good fit to the data. The
solid black curve is the thermal distribution $\mathcal{P}_{T}$ with a
temperature chosen so that the expectation value of the energy matches that
provided by the quench.
---
Figure 10:  Structural stability of caustics: here we investigate the effect
of unbalanced densities on caustics by tracking the same SG+ model dynamics as
those shown in Figure 5 except for an initial density imbalance of 0.1 in the
background of $\tilde{\rho}$ at each point z. We see that the cusp caustics in
the plots of $\langle\tilde{\rho}\rangle_{z}$ and
$\langle\tilde{\phi}\rangle_{z}$ versus time are distorted but still maintain
their basic structure. This is because the whorl in phase space is left intact
despite having a displaced centre. Caustics are resilient against
imperfections and perturbations and we expect them to be present under
realistic experimental conditions.

In Figure 9 we also include the thermal probability distribution

$\mathcal{P}_{T}(\tilde{\rho})=\frac{1}{Z}\int_{0}^{\infty}\mathcal{P}_{E}(\tilde{\rho})\
e^{-E/T}D(E)\ dE$ (46)

describing an ensemble of pendula at thermal equilibrium at temperature $T$
where $\mathcal{P}_{E}(\tilde{\rho})$ is the probability distribution at fixed
energy $E$, $D(E)$ is the density of states and $Z$ is a normalizing factor.
The details of our calculation of $\mathcal{P}_{T}(\tilde{\rho})$ are given in
Appendix F, where, for example, $\mathcal{P}_{E}(\tilde{\rho})$ is given in
Eq. (82). The temperature of this distribution is chosen such that the mean
energy of the thermal distribution $\langle E\rangle_{T}$ is equal to the mean
energy of the states excited by the quench. For a quench to $J=30$ Hz we show
in Appendix F that the effective temperature is 5.4 nK.

Clearly, the thermal distribution is very different to the circus tent
distribution: the thermal distribution takes the form of a smooth gaussian
with wings that extend beyond $\langle\tilde{\rho}\rangle_{z}^{\mathrm{max}}$
because the thermal Boltzmann factor allows for excitations with any energy
(albeit with exponentially small probability) including those involving
pendula undergoing rotation as well as libration, whereas the $J$-quench only
excites librational motion. The probability distribution for a thermal
pendulum is in fact quite delicate to compute because of the singularity in
the density of states between libration and rotation but the combined result
is smooth; see Appendix F for more details.

###  VI.2 Structural stability of caustics

The defining characteristic of the singularities described by catastrophe
theory is structural stability against perturbations and this ensures that
they occur generically. The same is not true of isolated singularities as can
be seen by comparing Figures 3 and 4 where it is shown that point foci do not
survive the introduction of nonlinearity. In two dimensions cusps are the
unique structurally stable catastrophe and from Figures 3 and 5 we see that
cusp-shaped caustics are indeed stable against random thermal fluctuations.
However, thus far we have imposed the symmetrical starting condition that the
initial number difference between the two quasicondensates is zero. One may
therefore wonder whether the caustics we see are a consequence of this
symmetry. To check that this is not the case we show in Figure 10 the dynamics
for the case where the initial background density $n_{1D}$ in the two
quasicondensates differs by $10\%$. We see that although the caustics in both
$\langle\tilde{\rho}\rangle_{z}$ and $\langle\tilde{\phi}\rangle_{z}$ are
distorted they maintain their basic cusp shape. Furthermore, the phase space
whorls still occur and this guarantees the existence of caustics.

###  VI.3 Coherence factor and relaxation towards equilibrium

Cold atom experiments have the ability to measure correlation functions in
nonequilibrium many-body states [116, 74, 117, 118]. As a simple example let
us consider the coherence factor

$\begin{split}\mathds{C}(\tilde{t})=\left\langle\langle\cos\tilde{\phi}\rangle_{z}\right\rangle\end{split}$
(47)

which depends on the spatial average of the phase difference field
$\tilde{\phi}(\tilde{z},\tilde{t})$ between points along the two
quasicondensates. The outer brackets indicate an ensemble average which means
averaging over many trajectories each sampled from the thermal distribution
discussed in Sec. IV. In the Vienna experiments, where one quasicondensate is
suddenly split into two, the coherence starts near unity and decays over time
as the two quasicondensates decohere [76, 77]. In the opposite case, where two
independent quasicondensates are suddenly coupled, one expects the converse
where the coherence starts at zero and grows. This situation has been
previously modelled by Horváth et al. using both the TWA and a truncated
conformal space approach [95]. They found that $\mathds{C}(\tilde{t})$
initially grows and then undergoes damped oscillations as it settles down
towards a finite constant value. The coherence factor therefore provides a
measure of how the system reaches equilibrium. In this context we note that
$\mathds{C}(\tilde{t})$ actually corresponds to an ensemble average of the
cosine term in the SG/SG+ Hamiltonian and thus gives information on the
exchange of energy between the different parts. In other words, since the
total energy is a constant of the motion, if the ‘potential’ part of the
energy settles down to a constant this suggests the ‘kinetic’ parts of the
energy are also constant, at least from an ensemble averaged point of view.
Our aim in this section is to see if the dynamics of $\mathds{C}(\tilde{t})$
is connected to the caustics.

In Figure 11 we plot $\mathds{C}(\tilde{t})$ for two models: the full SG+
model which is many-mode and nonlinear and a linearized version which obeys
the equations of motion

$\begin{split}\frac{d\tilde{\phi}}{d\tilde{t}}&=2\Gamma\tilde{\rho}-\frac{\Gamma}{2}\frac{\partial^{2}\tilde{\rho}}{\partial\tilde{z}^{2}}\\\
\frac{d\tilde{\rho}}{d\tilde{t}}&=2\epsilon\frac{\partial^{2}\tilde{\phi}}{\partial\tilde{z}^{2}}-2\mathcal{J}{\tilde{\phi}}.\end{split}$
(48)

This differs from the linearized two-mode approximation defined by Eq. (40)
because it describes an elongated multi-mode system. From Figure 11 we see
that $\mathds{C}(\tilde{t})$ for the SG+ model (dark blue curve) does indeed
initially grow, undergo damped oscillations and settle down to a non-zero
value (the fact that $\mathds{C}(\tilde{t})\neq 0$ at $\tilde{t}=0$ is due to
random fluctuations in the initial conditions: as we include more trajectories
we find that the initial value gets smaller). Meanwhile,
$\mathds{C}(\tilde{t})$ for the linear model (red dashed curve) executes
undamped oscillations and hence does not settle down to equilibrium. Both
models agree during the first oscillation but strongly differ after that.

It is clear that nonlinearity is important for reaching equilibrium at least
as far as global quantities such as $\mathds{C}(\tilde{t})$ are concerned. We
can understand this by interpreting the SG+ model as describing a chain of
coupled pendula. The nonlinearity of each pendulum means that its period
depends on the amplitude of its motion and hence an ensemble of pendula whose
motion is initiated together by the quench, but all with different degrees of
excitation, will dephase from one another over time so that collective
oscillations are damped out. By contrast, linear oscillators have a period
independent of their amplitudes of motion and hence remain in phase.

Apart from the ensemble averages shown by the darker curves in Figure 11, we
have also included the individual trajectories for
$\langle\cos\tilde{\phi}\rangle_{z}$ as fainter curves. The linear model
displays harmonic motion and hence perfect revivals whereas the trajectories
in the nonlinear model give rise to half-cusp caustics. These caustics overlap
in time such that averaging over them causes the coherence to strongly relax
after a single period. It is not so much that the caustics cause the
relaxation, but rather that both have a common origin in the nonlinearity of
the model and hence are generic features of dynamics in complex systems.

---
Figure 11:  The two dark lines give the time evolution of the coherence factor
$\mathds{C}(\tilde{t})$ defined in Eq. (47) for a linear model (dashed-dotted
red) and the SG+ model (solid blue). Both models are multi-mode (many
longitudinal modes along $\tilde{z}$) but the SG+ model is nonlinear. Also
included as faint lines are the raw trajectories
$\langle\cos\tilde{\phi}\rangle_{z}$ from which $\mathds{C}(\tilde{t})$ is
composed. As everywhere in this paper, $\langle\ldots\rangle_{z}$ indicates a
spatial average. This figure highlights that recurrences present in the linear
case are suppressed by nonlinearity in the SG+ system. The ensemble average
over trajectories with different periods causes $\mathds{C}(\tilde{t})$ to
relax towards an equilibrium value in the case of the SG+ model in line with
previous experimental observations [77, 76] and theory [95].

##  VII Summary and Conclusions

The sine-Gordon (SG) model is a nonlinear integrable field theory that can be
used to describe a wide range of systems from high energy physics to condensed
matter physics. A series of landmark experiments using two coupled 1D atomic
quasicondensates [63, 71, 72, 73, 74, 75, 76, 77] have realized the SG model
in a controllable quantum many body environment. The key parameters can be
varied in time allowing the implementation of sudden quenches that excite many
modes leading to nonequilibrium dynamics. This is the setting we adopt for the
current paper where we use experimentally realistic parameters and compute the
dynamics of the number and phase difference fields. However, in contrast to
the usual experimental protocol where the tunnel coupling $J$ is suddenly
switched off, we consider quenches where it is suddenly switched on. While the
former case is adapted to studying dephasing, decay and thermalization between
the two subsystems, the many body dynamics is governed by the Tomonaga-
Luttinger Hamiltonian describing independent 1D quasicondensates. If instead
$J$ is suddenly switched on then the dynamics is that of the full SG model.

Our calculations employ a thermal version of the semiclassical truncated
Wigner approximation (TWA) method. More specifically, we propagate a large
number of classical field configurations over time with initial conditions
sampled from a distribution at thermal equilibrium. The time evolved
configurations (trajectories) can be summed to obtain the probability
distributions for the observables and we find that these are dominated by
singular caustic patterns. The natural mathematical description of caustics is
catastrophe theory that predicts a hierarchy of structurally stable
singularities with characteristic shapes that depend on dimension. In two
dimensions (e.g. number or phase difference versus time) the structurally
stable catastrophes are fold lines that meet at cusps. This is exactly what we
find in both the number and phase differences following a $J$-quench, see
Figure 5. The probability distributions develop trains of caustics that are
born periodically as cusp points (located at the center of the distribution if
there is no tilt) at each plasma period and evolve into pairs of fold lines
that gradually move out to the wings where they accumulate. Fold catastrophes
manifest as strong non-gaussian fluctuations in the form of inverse square
root divergences in the intensity (probability density), as shown in Figure 6.

A special case is provided by the dynamics of a two mode system as shown in
Figure 3. Here the equations of motion are the Josephson equations given in
Eq. (39). The only fluctuations we include in this example are the quantum
fluctuations in the initial relative phase between the two condensates as
mandated by the uncertainty principle applied to systems in relative number
eigenstates. The two-mode case is relevant to small systems where the higher
modes are well above the temperature scale and so any spatial fluctuations are
suppressed. By contrast, the many-mode case shown in all the other figures
includes both quantum fluctuations and thermal fluctuations in the
longitudinal modes, i.e. thermal occupation of phonon modes in the 1D
quasicondensates. Despite the presence of the many longitudinal modes
(typically 50 in our calculations, as set by the parameter $N_{L}$) which give
rise to highly random looking phase and density profiles as seen in Figure 2,
we find that number and phase caustics survive for experimentally realistic
parameters. Furthermore, the qualitative features of the caustics are stable
against variations in quench strength and density imbalance, as seen in
Figures 8 and 10, respectively, and also against the details of the model (in
this paper we use the SG+ model which augments the SG model by including
longitudinal density gradients). All of these different examples confirm the
structural stability of caustics which is the reason why they occur
universally without the need for fine tuning.

The proliferation of caustics over time combined with their migration to the
edge of the probability distribution has important consequences for the long
time probability distribution. It takes on the shape of a circus tent
featuring a strong central peak due to the cusp tips which are the most
singular part of a caustic, flatter intermediate regions, and rapidly decaying
edges where the caustics pile up, see Figure 9. This shape is quite distinct
from a gaussian thermal distribution and can be derived assuming an ergodic
hypothesis in which individual pendula have equal probability to be anywhere
on their energy shell (see Appendix E). The approach to this equilibrium
distribution can be tracked over time using the coherence factor (Figure 11)
which is a spatial and ensemble average over the phase field and corresponds
to the cosine term in the Hamiltonian if the latter is ensemble averaged. The
attainment of equilibrium relies on the nonlinearity of the system to dephase
itself when ensemble averaged. The caustics also rely on the nonlinearity
without which they would reduce to nongeneric perfect revivals (point foci).
In this sense caustics are mutually exclusive to recurrences, at least in the
statistical sense in which caustics appear in this paper.

Caustics in the SG model could be observed experimentally by measuring the
probability density for either the phase difference or the number difference.
For example, the phase difference can be obtained by releasing the two
quasicondensates from their double well potential and letting them overlap
[80, 81, 82]. This process must be repeated many times and for as near
identical initial conditions and time evolution as possible in order to build
up a probability distribution, although due to the structural stability of
caustics they will not be particularly sensitive to differences in the
experimental setup from run to run. If the probability distribution is
obtained for a single time then we expect to see something like that shown in
Figure 6. In order to observe the time evolution of a caustic, one must then
repeat the whole process for a range of different evolution times. This is
laborious but technically possible, and since the first cusp caustic appears
at half the plasma period the experiment does not need to run for long.

The singular nature of caustics means that they dominate wave fields and are
well known in hydrodynamics and optics through phenomena such as tsunamis and
gravitational lensing. The results of this paper show that they also occur in
the nonequilibrium dynamics of 1D superfluids where a quench plays an
analogous role to an underwater earthquake by generating strong excitations
beyond the linear regime that are focused in this case by the cosine term in
the SG Hamiltonian. The universal properties of catastrophes imply caustics
likely also occur in the post-quench dynamics of other condensed matter
systems too: systems with more degrees of freedom will display higher
catastrophes beyond folds and cusps such as hyperbolic and elliptic umbilics
[34]. However, a special feature of the SG model is that it is integrable and
so one may ask if that property plays a crucial role in the existence of
caustics. In this context, we note that in classical mechanics caustics are
closely associated with the existence of tori in phase space upon which
trajectories live [58]. Tori are broken up by chaos, and thus caustics are not
expected to survive for long in systems which are deep in the chaotic regime.
Despite this, the Kolmogorov-Arnold-Moser (KAM) theorem shows that some tori
survive in moderately chaotic systems [119], which suggests caustics may also
survive in cases where the classical phase-space is mixed, which is the
typical case. Indeed, they survive in the three site Bose-Hubbard model [34]
which is known to be chaotic [120]. The important problem of extending the KAM
theorem to quantum mechanics [121] is thus intertwined with the analysis of
caustics in quantum systems and provides an interesting direction for
extending the present work.

## Acknowledgements

We thank Ryan Plestid for contributions on thermal field sampling in the early
stages of this project, Josh Hainge for suggesting the term ‘circus tent’, and
Igor Mazets for correspondence and advice about experiments. This work was
supported by the Mitacs Globalink research internship, by the Natural Sciences
and Engineering Research Council of Canada (NSERC), and Research at the
Perimeter Institute is supported in part by the Government of Canada, through
the Department of Innovation, Science and Economic Development Canada, and by
the Province of Ontario, through the Ministry of Colleges and Universities.
M.K. would like to acknowledge support from the project 6004-1 of the Indo-
French Centre for the Promotion of Advanced Research (IFCPAR), Ramanujan
Fellowship (SB/S2/RJN-114/2016), SERB Early Career Research Award
(ECR/2018/002085) and SERB Matrics Grant (MTR/2019/001101) from the Science
and Engineering Research Board (SERB), Department of Science and Technology
(DST), Government of India. M.K. acknowledges support from the Infosys
Foundation International Exchange Program at ICTS. M.K acknowledges support of
the Department of Atomic Energy, Government of India, under Project No.
19P1112R&D.

##  Appendix A Derivation of the sine-Gordon Hamiltonian

In this appendix we derive the Hamiltonian $H_{\mathrm{SG}}$ as the effective
low energy description of two cigar shaped tunnel-coupled quasicondensates
[50, 74] within a classical field description (Gross-Pitaevskii theory). Along
the way we also obtain a slightly enhanced Hamiltonian $H_{\mathrm{SG}+}$ that
includes contributions from the gradient of density fluctuations that are not
included in the sine-Gordon (SG) Hamiltonian. These contributions are not very
important for our parameters but play an important conceptual role by
introducing an energetic price for a rapidly varying density and hence
effectively cut off these fluctuations.

Assuming tight radial trapping such that each quasicondensate is in its radial
ground state, meaning that only longitudinal excitations are taken into
account, the second quantized Hamiltonian for the total system be written

$\begin{split}H&=\int_{-\infty}^{\infty}dz\
\Bigg{\\{}\sum_{j=1,2}\bigg{[}-\frac{\hbar^{2}}{2m}\hat{\psi}_{j}^{\dagger}(z)\frac{\partial^{2}\hat{\psi}_{j}(z)}{\partial
z^{2}}+\\\
&U(z)\hat{\psi}_{j}^{\dagger}(z)\hat{\psi}_{j}(z)+\frac{g_{1D}}{2}\hat{\psi}_{j}^{\dagger}(z)\hat{\psi}_{j}^{\dagger}(z)\hat{\psi}_{j}(z)\hat{\psi}_{j}(z)\bigg{]}\\\
&-\hbar
J\Big{[}\hat{\psi}_{1}^{\dagger}(z)\hat{\psi}_{2}(z)+\hat{\psi}_{2}^{\dagger}(z)\hat{\psi}_{1}(z)\Big{]}\Bigg{\\}}\
.\end{split}$ (49)

The quantum field operator $\hat{\psi}_{j}(z)$ annihilates a particle at the
point $z$ in the $j^{\mathrm{th}}$ well, where $z$ is the coordinate along the
longitudinal direction (long axis of the system). $m$ is the mass of the
particles, $U(z)$ is a possible external potential (in this paper it will be
set to zero), $g_{1D}$ controls the interparticle interaction strength, and
$J$ is the tunneling frequency between the two wells. In the classical field
approximation we replace the field operators by complex functions

$\hat{\psi}_{j}(z)\rightarrow\psi_{j}(z)=e^{i\phi_{j}(z)}\sqrt{n_{1D}+\rho_{j}(z)}\
.$ (50)

Note that $\phi_{j}$ and $\rho_{j}$ are the phase and density variables for
each well rather than their antisymmetric versions which are used extensively
in the main text.

Let us start by manipulating the kinetic energy term

$\displaystyle-\sum_{j=1,2}\int_{-\infty}^{\infty}$ $\displaystyle
dz\,\frac{\hbar^{2}}{2m}\hat{\psi}_{j}^{\dagger}(z)\frac{\partial^{2}\hat{\psi}_{j}(z)}{\partial
z^{2}}$ (51) $\displaystyle=\int_{-\infty}^{\infty}dz$
$\displaystyle\sum_{j=1,2}\frac{\hbar^{2}}{2m}\bigg{[}\left(\frac{\partial}{\partial
z}e^{-i\phi_{j}(z)}\sqrt{n_{1D}+\rho_{j}(z)}\right)$
$\displaystyle\times\left(\frac{\partial}{\partial
z}e^{+i\phi_{j}(z)}\sqrt{n_{1D}+\rho_{j}(z)}\right)\bigg{]}$
$\displaystyle=\int_{-\infty}^{\infty}dz$
$\displaystyle\sum_{j=1,2}\frac{\hbar^{2}}{2m}\bigg{(}-i\frac{\partial\phi_{j}}{\partial
z}\hat{\psi}_{j}^{\dagger}+\frac{e^{-i\phi_{j}}\frac{\partial\rho_{j}}{\partial
z}}{2\sqrt{n_{1D}+\rho_{j}}}\bigg{)}$
$\displaystyle\times\bigg{(}i\frac{\partial\phi_{j}}{\partial
z}\hat{\psi}_{j}+\frac{e^{i\phi_{j}}\frac{\partial\rho_{j}}{\partial
z}}{2\sqrt{n_{1D}+\rho_{j}}}\bigg{)}$
$\displaystyle=\int_{-\infty}^{\infty}dz$
$\displaystyle\sum_{j=1,2}\frac{\hbar^{2}}{2m}\bigg{\\{}\hat{\psi}_{j}^{\dagger}\hat{\psi}_{j}\left(\frac{\partial\phi_{j}}{\partial
z}\right)^{2}+\frac{(\frac{\partial\rho_{j}}{\partial
z})^{2}}{4(n_{1D}+\rho_{j})}$
$\displaystyle+i\frac{\frac{\partial\rho_{j}}{\partial
z}\frac{\partial\phi_{j}}{\partial
z}}{2\sqrt{n_{1D}+\rho_{j}}}[\hat{\psi}_{j}e^{-i\phi_{j}}-\hat{\psi}_{j}^{\dagger}e^{i\phi_{j}}]\bigg{\\}}$
$\displaystyle=\int_{-\infty}^{\infty}dz$
$\displaystyle\sum_{j=1,2}\frac{\hbar^{2}}{2m}\bigg{\\{}\hat{\psi}_{j}^{\dagger}\hat{\psi}_{j}\left(\frac{\partial\phi_{j}}{\partial
z}\right)^{2}+\frac{(\frac{\partial\rho_{j}}{\partial
z})^{2}}{4(n_{1D}+\rho_{j})}\bigg{\\}}$
$\displaystyle\approx\int_{-\infty}^{\infty}dz$ $\displaystyle\
\frac{\hbar^{2}}{2m}\Bigg{\\{}\frac{n_{1D}}{2}\left[\left(\frac{\partial\phi_{s}}{\partial
z}\right)^{2}+\left(\frac{\partial\phi_{a}}{\partial z}\right)^{2}\right]$
$\displaystyle+\frac{1}{2n_{1D}}\left[\left(\frac{\partial\rho_{s}}{\partial
z}\right)^{2}+\left(\frac{\partial\rho_{a}}{\partial
z}\right)^{2}\right]\Bigg{\\}}$ (52)

where

$\displaystyle\phi_{a}=\phi_{1}-\phi_{2},\quad\phi_{s}=\phi_{1}+\phi_{2}$ (53)
$\displaystyle\rho_{a}=\frac{\rho_{1}-\rho_{2}}{2},\quad\rho_{s}=\frac{\rho_{1}+\rho_{2}}{2},$
(54)

and we assume that $n_{1D}\gg\rho_{j}$. Next we consider the interactions

$\displaystyle\sum_{j=1,2}$
$\displaystyle\frac{g_{1D}}{2}\psi_{j}^{\dagger}\psi_{j}^{\dagger}\psi_{j}\psi_{j}=\sum_{j=1,2}\frac{g_{1D}}{2}[n_{1D}+\rho_{j}(z)]^{2}$
$\displaystyle=$
$\displaystyle\sum_{j=1,2}\left(\frac{g_{1D}n_{1D}^{2}}{2}+\frac{g_{1D}\rho_{j}^{2}}{2}+g_{1D}n_{1D}\rho_{j}\right)$
$\displaystyle=$ $\displaystyle
g_{1D}n_{1D}^{2}+g_{1D}(\rho_{s}^{2}+\rho_{a}^{2})+2g_{1D}n_{1D}\rho_{s}\ .$
(55)

Finally, we consider the tunneling term

$\displaystyle-\hbar$ $\displaystyle
J\Big{[}\psi_{1}^{\dagger}(z)\psi_{2}(z)+\psi_{2}^{\dagger}(z)\psi_{1}(z)\Big{]}$
(56) $\displaystyle=$ $\displaystyle-\hbar
J\Big{[}(e^{-i(\phi_{1}-\phi_{2})}+e^{-i(\phi_{2}-\phi_{1})})\sqrt{n_{1D}+\rho_{1}}\sqrt{n_{1D}+\rho_{2}}\Big{]}$
$\displaystyle=$ $\displaystyle-2\hbar
J\cos(\phi_{a})\sqrt{n_{1D}+\rho_{1}}\sqrt{n_{1D}+\rho_{2}}$ $\displaystyle=$
$\displaystyle-2\hbar
J\cos(\phi_{a})\sqrt{n_{1D}^{2}+2n_{1D}\rho_{s}+\rho_{s}^{2}-\rho_{a}^{2}}$
$\displaystyle\approx$ $\displaystyle-2\hbar
J\cos(\phi_{a})(n_{1D}+\rho_{s})\approx-2\hbar n_{1D}J\cos(\phi_{a})\ .$ (57)

At very low temperatures the symmetric and antisymmetric components decouple
and hence can be treated separately. The lower energy terms are the
antisymmetric ones and we obtain the following Hamiltonian

$\begin{split}H_{\mathrm{SG}+}=\int_{-\infty}^{\infty}dz\bigg{[}g_{1D}&\,\rho_{a}(z)^{2}+\frac{\hbar^{2}n_{1D}}{4m}\left(\frac{\partial\phi_{a}}{\partial
z}\right)^{2}\\\
&+\frac{\hbar^{2}}{4mn_{1D}}\left(\frac{\partial\rho_{a}}{\partial
z}\right)^{2}\bigg{]}\\\ -\int_{-\infty}^{\infty}dz&\ 2\hbar
Jn_{1D}\cos\left[\phi_{a}(z)\right]\ .\end{split}$ (58)

When the higher wavelength $\rho$ modes are suppressed this reduces to the
sine-Gordon model

$\begin{split}H_{\mathrm{SG}}=&\int_{-\infty}^{\infty}dz\bigg{\\{}g_{1D}\,\rho_{a}(z)^{2}+\frac{\hbar^{2}n_{1D}}{4m}\left(\frac{\partial\phi_{a}}{\partial
z}\right)^{2}\\\ &-2\hbar J\,n_{1D}\,\cos\left[\phi_{a}(z)\right]\bigg{\\}}\
.\end{split}$ (59)

Eq. (59) is the finally obtained SG Hamiltonian $H_{\mathrm{SG}}$ which is the
low energy description of two cigar shaped tunnel-coupled quasicondensates
[50, 74].

##  Appendix B Derivation of the Tomonaga-Luttinger (TL) Hamiltonian in
Fourier space

In this appendix we derive the Fourier space version of the Tomonaga-Luttinger
(TL) Hamiltonian. Starting from Eq. (25), and applying the discrete Fourier
decompositions given in Eq. (26) and Eq. (27), we have

$\displaystyle H_{\mathrm{TL+}}(ra)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}dz\,\frac{g_{1D}}{N_{L}+1}\left(\sum_{k=-N_{L}/2}^{N_{L}/2}\varrho_{k}e^{i\frac{2\pi
kr}{N_{L}+1}}\right)\times\left(\sum_{l=-N_{L}/2}^{N_{L}/2}\varrho_{l}e^{i\frac{2\pi
lr}{N_{L}+1}}\right)$ (60) $\displaystyle+$
$\displaystyle\int_{-\infty}^{\infty}dz\,\frac{\hbar^{2}n_{1D}}{4ma^{2}(N_{L}+1)}\frac{\partial}{\partial
r}\left(\sum_{k=-N_{L}/2}^{N_{L}/2}\varphi_{k}e^{i\frac{2\pi
kr}{N_{L}+1}}\right)\times\frac{\partial}{\partial
r}\left(\sum_{l=-N_{L}/2}^{N_{L}/2}\varphi_{l}e^{i\frac{2\pi
lr}{N_{L}+1}}\right)$ $\displaystyle+$
$\displaystyle\int_{-\infty}^{\infty}dz\,\frac{\hbar^{2}}{4mn_{1D}a^{2}(N_{L}+1)}\frac{\partial}{\partial
r}\left(\sum_{k=-N_{L}/2}^{N_{L}/2}\varrho_{k}e^{i\frac{2\pi
kr}{N_{L}+1}}\right)\times\frac{\partial}{\partial
r}\left(\sum_{l=-N_{L}/2}^{N_{L}/2}\varrho_{l}e^{i\frac{2\pi
lr}{N_{L}+1}}\right)$ $\displaystyle=$ $\displaystyle
a\sum_{r=-N_{L}/2}^{N_{L}/2}\sum_{k=-N_{L}/2}^{N_{L}/2}\sum_{l=-N_{L}/2}^{N_{L}/2}\left(\frac{g_{1D}\varrho_{k}\varrho_{l}e^{i\frac{2\pi(k+l)r}{N_{L}+1}}}{N_{L}+1}\right)$
$\displaystyle-$ $\displaystyle
a\sum_{r=-N_{L}/2}^{N_{L}/2}\sum_{k=-N_{L}/2}^{N_{L}/2}\sum_{l=-N_{L}/2}^{N_{L}/2}\frac{\hbar^{2}n_{1D}}{4ma^{2}(N_{L}+1)}\times\left(\frac{2\pi}{N_{L}+1}\right)^{2}kl\varphi_{k}\varphi_{l}e^{i\frac{2\pi(k+l)r}{N_{L}+1}}$
$\displaystyle-$ $\displaystyle
a\sum_{r=-N_{L}/2}^{N_{L}/2}\sum_{k=-N_{L}/2}^{N_{L}/2}\sum_{l=-N_{L}/2}^{N_{L}/2}\frac{\hbar^{2}}{4mn_{1D}a^{2}(N_{L}+1)}\times\left(\frac{2\pi}{N_{L}+1}\right)^{2}kl\varrho_{k}\varrho_{l}e^{i\frac{2\pi(k+l)r}{N_{L}+1}}$

where we have split the $z$ coordinate into $N_{L}+1$ grid points separated by
distance $a$ so that $z=r\,a$ where $r$ in an integer lying in the range
specified by Eq. (11). Using the fact that $N_{L}a=L$, and applying the
identity
$\sum_{r=-N_{L}/2}^{N_{L}/2}e^{i\frac{2\pi(k+l)r}{N_{L}+1}}=(N_{L}+1)\delta_{k,-l}$
we obtain

$\begin{split}H_{\mathrm{TL}+}\approx&a\sum_{k}\sum_{l}g_{1D}\varrho_{k}\varrho_{l}\delta_{k,-l}\\\
&-a\sum_{k}\sum_{l}\left(\frac{\hbar^{2}n_{1D}\pi^{2}}{mL^{2}}\right)kl\varphi_{k}\varphi_{l}\delta_{k,-l}\\\
&-a\sum_{k}\sum_{l}\left(\frac{\hbar^{2}\pi^{2}}{mn_{1D}L^{2}}\right)kl\varrho_{k}\varrho_{l}\delta_{k,-l}\end{split}$
(61)

where in the second term we have also replaced $a^{2}(N_{L}+1)^{2}$ by $L^{2}$
which holds when $N_{L}\gg 1$. The limits of the summation in Eq. (61) has
been omitted for the sake of brevity. We therefore find

$\begin{split}H_{\mathrm{TL}+}\approx\sum_{k}\bigg{[}ag_{1D}\varrho_{k}\varrho_{-k}+&\frac{a\hbar^{2}n_{1D}\pi^{2}k^{2}}{mL^{2}}\varphi_{k}\varphi_{-k}\\\
&+\frac{a\hbar^{2}\pi^{2}k^{2}}{mn_{1D}L^{2}}\varrho_{k}\varrho_{-k}\bigg{]}\\\
=\sum_{k}\bigg{[}ag_{1D}|\varrho_{k}|^{2}+&\frac{a\hbar^{2}n_{1D}\pi^{2}k^{2}}{mL^{2}}|\varphi_{k}|^{2}\\\
&+\frac{a\hbar^{2}\pi^{2}k^{2}}{mn_{1D}L^{2}}|\varrho_{k}|^{2}\bigg{]}\end{split}$
(62)

where we used the property of real fields that

$\varphi_{-k}=\varphi^{\star}_{k},\quad\mbox{and}\quad\varrho_{-k}=\varrho^{\star}_{k}\
.$ (63)

Hence the Hamiltonian takes the form given in Eq. (28) of the main text.

##  Appendix C Bench marking of the numerical method

The results given in this paper rely on numerically evolving the equations of
motion over time for various models [e.g. for the full SG+ model the equations
of motion are given in Eq. (22)], which we accomplish using the Julia package
DifferentialEquations.jl [114]. This implements a Runge-Kutta solver with a
user-defined time step. As a measure of the accuracy of our numerical method
we use the deviation of the Hamiltonian from its initial value. Since the
Hamiltonian should be a constant of motion this gives an indication of the
size of the numerical errors.

In Figures 12 and 13 we plot the relative error in the SG+ Hamiltonian given
in Eq. (18) for different time and spatial resolutions. More precisely, Figure
12 shows the effect of varying the time step $d\tilde{t}$, whereas Figure 13
shows the effect of varying the number of grid points $N_{L}$ which sets the
spatial step $d\tilde{z}$. In both cases we have evolved the system for a
total elapsed time of $\tilde{t}=1000$ which corresponds to the longest times
we use in this paper (for the calculation of the long-term distribution shown
in Figure 9), and also taken an ensemble average over 100 different
trajectories similar to those in Figure 5. Furthermore, we also performed a
moving time average of 30-time steps around $\tilde{t}=1000$ to average out
the effect of fast oscillations.

As expected, the relative error decreases as $d\tilde{t}$ and $d\tilde{z}$
decrease. For all the calculations in this paper we chose $d\tilde{t}=0.2$ and
$N_{L}=50$ because this keeps the relative error below 10 % and does not
significantly slow down the simulations.

---
Figure 12:  The relative error in the SG+ Hamiltonian is plotted here as a
function of the time step $d\tilde{t}$. The definition of the SG+ Hamiltonian
is given in Eq. 18 and should be a constant of the motion were it not for
numerical errors. The moving time average of relative error is evaluated after
propagating the equations of motion for a total elapsed time of
$\tilde{t}=1000$. All parameter values are the same as in Figure 5 including
$N_{L}=50$.
---
Figure 13:  The relative error in the SG+ Hamiltonian is plotted here as a
function of the number of lattice points $N_{L}$ on the numerical spatial
lattice. Like in Figure 12, the Hamiltonian is evaluated after evolving the
equations of motion for a total elapsed time of $\tilde{t}=1000$. The moving
time average of the relative error fluctuates (at around 10 %) but does
decrease as d$\tilde{z}$ decreases (or $N_{L}$ increases). All other parameter
values are the same as in Figure 5 with $d\tilde{t}=0.2$

##  Appendix D Caustic curve

In this appendix we use the exact solution for the motion of a pendulum to
calculate the caustic curve plotted as the solid black line in Figure 3. The
caustic is in fact the envelope of a whole family of trajectories. To begin,
we take the equations of motion for the SG model given in Eq. (22) and drop
the second order derivative term proportional to $\epsilon$ which couples the
different pendula. Next, we make the change of variables

$\tilde{t}=At,\quad\tilde{\rho}=Bp,\quad\tilde{\phi}=2y$ (64)

where

$A=\frac{1}{2}\frac{1}{\sqrt{\mathcal{J}\Gamma}}\quad,\quad
B=2\sqrt{\frac{\mathcal{J}}{\Gamma}}$ (65)

so the equations of motion simplify to

$\displaystyle\frac{dy}{dt}$ $\displaystyle=$ $\displaystyle p$ (66)
$\displaystyle\frac{dp}{dt}$ $\displaystyle=$ $\displaystyle-\frac{1}{2}\sin
2y\ .$ (67)

These equations are Hamilton’s equations obtained from a standard pendulum
hamiltonian of the form

$H(y,p)=\frac{p^{2}}{2}+\frac{1}{2}\sin^{2}y\ .$ (68)

The equations of motion given in Eqns. (66) and (67) have exact solutions in
terms of the Jacobi elliptic functions $\mathrm{sn}[u|m]$ and
$\mathrm{cn}[u|m]$ [122]. For the case relevant to us where the pendulum
starts at angle $y_{0}$, with zero initial angular momentum, they are

$\displaystyle y(t,y_{0})$ $\displaystyle=$ $\displaystyle\arcsin\\{\sin
y_{0}\ \mathrm{sn}[t+K(\sin y_{0})|\sin y_{0}]\\}$ (69) $\displaystyle
p(t,y_{0})$ $\displaystyle=$ $\displaystyle\sin(y_{0})\ \mathrm{cn}[t+K(\sin
y_{0})|\sin y_{0}]$ (70)

where $K(m)=\int_{0}^{\pi/2}d\theta/\sqrt{1-m^{2}\sin^{2}\theta}$ is the
complete elliptic integral of the first kind [122] (we caution the reader that
some computer packages such as Mathematica use the syntax $K(m^{2})$ for this
integral).

Caustics occur when trajectories are focused, in other words they are the
places where the trajectory does not change (to first order) when the initial
conditions are varied. Thus, caustics in the momentum variable $p$ occur when
$dp/dy_{0}=0$ since the initial condition here is specified by $y_{0}$. By
differentiating Eq. (70) an implicit expression for the position of the
caustics can be found [123]

$\displaystyle\mathrm{sn}(u|m)\mathrm{dn}(u|m)\left(\frac{\mathrm{E}(\mathrm{am}(-t|m)\,|m)}{\cos(y_{0})}+t\cos(y_{0})\right)$
$\displaystyle-\cos(y_{0})\mathrm{cn}(u|m)=0$ (71)

where $u=t+K(\sin y_{0})$, $m=\sin y_{0}$, $\mathrm{E}(u|m)$ is an elliptic
integral of the second kind, $\mathrm{dn}(u|m)$ is another Jacobi elliptic
function, and $\mathrm{am}(u|m)=\arcsin[\sin(\phi)/m]$ is the Jacobi amplitude
[122]. Finding the roots $y_{0}$ of Eq. (71) numerically at each value of the
time gives pairs of values $(y_{0},t)$ that can then be put back into Eq. (70)
to yield the black curve for the caustic shown in Figure 3. The match to the
numerics is very good.

##  Appendix E Derivation of ergodic (“circus tent”) probability distribution
at long times

In this appendix we outline the derivation of an analytic approximation to the
probability distribution for the number difference at long times, as shown in
Figure 9. This derivation is based upon a calculation given in Ref. 124 and
assumes that the average behaviour of a continuous chain of coupled pendula
(the mechanical system that underlies the sine-Gordon model) can be described
by a suitably ‘ergodized’ single pendulum.

To keep the calculation general we use the pendulum Hamiltonian in standard
form as given in Eq. (68). With this hamiltonian we define a microcanonical
probability density in phase space:

$d_{m}(y,p;y_{0})=\frac{\delta[H(y,p)-H(y_{0},p)]}{\int\int dy\ dp\
\delta[H(y,p)-H(y_{0},p)]}$ (72)

where $y_{0}$ is the initial angle of the pendulum which fixes its total
energy to be $E=(1/2)\sin^{2}y_{0}$ if the the initial angular momentum is
zero (this is the appropriate initial condition for the tunneling quench
considered in this paper where the initial number difference is taken to be
zero), and the denominator ensures that $d_{m}$ is normalized to unity. A
microcanonical distribution has equal probability to be anywhere on its energy
shell (in this case a closed curve in $y,p$ phase space) and thus by adopting
Eq. (72) we are making an ergodic hypothesis. This does not hold for a single
pendulum starting at position $y_{0}$ since it will spend the most time at its
turning points $y=\pm y_{0}$, but when averaged over $y_{0}$ and $y$ (see
below) it gives a very good approximation at long times, as can be seen in
Figure 9.

The normalization integral can be evaluated exactly by re-expressing the delta
function using the relation
$\delta[g(x)]=\sum_{i}\delta(x-x_{i})/|g^{\prime}(x_{i})|$, where $x_{i}$ are
the roots of $g(x)$. In the present case this gives

$\begin{split}\delta[(p^{2}+\sin^{2}y-\sin^{2}y_{0})/2]=&\frac{\delta[p-p_{1}]}{|p_{1}|}+\frac{\delta[p-p_{2}]}{|p_{2}|}\\\
=&2\frac{\delta[p-p_{1}]}{|p_{1}|}\end{split}$ (73)

where $|p_{1}|=|p_{2}|=\sqrt{\sin^{2}y_{0}-\sin^{2}y}$. In obtaining this
expression we have used the fact that for values of $y$ within the range
accessed by the pendulum, there are two values of $p$ where the integral
crosses the energy shell. The integral over $p$ is now trivial due to the
delta function and the integral over $y$ can be performed by putting $\sin
y=\sin y_{0}\sin\zeta$ so that

$\begin{split}2\int_{-y_{0}}^{y_{0}}\frac{dy}{|p(y,y_{0})|}=2\int_{-y_{0}}^{y_{0}}\frac{dy}{\sqrt{\sin^{2}y_{0}-\sin^{2}y}}=2\int_{-\pi/2}^{\pi/2}\frac{d\zeta}{\sqrt{1-\sin^{2}y_{0}\sin^{2}\zeta}}=4\int_{0}^{\pi/2}\frac{d\zeta}{\sqrt{1-\sin^{2}y_{0}\sin^{2}\zeta}}=4K(\sin
y_{0})\ .\end{split}$ (74)

Therefore, the normalized microcanonical probability density can be written as

$\begin{split}d_{m}(y,p;y_{0})=&\frac{1}{4K(\sin
y_{0})}\delta[(p^{2}+\sin^{2}y-\sin^{2}y_{0})/2]\\\ =&\frac{1}{2K(\sin
y_{0})}\delta(p^{2}+\sin^{2}y-\sin^{2}y_{0})\end{split}$ (75)

where we have used the property of delta functions that $\delta(\alpha
x)=(1/\alpha)\delta(x)$.

The initial condition for our dynamics is such that the number difference is
well defined but the phase difference is completely undefined. We must
therefore average the microcanonical probability density over all $y_{0}$.
This gives the phase space probability density relevant to $J$-quenches as
being

$W(y,p)=\frac{1}{\pi}\int_{-\pi/2}^{\pi/2}dy_{0}\ d_{m}(y,p;y_{0})$ (76)

where we employ the notation $W$ to indicate that this is a classical version
of the Wigner function. The properties of the delta function can once more be
used to write

$\delta(p^{2}+\sin^{2}y-\sin^{2}y_{0})=\sum_{i}\frac{\delta(y_{0}-y_{0i})\theta(\cos
y-|p|)}{2\sqrt{p^{2}+\sin^{2}y}\sqrt{\cos^{2}y-p^{2}}}$ (77)

where $\theta(x)$ is the Heaviside step function. The integral over $y_{0}$
can now be evaluated exactly to give

$W(y,p)=\frac{2}{4\pi}\frac{\theta(\cos
y-|p|)}{K(\sqrt{p^{2}+\sin^{2}y})\sqrt{p^{2}+\sin^{2}y}\sqrt{\cos^{2}y-p^{2}}}.$
(78)

The final step is to integrate out the $y$ coordinate to obtain the
probability distribution $\mathcal{P}_{\mathrm{CT}}(p)$ for $p$ alone

$\mathcal{P}_{\mathrm{CT}}(p)=\int_{-\pi/2}^{\pi/2}dy\ W(y,p)\ ,$ (79)

where “CT” stands for circus tent. Although this integral cannot be done
analytically, it can be put in a form which is convenient to evaluate
numerically. Denoting $m=\sin y_{0}=\sqrt{p^{2}+\sin^{2}y}$, one finds that

$\mathcal{P}_{\mathrm{CT}}(\tilde{\rho})=\frac{1}{2\pi
B}\int_{\tilde{\rho}^{2}/B^{2}}^{1}\frac{dm}{K(m)\sqrt{m(1-m)(m-\tilde{\rho}^{2}/B^{2})(1+\tilde{\rho}^{2}/B^{2}-m)}}$
(80)

where we have also converted back from angular momentum $p$ to number
difference $\tilde{\rho}$ using Eq. (64). This equation is given in the main
text as Eq. (44) and is plotted in Figure 9 where it is compared against the
long-time spatially and temporally averaged numerical data for the various
nonlinear models considered in this paper. As can be seen in Figure 9,
$\mathcal{P}_{\mathrm{CT}}$ is characterized by a diverging (yet normalizable)
peak at the center and then relatively flat wings until it drops sharply to
zero at the edges. In Ref. 124 it is shown that
$\mathcal{P}_{\mathrm{CT}}(\tilde{\rho})$ diverges logarithmically at the
origin $\tilde{\rho}=0$ and also tends suddenly to zero with logarithmic
singularities at $\tilde{\rho}=\pm B$. Both these non-thermal features can be
attributed to the presence of caustics.

##  Appendix F Pendulum at thermal equilibrium

In Figure 9 the long time probability distribution for the number difference
is compared against the ergodic prediction derived in Appendix E, and also
against the thermal equilibrium prediction. In this Appendix we explain how to
calculate the latter case. In order to make the calculation tractable we make
the assumption that the SG+ model can be approximated by a thermal ensemble of
independent pendula. We also adopt the same notation as Appendix E and hence
work with a pendulum Hamiltonian in the standard form
$H=(1/2)(p^{2}+\sin^{2}y)$. This is related to the two mode Hamiltonian
$H_{2M}=\Gamma\tilde{\rho}^{2}-2\mathcal{J}\cos\phi$ by
$H=H_{2M}/8\mathcal{J}+1/4$.

We proceed in two steps: we first calculate the probability distribution
$\mathcal{P}_{E}(p)$ for the momentum variable $p$ (that here plays the role
of the number difference) for a fixed energy $E$. Secondly, we assume our
system is at thermal equilibrium with a bath at temperature $T$ such that the
relative probability of any energy is given by the Boltzmann factor
$\exp[-E/T]$. Thus the thermal probability distribution is

$\mathcal{P}_{T}(p)=\frac{1}{Z}\int_{0}^{\infty}\mathcal{P}_{E}(p)\
e^{-E/T}D(E)\ dE$ (81)

where $Z$ is a normalizing factor (found numerically) and $D(E)$ is the
density of states.

The probability distribution $\mathcal{P}_{E}(p)$ at fixed $E$ is proportional
to $1/\dot{p}$ as this determines how long the pendulum spends at each value
of $p$. According to Hamilton’s equation $\dot{p}=-\partial H/\partial
x=-(1/2)\sin 2y$, and using the fact that $\sin y=\sqrt{2E-p^{2}}$, we find
that this probability distribution for a fixed value of $E$ is

$\mathcal{P}_{E}(p)=\frac{\mathcal{N}}{(1/2)\sin(2\arcsin{\sqrt{2E-p^{2}}})},$
(82)

where $\mathcal{N}$ is a normalization factor given by the period of the
motion. Two cases must be distinguished: for $E<1/2$ the energy is less than
the separatrix and the pendulum undergoes vibrational motion (also known as
librational motion in some literature). Conversely, when $E>1/2$ the energy is
above the separatrix and the pendulum undergoes rotational motion.

For motion below the separatrix we have $|p|<p_{\mathrm{max}}=\sqrt{2E}$. We
must therefore supplement the expression for $\mathcal{P}_{E}(p)$ with the
condition that it is zero if $|p|>p_{\mathrm{max}}$ and this ensures that
$\mathcal{P}_{E}(p)$ is real. $\mathcal{N}$ is given in this case by

$\mathcal{N}=\frac{1}{2\ K(\sqrt{2E})}$ (83)

where, as in Appendix E, $K$ is the complete elliptic integral of the first
kind.

For motion above the separatrix we have $\sqrt{2E-1}<|p|<\sqrt{2E}$ and
$\mathcal{P}_{E}(p)$ is zero outside this range. $\mathcal{N}$ is now given by

$\mathcal{N}=\frac{\sqrt{2E}}{4\ K(1/\sqrt{2E})}.$ (84)

To obtain the total thermal probability distribution $\mathcal{P}_{T}(p)$
given in Eq. (81) we need the density of states $D(E)\equiv dn/dE$, where $n$
is the number of states below energy $E$. According to the Bohr-Sommerfeld
rule $n=S(E)/(2\pi\hbar)$, where the action $S(E)=\oint p\,dy$ is the area in
phase space enclosed by the energy contour $E$. However, assuming that our
Hamiltonian $H$ is in units $\hbar\omega$ then the $2\pi\hbar$ factor is
absorbed into the definitions of $p$ and $y$ and we have $D(E)=(d/dE)\oint
p\,dy$. Below the separatrix we have

$\oint p(y)dy=4\int_{0}^{\arcsin\sqrt{2E}}\sqrt{2E-\sin^{2}y}\ dy$ (85)

and putting $2E=\sin^{2}y_{0}$ we find

$\displaystyle D_{<}(E)=$ $\displaystyle\frac{d}{dE}$ $\displaystyle\oint
p(y)dy$ (86) $\displaystyle=$ $\displaystyle
4\int_{0}^{\arcsin\sqrt{2E}}\frac{dy}{\sqrt{\sin^{2}y_{0}-\sin^{2}y}}$
$\displaystyle=$ $\displaystyle 4K(\sqrt{2E})$

where the integral is performed in a similar fashion to the one in Eq. (74)
and the subscript “$<$” indicates that this is the expression valid below the
separatrix. Above the separatrix we find that the area enclosed in phase space
between two oppositely rotating states of the same energy is

$\oint p(y)dy=2\int_{-\pi/2}^{\pi/2}\sqrt{2E-\sin^{2}y}\ dy$ (87)

and thus

$\displaystyle D_{>}(E)=$ $\displaystyle\frac{d}{dE}$ $\displaystyle\oint
p(y)dy$ (88) $\displaystyle=$ $\displaystyle
2\int_{-\pi/2}^{\pi/2}\frac{dy}{\sqrt{2E-\sin^{2}y}}$ $\displaystyle=$
$\displaystyle\frac{4}{\sqrt{2E}}K\left(\frac{1}{\sqrt{2E}}\right)\ .$

Due to the fact that above the separatrix $2E>\sin^{2}y$ we no longer need to
make the substitutions $2E=\sin^{2}y_{0}$ and $\sin y=\sin y_{0}\sin\zeta$,
and the integral is straightforward. The subscript “$>$” indicates that this
expression holds above the separatrix.

We now have all the necessary ingredients to perform the integral for
$\mathcal{P}_{T}(p)$ which we do numerically. The two contributions, one from
below the separatrix and one from above, are added together to get the total.
Interestingly, both density of states factors, Eqns. (86) and (88), diverge at
the separatrix such that the two contributions individually display singular
features but remarkably these cancel out when the two parts are added and
result in the smooth gaussian curve plotted in Figure 9.

In order to compare the thermal distribution against the quenched (followed by
integrable SG evolution) distribution derived in Appendix E we need to choose
a temperature $T$ for the thermal distribution $\mathcal{P}_{T}$. We do this
by matching the expectation value of the energy $\langle E\rangle$ for both
distributions. In the quenched case the initial state corresponds to an
ensemble of pendula with different starting angles $y_{0}$ and zero kinetic
energy. Each starting angle in the range $-\pi/2<y_{0}\leq\pi/2$ is equally
probable in our $J$-quench. Therefore

$\langle
E\rangle_{\mathrm{quench}}=\frac{1}{\pi}\int_{-\pi/2}^{\pi/2}\frac{1}{2}\sin^{2}y_{0}\
dy_{0}=\frac{1}{4}\ .$ (89)

To calculate $\langle E\rangle$ in the thermal case we compute

$\langle E\rangle_{T}=\frac{1}{\zeta}\int_{0}^{\infty}E\,e^{-E/T}\,D(E)\,dE$
(90)

numerically for a large number of different values of $T$, performing the
integrals below and above the separatrix separately and adding the results.
Here $\zeta=\int_{0}^{\infty}e^{-E/T}\,D(E)\,dE$ gives the normalization
factor. We then fit a curve to the results and find the value of $T$ that best
matches the result given in Eq. (89). We find that $T=0.184$ gives the best
match. Putting back the units this result is

$\frac{k_{B}T}{8\mathcal{J}\hbar c/\xi_{h}}=\frac{k_{B}T}{16J\hbar
K/\pi}=0.184$ (91)

where $c$ is the speed of sound and $K$ is the Luttinger parameter and $J$ is
the tunnel coupling rate between the two wells. In this paper we take $K=25$
and $J=30$ Hz (see Table 1) giving a temperature in SI units of 5.4 nK.

## References

  * Nye [1999] J. F. Nye,  _Natural Focusing and Fine Structure of Light: Caustics and Wave Dislocations_ (Institute of Physics Publishing: Bristol and Philadelphia, 1999).
  * Lord Kelvin [1905] Lord Kelvin, Deep water ship-waves, Phil. Mag. 9, 733 (1905).
  * Titov  _et al._ [2005] V. V. Titov, A. B. Rabinovich, H. O. Mofjeld, and F. I. Thomson R. E. González, The global reach of the 26 December 2004 Sumatra tsunami, Science 309, 2045 (2005).
  * Berry [2007] M. V. Berry, Focused tsunami waves, Proc. R. Soc. London A 463, 3055 (2007).
  * Degueldre  _et al._ [2016] H. Degueldre, J. J. Metzger, T. Geisel, and R. Fleischmann, Random focusing of tsunami waves, Nat. Phys. 12, 259 (2016).
  * Berry [2018] M. V. Berry, Minimal analytical model for undular tidal bore profile; quantum and hawking effect analogies, New J. Phys. 20, 053066 (2018).
  * Heller  _et al._ [2021] E. J. Heller, R. Fleischmann, and T. Kramer, Branched flow, Physics Today 74, 44 (2021).
  * Peebles [1980] P. J. E. Peebles,  _The Large-Scale Structure of the Universe_ (Princeton University Press, Princeton, NJ, 1980).
  * Arnold  _et al._ [1982] V. I. Arnold, S. F. Shandarin, and Y. B. Zeldovich, The large scale structure of the universe I. General properties. One- and two-dimensional models, Geophys. Astrophys. Fluid Dyn. 20, 111 (1982).
  * Gurbatov  _et al._ [2012] S. N. Gurbatov, A. I. Saichev, and S. F. Shandarin, Large-scale structure of the universe. The Zeldovich approximation and the adhesion model, Phys.-Usp. 55, 223 (2012).
  * Feldbrugge  _et al._ [2018] J. Feldbrugge, R. van de Weygaert, J. Hidding, and J. Feldbrugge, Caustic skeleton & cosmic web, J. Cosmol. Astropart. Phys. 5, 27.
  * Berry [1977] M. V. Berry, Regular and irregular semiclassical wavefunctions, J. Phys. A: Math. Gen. 10, 2083 (1977).
  * Berry [1981] M. Berry, Singularities in waves and rays, in  _Physics of Defects (Les Houches Session XXXV)_, edited by R. Balian, M. Kléman, and J.-P. Poirier (North-Holland Publishing, Amsterdam, 1981).
  * Thom [1975] R. Thom,  _Structural Stability and Morphogenesis_ (Benjamin, Reading, MA, 1975).
  * Arnol’d [1975] V. I. Arnol’d, Critical points of smooth functions and their normal forms, Russ. Math. Survs. 30, 1 (1975).
  * Zeeman [1977] E. C. Zeeman,  _Catastrophe Theory: Selected Papers 1972-1977_ (Addison-Wesley, Reading, MA, 1977).
  * Mumford  _et al._ [2017] J. Mumford, W. Kirkby, and D. H. J. O’Dell, Catastrophes in non-equilibrium many-particle wave functions: universality and critical scaling, J. Phys. B: At. Mol. Opt. Phys. 50, 044005 (2017).
  * Nesvizhevsky  _et al._ [2002] V. V. Nesvizhevsky, H. G. Börner, A. K. Petukhov, H. Abele, S. Baeßler, F. J. Rueß, T. Stöferle, A. Westphal, A. M. Gagarski, G. A. Petrov, and A. V. Strelkov, Quantum states of neutrons in the Earth’s gravitational field, Nature 415, 297 (2002).
  * Jenke  _et al._ [2011] T. Jenke, P. Geltenbort, H. Lemmel, and H. Abele, Realization of a gravity-resonance-spectroscopy technique, Nature Physics 7, 468 (2011).
  * Petersen  _et al._ [2013] T. C. Petersen, M. Weyland, D. M. Paganin, T. P. Simula, S. A. Eastwood, and M. J. Morgan, Electron vortex production and control using aberration induced diffraction catastrophes, Phys. Rev. Lett. 110, 033901 (2013).
  * Rooijakkers  _et al._ [2003] W. Rooijakkers, S. Wu, P. Striehl, M. Vengalattore, and M. Prentiss, Observation of caustics in the trajectories of cold atoms in a linear magnetic potential, Phys. Rev. A 68, 063412 (2003).
  * Huckans  _et al._ [2009] J. H. Huckans, I. B. Spielman, B. L. Tolra, W. D. Phillips, and J. V. Porto, Quantum and classical dynamics of a Bose-Einstein condensate in a large-period optical lattice, Phys. Rev. A 80, 043609 (2009).
  * Rosenblum  _et al._ [2014] S. Rosenblum, O. Bechler, I. Shomroni, R. Kaner, T. Arusi-Parpar, O. Raz, and B. Dayan, Demonstration of fold and cusp catastrophes in an atomic cloud reflected from an optical barrier in the presence of gravity, Phys. Rev. Lett. 112, 120403 (2014).
  * Mossman  _et al._ [2021] M. E. Mossman, T. M. Bersano, M. McNeil Forbes, and P. Engels, Gravitational caustics in an atom laser, Nature Communications 12, 7226 (2021).
  * Simula  _et al._ [2013] T. P. Simula, T. C. Petersen, and D. M. Paganin, Diffraction catastrophes threaded by quantized vortex skeletons caused by atom-optical aberrations induced in trapped Bose-Einstein condensates, Phys. Rev. A 88, 043626 (2013).
  * Mumford  _et al._ [2019] J. Mumford, E. Turner, D. W. L. Sprung, and D. H. J. O’Dell, Quantum spin dynamics in Fock space following quenches: Caustics and vortices, Phys. Rev. Lett. 122, 170402 (2019).
  * Kirkby  _et al._ [2019] W. Kirkby, J. Mumford, and D. H. J. O’Dell, Quantum caustics and the hierarchy of light cones in quenched spin chains, Phys. Rev. Research 1, 033135 (2019).
  * Berry and Dennis [2004] M. V. Berry and M. R. Dennis, Quantum cores of optical phase singularities, J. Opt. A 6, S178 (2004).
  * Berry [2008] M. V. Berry, Three quantum obsessions, Nonlinearity 21, T19 (2008), publisher: Institute of Physics Publishing.
  * Leonhardt [2002] U. Leonhardt, A laboratory analogue of the event horizon using slow light in an atomic medium, Nature (London) 415, 406 (2002).
  * O’Dell [2012] D. H. J. O’Dell, Quantum catastrophes and ergodicity in the dynamics of bosonic Josephson junctions, Phys. Rev. Lett. 109, 150406 (2012).
  * Goldberg  _et al._ [2019] A. Z. Goldberg, A. Al-Qasimi, J. Mumford, and D. H. J. O’Dell, Emergence of singularities from decoherence: Quantum catastrophes, Phys. Rev. A 100, 063628 (2019).
  * Plestid  _et al._ [2018] R. Plestid, P. Mahon, and D. H. J. O’Dell, Violent relaxation in quantum fluids with long-range interactions, Phys. Rev. E 98, 012112 (2018).
  * Kirkby  _et al._ [2022] W. Kirkby, Y. Yee, K. Shi, and D. H. J. O’Dell, Caustics in quantum many-body dynamics, Phys. Rev. Res. 4, 013105 (2022).
  * Berry and Wright [1980] M. V. Berry and F. J. Wright, Phase-space projection identities for diffraction catastrophes, J. Phys. A: Math. Gen. 13, 149 (1980).
  * Coleman [1975] S. Coleman, Quantum sine-Gordon equation as the massive Thirring model, Phys. Rev. D 11, 2088 (1975).
  * Rajaraman [1975] R. Rajaraman, Some non-perturbative semi-classical methods in quantum field theory (a pedagogical review), Phys. Rep. 21, 227 (1975).
  * Fogel  _et al._ [1977] M. B. Fogel, S. E. Trullinger, A. R. Bishop, and J. A. Krumhansl, Dynamics of sine-Gordon solitons in the presence of perturbations, Phys. Rev. B 15, 1578 (1977).
  * Malomed [2014] B. A. Malomed, The sine-Gordon model. General background, physical motivations, inverse scattering, and solitons, in  _The sine-Gordon model and its applications: from pendula and Josephson junctions to gravity and high-energy physics_, edited by J. Cuevas-Maraver, P. G. Kevrekidis, and F. Williams (Springer International Publishing, Heidelberg, 2014) pp. 1–30.
  * Frenkel and Kontorova [1939] J. Frenkel and T. Kontorova, On the theory of plastic deformation and twinning, Izvestiya Akademii Nauk SSSR, Seriya Fizicheskaya 1, 137 (1939).
  * Enz [1964] U. Enz, The motion of Bloch magnetic walls in magnetic crystals, Helv. Phys. Acta 37, 245 (1964).
  * Gallemí  _et al._ [2019] A. Gallemí, L. P. Pitaevskii, S. Stringari, and A. Recati, Decay of the relative phase domain wall into confined vortex pairs: The case of a coherently coupled bosonic mixture, Phys. Rev. A 100, 023607 (2019).
  * Oshikawa and Affleck [1997] M. Oshikawa and I. Affleck, Field-induced gap in $s=1/2$ antiferromagnetic chains, Phys. Rev. Lett. 79, 2883 (1997).
  * Affleck and Oshikawa [1999] I. Affleck and M. Oshikawa, Field-induced gap in Cu benzoate and other $s=1/2$ antiferromagnetic chains, Phys. Rev. B 60, 1038 (1999).
  * Cortés Cubero and Schuricht [2017] A. Cortés Cubero and D. Schuricht, Quantum quench in the attractive regime of the sine-Gordon model, J. Stat. Mech. , 103106 (2017).
  * Haldane [1981] F. D. M. Haldane, Effective harmonic-fluid approach to low-energy properties of one-dimensional quantum fluids, Phys. Rev. Lett. 47, 1840 (1981).
  * Giamarchi [2004] T. Giamarchi,  _Quantum Physics in One Dimension_ (Oxford University Press, New York, 2004).
  * Hadzibabic  _et al._ [2006] Z. Hadzibabic, P. Kruger, M. Cheneau, B. Battelier, and J. Dalibard, Berezinskii–Kosterlitz–Thouless crossover in a trapped atomic gas, Nature 441, 1118 (2006).
  * Chelpanova  _et al._ [2022] O. Chelpanova, S. P. Kelly, G. Morigi, F. Schmidt-Kaler, and J. Marino, Injection and nucleation of topological defects in the quench dynamics of the Frenkel-Kontorova model,  (2022), arXiv:2210.14904v2 [cond.mat] .
  * Bouchoule [2005] I. Bouchoule, Modulational instabilities in Josephson oscillations of elongated coupled condensates, Eur. Phys. J. D 35, 147 (2005).
  * Gritsev  _et al._ [2007] V. Gritsev, A. Polkovnikov, and E. Demler, Linear response theory for a pair of coupled one-dimensional condensates of interacting atoms, Phys. Rev. B 75, 174511 (2007).
  * Iucci and Cazalilla [2010] A. Iucci and M. A. Cazalilla, Quantum quench dynamics of the sine-Gordon model in some solvable limits, New J. Phys. 12, 055019 (2010).
  * Dalla Torre  _et al._ [2013] E. G. Dalla Torre, E. Demler, and A. Polkovnikov, Universal rephasing dynamics after a quantum quench via sudden coupling of two initially independent condensates, Phys. Rev. Lett. 110, 090404 (2013).
  * Foini and Giamarchi [2015] L. Foini and T. Giamarchi, Nonequilibrium dynamics of coupled luttinger liquids, Phys. Rev. A 91, 023627 (2015).
  * Schmiedmayer [2018] J. Schmiedmayer, One-dimensional atomic superfluids as a model system for quantum thermodynamics, in  _Thermodynamics in the Quantum Regime_, Vol. 195 (Springer Nature, 2018) p. 823.