Mechanisms for the emergence of Gaussian correlations

M. Gluza1,2*, T. Schweigler3,4, M. Tajik3, J. Sabino3,5,6, F. Cataldini3, F.
Møller3, S.-C. Ji3, B. Rauer3,7, J. Schmiedmayer3, J. Eisert1,8, S.
Sotiriadis1,9

1 Dahlem Center for Complex Quantum Systems, Freie Universität Berlin, 14195
Berlin, Germany
2 School of Physical and Mathematical Sciences, Nanyang Technological
University, 21 Nanyang Link, 637371 Singapore, Republic of Singapore
3 Vienna Center for Quantum Science and Technology, Atominstitut, TU Wien,
1020 Vienna, Austria
4 JILA, University of Colorado, Boulder, Colorado 80309-0440, USA
5 Instituto de Telecomunicações, Physics of Information and Quantum
Technologies Group, Av. Rovisco Pais 1, 1049-001, Lisbon, Portugal
6 Instituto Superior Técnico, Universidade de Lisboa, Av. Rovisco Pais 1,
1049-001, Lisbon, Portugal
7 Laboratoire Kastler Brossel, Ecole Normale Supérieure, 24 rue Lhomond,
F-75231 Paris, France
8 Helmholtz Center Berlin, 14109 Berlin, Germany
9 Department of Physics, University of Ljubljana, 1000 Ljubljana, Slovenia
*<EMAIL_ADDRESS>
## Abstract

We comprehensively investigate two distinct mechanisms leading to memory loss
of non-Gaussian correlations after switching off the interactions in an
isolated quantum system undergoing out-of-equilibrium dynamics. The first
mechanism is based on spatial scrambling and results in the emergence of
locally Gaussian steady states in large systems evolving over long times. The
second mechanism, characterized as ‘canonical transmutation’, is based on the
mixing of a pair of canonically conjugate fields, one of which initially
exhibits non-Gaussian fluctuations while the other is Gaussian and dominates
the dynamics, resulting in the emergence of relative Gaussianity even at
finite system sizes and times. We evaluate signatures of the occurrence of the
two candidate mechanisms in a recent experiment that has observed
Gaussification in an atom-chip controlled ultracold gas and elucidate evidence
that it is canonical transmutation rather than spatial scrambling that is
responsible for Gaussification in the experiment. Both mechanisms are shown to
share the common feature that the Gaussian correlations revealed dynamically
by the quench are already present though practically inaccessible at the
initial time. On the way, we present novel observations based on the
experimental data, demonstrating clustering of equilibrium correlations,
analyzing the dynamics of full counting statistics, and utilizing tomographic
reconstructions of quantum field states. Our work aims at providing an
accessible presentation of the potential of atom-chip experiments to explore
fundamental aspects of quantum field theories in quantum simulations.

###### Contents

  1. 1 Introduction
    1. 1.1 Two mechanisms for the emergence of Gaussian correlations
    2. 1.2 Experimental observation of the decay and recurrence of non-Gaussian correlations
  2. 2 Theoretical discussion of two readings of Gaussification
    1. 2.1 Gaussification by spatial scrambling
    2. 2.2 Gaussification by canonical transmutation
    3. 2.3 Similarities and differences of the two mechanisms
  3. 3 Characteristics of the initial correlations in the experiment
    1. 3.1 Clustering and scaling in initial states
      1. 3.1.1 Quantum field simulation of the relativistic field theory models
      2. 3.1.2 Velocity field correlations in the initial state
    2. 3.2 Canonical decoupling into Gaussian and non-Gaussian sectors
    3. 3.3 Dominance of density over phase initial correlations in the experimental states
      1. 3.3.1 Intuitive understanding of the experimental findings
      2. 3.3.2 Estimation of density fluctuations via phonon tomography
  4. 4 Characteristics of the dynamics in the experiment
    1. 4.1 Signatures of light-cone dynamics
    2. 4.2 Recurrences of the initial state
    3. 4.3 Deviations from the non-interacting phonon dynamics in the experiment
  5. 5 Characteristics of Gaussification in the experiment
    1. 5.1 Gaussification in the dynamics of full counting statistics
    2. 5.2 Temporal properties of the decaying connected correlations
    3. 5.3 Interpretation of the experimental data
  6. 6 Conclusions
  7. A Tomonaga-Luttinger liquid description of the experimental system
  8. B Tomonaga-Luttinger liquid dynamics of correlations in a homogeneous infinite system
  9. C Deviations from the Tomonaga-Luttinger liquid dynamics
    1. C.1 Effects of a non-uniform trap potential
    2. C.2 Effects of a non-linear phonon dispersion
  10. D Further tomographic reconstruction plots

##  1 Introduction

By appropriately choosing the effective degrees of freedom, it is frequently
possible to capture complex collective behavior of an interacting quantum
many-body system using a simple Gaussian model. There is an abundance of exact
analytical tools for computing physical properties using such Gaussian
effective descriptions. Most crucially, due to Wick’s theorem the second
moments are decisive for higher order correlation functions which do not yield
further information. Thanks to such concise features in sync with a high
predictive power, it is fair to say that Gaussian models are the bread and
butter of many physicists, independent of their focus, be it experimental or
theoretical.

The study of far from equilibrium quantum dynamics, however, hints at a
fundamental question concerning the applicability of such Gaussian models.
Whenever many-body interactions are sufficiently strong so that they can no
longer be meaningfully neglected, their properties become imprinted onto the
correlations in the system. But what will happen if these interactions are
suddenly removed? After such a quench protocol – a ubiquitous one actually
considered in a quite large body of literature [1, 2, 3, 4] – the ensuing time
evolution will be Gaussian and governed by a non-interacting quadratic
Hamiltonian. Will the state over time lose memory of those initial
interactions? And if so, in what precise way?

This is quite a complex physics question in that immediately after the
interaction quench, there is an apparent discrepancy between the state (which
is non-Gaussian) and the (Gaussian) Hamiltonian closed system dynamics. For
quantum states with Gaussian correlations, all correlations can be determined
from first and second moments of canonical coordinates using Wick’s theorem. A
natural expectation may be that for long times, the discrepancy will become
less stark, eventually the state becoming approximately Gaussian, in
accordance with the Hamiltonian governing the dynamics. A complication arises,
however, because quadratic Hamiltonians which give rise to Gaussian dynamics
feature a large number of conserved charges, which may be relevant in a given
quench scenario leading to intricate memory effects and hence complicating the
task of understanding the dynamics following an interaction quench.

More specifically, conserved charges may enable a physical realization of the
persistence of the discrepancy between the initial non-Gaussian state and the
quadratic Hamiltonian. Moreover, even independently of the presence of
conserved charges, the nature of the dynamics may influence the timescale for
the decay of the discrepancy between the non-Gaussian initial state and the
Gaussian quench dynamics. Thanks to the Gaussian character of the quench
dynamics, such memory loss of the signatures of interactions present prior to
the quench can be characterized theoretically for certain scenarios [5, 6, 7].
Thus, it is natural to ask whether the existing theoretical results can be
used to explain the decay of non-Gaussianity of a state whenever it appears
experimentally. As a wider perspective, the question that is on the desk is a
particular reading of the question of the emergence of a _generalized Gibbs
ensemble_ [8, 9, 10, 5, 11, 12]. Understanding the mechanism and conditions
for the emergence of Gaussian correlations is a good opportunity to understand
more generally the nature of equilibration in isolated quantum systems, since
Gaussian dynamics is an exceptionally convenient case for an in-depth theory-
experiment comparison. This is due to two main reasons. On the theoretical
side, exact analytical methods are available for the study of not only the
final steady state, but also for the full dynamics of correlations of any
order and for arbitrary initial states [13, 5, 6, 14, 15, 16, 17, 18, 7, 19,
20, 21]. From an experimental perspective, the possibility to measure and
characterize the factorization properties of higher-order correlations first
achieved in Ref. [22, 23] offers a practically complete description of quantum
states and their dynamics.

The main aim of this work is to elaborate on both theoretical and experimental
viewpoints on the question how quantum states become Gaussian over time in the
context of the experimental findings recently presented in Ref. [24]. In that
experimental work, non-Gaussian initial states have been prepared through the
coupling between two adjacent one-dimensional ultra-cold gases, and after a
fast decoupling that corresponds effectively to an interaction switch-off
quench, a decay and subsequently a revival of non-Gaussianity has been
observed. In the remainder of this introductory section, we will present the
essential ideas behind two distinct theoretical mechanisms which yield memory
loss of non-Gaussianity in the system. The first of the two mechanisms is an
instance of those studied in earlier theoretical works, while the second one
has been introduced in Ref. [24]. We will then lay out in more detail the
precise experimental observations presented in Ref. [24]. In Section 2 the two
mechanisms, which can be considered as candidates to explain the experimental
observations, will be presented in detail. Each of the two mechanisms relies
on the presence of different properties of the initial states and of the
dynamics. The characteristics of the initial states in the experiment relevant
to corroborate the two candidate mechanisms will be presented in Section 3 and
those related to dynamical aspects of the experimental system will be
presented in Section 4. Further observations based on the experimental data
are discussed in Section 5. By combining the conclusions of all previous
sections we finally evaluate the role of each of the two mechanisms in
Subsection 5.3 and conclude that the one that is mainly responsible for the
emergence of Gaussianity in the experiment is the second. At the end of the
manuscript we provide conclusions and outlook in Section 6.

Our results are to a large extent based on the use of a specific reading of
quantum field tomography. In the present context, this term refers to an
indirect recovery of all second moments of quadratures [25] from an ensemble
of identically prepared quantum states, exploiting suitable time evolution
under non-interacting Hamiltonians. It is a recovery method reminiscent of
full quantum state recovery based on data [26]. After all, a quantum state of
a quantum many-body system can be characterized by the collection of all
moments of observables [22]. It is important to stress that the experimental
prescription allows for the reconstruction of higher moments of quadratures as
well [22], an insight that allows us to monitor the process of Gaussification
in time. The tomographic method allows us to verify that the conditions of the
second Gaussification mechanism are satisfied to a sufficient degree in the
experiment by direct analysis of the experimental data (Subsection 3.3.2).
Preliminary evidence for the validity of these requirements has been presented
in Ref. [24] using simulations based on the classical fields approximation
[27].

In the course of our analysis, we derive further results of independent
interest, including an analysis of the scaling of correlations in the initial
equilibrium states (Subsection 3.1), a demonstration of light-cone propagation
of correlations (Subsection 4.1) and a study of the dynamics of full counting
statistics (Subsection 5.1). Lastly, we summarize the theoretical description
of the experimental system by means of TLL theory in the Appendix, justifying
why deviations from this theoretical model are negligible based directly on
experimental observations. Overall, our work provides a detailed comparison
between theory and experiment that aims to contribute to our understanding of
the physics of one-dimensional coupled atomic condensates and their dynamics.

###  1.1 Two mechanisms for the emergence of Gaussian correlations

Figure 1: After an interaction quench in an isolated system, the memory of
interactions as encoded in the initial state through non-Gaussian correlations
can gradually decay by means of two distinct mechanisms (broadly introduced in
Subsection 1.1 and discussed in detail in Section 2): a) Gaussification by
spatial scrambling (Subsection 2.1) rests on two pillars. The first is a
property of the initial state to yield independent results in measurements
between distant points, which is briefly called clustering correlations. Its
validity in the experiment is discussed in Subsection 3.1. The second pillar
is a property of the dynamics to make initially local fluctuations spread over
large portions of the system during the time evolution, and its validity is
tested in Subsections 4.1 and 4.2. b) Gaussification by canonical
transmutation is the second mechanism discussed in this work (Subsection 2.2).
It relies on the initial state having at least one Gaussian canonical sector
(Subsection 3.2) that dominates the other sector which is non-Gaussian
(Subsection 3.3). In this mechanism, quadrature rotation during the dynamics
results in the initially non-Gaussian quadrature turning into Gaussian. In
coordinate space, this gives rise to dephasing and rephasing dynamics which we
will broadly refer to as canonical transmutation dynamics.

We aim to precisely understand what happens when a system lies initially in a
non-Gaussian state and subsequently evolves under Gaussian dynamics. Any
process wherein non-Gaussian correlations in an isolated quantum system decay
over time rendering the resulting state effectively Gaussian will be called in
this work ‘Gaussification’.

In nature, we often find that a description of the system by means of
statistical mechanics is instructive and accurate. For this to be true when
starting out of equilibrium, some sort of scrambling of information encoded in
the initial state must occur during the dynamics in one way or the other. We
expect this to occur generically, under interacting or weakly interacting
dynamics, in accordance with classical intuition built around Boltzmann’s
$H$-theorem and phenomenology surrounding kinetic equations. Such scrambling
of information and memory loss effects in isolated quantum systems can be
induced by simple Gaussian dynamics. The induced memory loss can be sufficient
to give rise to an agreement between local reductions of the unitarily evolved
density matrix and the respective marginals of a Gaussian steady state [5, 6,
14, 16, 17, 18, 7, 19].

One such process of Gaussification – as we will argue – occurs in conjunction
with a notion of what we call spatial scrambling. Intuitively, in a large
system as the elapsed time becomes long, local observables depend on larger
and larger amounts of incoherent initial information originating from distant
points. If distant points have initially been only weakly correlated then this
results in an elimination of non-Gaussian features in correlation functions.
The unitary Gaussian dynamics implements it in a way that is arguably in
reminiscence of classical central limit theorems [14, 7]. In fact, a precise
connection to mathematical proofs of Lindeberg central limit theorems on the
level of characteristic functions can be established [14].

The mechanism of Gaussification by _spatial scrambling_ rests on two essential
physical ingredients. Firstly, the initial correlations of the effective
fields in terms of which the dynamics is Gaussian must satisfy the condition
of clustering, i.e., their correlations between distant points must factorize
as for independent variables, or, equivalently, their connected correlations
must decay with distance. The real-space distance scale for this is commonly
set by the correlation length. Weaker conditions (like algebraic instead of
exponential clustering or sub-extensivity of the initial fluctuations of
macroscopic conserved quantities [21]) can also play the same role, with the
main physical condition on the initial state being that it has local
characteristics like the equilibrium states of systems with local
interactions. Such conditions are quite ubiquitously valid in nature. This is
important because statistical mechanics has wide applicability so the
prerequisites for its emergence in isolated systems should be broadly
fulfilled. Secondly, the dynamics must induce delocalization, i.e., an
initially localized fluctuation of the effective field must spread with time,
not remain close to its original position or just rigidly move through the
system during the evolution. In Gaussian dynamics such a behavior is quite
typical and can be linked to properties of the energy dispersion relation,
specifically its non-linearity. These two broadly valid conditions in essence
imply a Gaussification process [16, 17, 18, 5, 7, 19, 20]. Fig. 1 provides an
illustration of these ‘pillars’ and summary for reference during the reading
of this work, including information on which section discusses each of the
‘pillars’.

While the pillars on which Gaussification by spatial scrambling rests are
quite general, their applicability cannot be taken for granted. A particularly
interesting case where this is not given a priori occurs when the collective
fields providing the Gaussian effective description of the system’s dynamics
are non-locally related to the physical local degrees of freedom (e.g.
particles), or if the spectrum of the effective fields is to a good
approximation linear. A non-local relation between local physical fields and
the collective fields of the effective description might mean that the initial
clustering condition is not guaranteed for the latter, even if valid for the
former: indeed, while clustering is a typical physical requirement for local
fields and observables, it does not have to be satisfied by non-local
collective fields. Moreover, when the dispersion relation of the collective
field excitations is linear, then, even if these fields are genuinely local,
the dynamical delocalization condition is broken because the time evolution
does not induce spreading of initially localized wave-packets which travel
instead pinned at two moving points. Both these aspects in question are
directly motivated by the experimental study of a quench of the effective
interaction between phononic excitations in coupled one-dimensional ultra-cold
gases.

Specifically, in Ref.[24], a decay of non-Gaussian correlations in time has
been observed and a novel mechanism for its explanation has been proposed, a
mechanism which we will call in this work Gaussification by _canonical
transmutation_. This mechanism is at play when the initial state capturing the
correlations of two canonically conjugate fields yields Gaussian correlations
for one canonical variable and non-Gaussian for the other. A simple harmonic
phase-space rotation of the system’s eigen-modes after switching off the
interaction may then lead to the decay of non-Gaussianity due to dilution of
the non-Gaussian into the Gaussian component if the latter dominates in the
mixing process. Such a change in the internal make-up of an object changing
drastically its overall properties seems to agree with the general meaning of
the word ‘transmutation’, so for lack of a better term we will consistently
employ it in this work to make clear which of the two mechanisms we are
referring to. As we will discuss in detail later, this mechanism rests on
three ‘pillars’ as summarized and illustrated in Fig. 1, and the individual
ingredients will be discussed one-by-one in the following. Again, we speak of
pillars in the sense that they seem to be necessary for the effect as the
absence of each one of them breaks down the mechanism and leads to
preservation of the memory of non-Gaussian correlations.

###  1.2 Experimental observation of the decay and recurrence of non-Gaussian
correlations

Having introduced the two theoretical mechanisms of Gaussification, our goal
is to investigate whether they can explain the observed decay of non-Gaussian
correlations in the experiment of Ref. [24]. Let us look into the context and
findings of this experiment in more detail. We will give a description of the
experiment and of the analysis of measured data and summarize the main
results.

In principle, an overall intuition regarding the main experimental observation
can be obtained based on a purely statistical consideration: The experiment
yields outcomes which differ from one experimental realization to another. The
outcomes of measurements are treated as instances of random variables so that
from this sample estimates of statistical moments of these variables can be
extracted. If the fourth and higher order cumulants are negligible then we
speak of a Gaussian state; in the opposite case we have a non-Gaussian state.

The experimentally measured variables (phase and density of the atomic gas)
constitute collective degrees of freedom in an effective description of the
actual many-body system. By varying a parameter of the system and measuring
the statistical moments when the system is at equilibrium, we can study how
the parameter controls the non-Gaussianity of equilibrium states in the
effective description. From the observation that the size of non-Gaussianity
depends on the parameter under consideration, just from statistical
considerations we conclude that a non-linear _effective interaction_ must be
present and can be controlled by the above parameter.

In addition, by rapidly changing the parameter and measuring the non-
Gaussianity at various subsequent times, we can study how this changes over
time. The experiment has displayed a decay of non-Gaussianity as a function of
the time elapsed after switching the parameter from the non-Gaussian to the
Gaussian regime. This may sound quite reasonable intuitively: the system
dynamically adapts to the change in the external parameter relaxing to the
corresponding equilibrium state. However, one aspect of the dynamics studied
in Ref. [24] evades a simple intuitive interpretation. How should one
interpret the fact that after a monotonous decay of non-Gaussianity the
experiment shows a revival of the non-Gaussian correlations at a later time?
To answer this question we first need to know a crucial piece of information
about the dynamics. To a very good approximation the system evolves as a
_closed_ one: interaction with the environment is strongly suppressed over the
time scale of the experiment, so that the system is practically isolated.
Under such settings the _information_ about the system being initially non-
Gaussian is never lost but gets first hidden from view, resulting in a
Gaussian state at intermediate times, and is subsequently retrieved again at
the time of the revival of non-Gaussianity. Hence, given that the dynamics is
essentially unitary, the information about the initial state has never left
the system and has not been irreversibly scrambled. The emergence of a
Gaussian equilibrium-like state is now somewhat less obvious and the mechanism
leading to it deserves investigation.

Figure 2:  a) Illustration of the experimental setup: A gas of ultra-cold
atoms confined into a one-dimensional trap can be split in two adjacent
parallel traps controlled by an atom chip. By controlling the barrier height
of the transverse double well potential which plays the role of a Josephson
junction, one can engineer a many-body interaction of the cosine type. Ramping
up the barrier height amounts to switching the interaction off. b) Preparing
the gas at an equilibrium state and releasing it from the trap, the
interference between the two atomic clouds reveals information about the
spatial profile $\phi(z,t)$ of their relative phase at the time of the
release. Repeating the experiment many times results in an ensemble of phase
profiles from which one can extract the probability distribution of various
physical observables derivable from the phase field. In the presence of
interaction and at sufficiently high temperature, the distribution of the
relative phase is highly non-Gaussian, exhibiting a central peak at $\phi=0$
and accompanied by side peaks at $\phi=\pm 2\pi$ with a lower height. This
pattern clearly indicates the presence of solitons in the phase profiles and
the concentration of the phase at the minima of the cosine potential. If the
coupling $J$ is very large, the interaction potential is approximately a steep
parabola and phase fluctuations are Gaussian. Changing the interaction
strength $J$ from a nonzero value to zero and monitoring the subsequent time
evolution of the phase profiles, we observe the transition from an initial
non-Gaussian to an emergent Gaussian distribution. The plots show typical
examples of equal-time ensembles of phase profiles and the corresponding
histograms of phase values in all points $z$, as measured in the experiment.
The three cases correspond to two different types of initial states (at
intermediate or large $J$) and a typical state that emerges after quenching to
$J=0$. (c) Illustration of the time evolution of the phase distribution.

We will now proceed to explaining how this physical picture has been realized
in the experiment. A description of the experiment and data analysis is
illustrated in Fig. 2. The system considered consists of two parallel and
adjacent one-dimensional quasi-condensates of a few thousand rubidium atoms
trapped, cooled and controlled by means of an _atom-chip_ setup (see Fig.
2.a). The experimental measurements are interference pictures generated when
the two-component gas is released from the trap and let to expand freely. From
the interferometric measurements we extract phase profiles $\phi(z,t)$
corresponding to the relative phase between the two components of the gas as a
function of the longitudinal dimension $z$ at the time of the release $t$. By
repeating the experiment many (several hundreds of) times we obtain an
ensemble of phase profiles from which we can derive statistical measures of
the phase field distribution: correlation functions between different points,
moments, cumulants, or the full distribution function of the phase.

The relative phase profiles play the role of the fundamental random variables
in our statistical consideration. In practice, a phase profile is a vector
whose entries correspond to the estimate of the relative phase between the two
quasi-condensates at the position corresponding to a given pixel of the read-
out camera [28, 29, 30, 31, 32, 22, 33]. The phase observable is modelled by a
bosonic field $\hat{\varphi}(.)$ ranging over the extension of the one-
dimensional quasi-condensate, whose excitations have the physical
interpretation of phonons.

The two quasi-condensates are trapped in a transversal double well as
illustrated in Fig. 2.a with a barrier that controls the tunneling strength
$J$ between the two wells. This tunnel coupling is quite well characterized:
For intermediate values of $J$ (relative to the temperature) it leads to an
effective many-body interaction of the phonons and hence the state of the
system that is prepared at equilibrium at such couplings is non-Gaussian. This
interaction term is known to agree well with a _sine-Gordon potential_ , as
predicted in Refs. [34, 35] and verified experimentally in Ref. [22]. The
latter work has also shown that in the limit cases of small or large $J$ a
Gaussian distribution is obtained. One simple way to assess the non-
Gaussianity of the system’s state is to disregard the spatial variation of the
phase profiles in the $z$-direction and consider a cumulative histogram. The
question is then how far the obtained distribution is from a Gaussian. Fig.
2.b illustrates two typical cases. In the case of intermediate $J$ discussed
above the histogram exhibits long tails or even shoulder peaks that signify
deviations from a regular Gaussian distribution. In the case of small or large
$J$ on the contrary the distribution is plain Gaussian.

Ramping up the barrier height of the double-well between the two adjacent
quasi-condensates, as illustrated in Fig. 2.a, reduces the tunneling of atoms
between the two wells, hence reducing the $J$ parameter to zero. Doing this
from the intermediate regime where interactions play initially a substantial
role corresponds to switching off the effective interaction, which is the
quench situation that we want to consider. The resulting dynamics of the
histograms is illustrated in Fig. 2.c: An initial distribution with shoulders
that constitute deviations from a Gaussian evolves so that the deviations
decay over time.

Fig. 2 is at this stage an illustration of a subset of questions one can study
experimentally using the atom-chip platform. We will next discuss a more
space-resolved way of assessing the non-Gaussianity of the system using
correlation functions that has been presented in Ref. [24]. Nonetheless, the
dynamical behavior illustrated in Fig. 2.c will make an appearance towards the
end of the manuscript in Subsection 5.1 where we will present how histograms
of the type illustrated in Fig. 2.b vary in time.

The measurement outcomes of the phononic phase field $\hat{\varphi}(.)$ can be
viewed as spatially resolved values of the relative phase between the two
quasi-condensates. At this point, it should be noted that given that the
relative phase is an angular variable its value at any point can only be
measured with an ambiguity of a $2\pi n$ shift where $n$ is an integer. The
measured phase profiles are derived by imposing the condition that the phase
at some reference point is within the interval $[-\pi,+\pi)$ and then
‘unwrapping’ the phase profile so that values at any two neighboring points
differ no more than $\pi$. The ensembles of phase profiles shown in Fig. 2.b
are extracted in precisely this way, where the reference point is fixed to the
middle. This procedure still involves however an arbitrary choice of an
overall $2\pi n$ phase shift. To completely remove this ambiguity we can
restrict our analysis to observables calculated strictly on phase differences
between any point and some (arbitrarily chosen) reference point. For this
reason, we define the observable phase difference with respect to the
reference point $z_{0}$ as

$\displaystyle\Delta\hat{\varphi}(z):=\hat{\varphi}(z)-\hat{\varphi}(z_{0}).$
(1)

By repetition of the experiment and interferometric measurement under
identical conditions, the full correlation function can be estimated as

$\displaystyle\Phi(z_{1},\ldots
z_{n})=\left\langle\prod_{i=1}^{n}\Delta\hat{\varphi}(z_{i})\right\rangle\ .$
(2)

This is the statistical moment that we have defined above. Taking $n=1$ gives
just the average phase that typically vanishes due to effective symmetries of
the physical processes involved in the state preparations. For $n=2$, we
obtain the second moments which allow us to parametrize a Gaussian
distribution and to build higher-order correlations using Wick’s theorem.

Noticing that correlation functions of odd order are negligible, we find that
for $n=4$, the _connected_ correlation functions take the form

$\displaystyle\Phi_{\text{\rm con}}(z_{1},\ldots z_{4})=\Phi(z_{1},\ldots
z_{4})-\Phi_{\text{Wick}}(z_{1},\ldots z_{4})\ ,$ (3)

where we just subtract the standard _Wick_ decomposition

$\displaystyle\Phi_{\text{Wick}}(z_{1},z_{2},z_{3},z_{4})=~{}$
$\displaystyle\Phi(z_{1},z_{2})\Phi(z_{3},z_{4})+\Phi(z_{1},z_{3})\Phi(z_{2},z_{4})+\Phi(z_{1},z_{4})\Phi(z_{2},z_{3}),$
(4)

an expression in which the sum ranges over all possible pairings of indices.
We see that such a 4-point connected function vanishes for all Gaussian
distributions because then $\Phi(z_{1},\ldots
z_{4})=\Phi_{\text{Wick}}(z_{1},\ldots z_{4})$ essentially by definition.

Such an analysis, i.e., extraction of various moments and evaluation of
connected functions, depends only on the measured phase profiles and hence
until now we did not indicate any time dependence. The time dependence is
encoded in the quantum state at the time of the measurement from which the
phase profiles are being sampled. That is, each time $t$ we measure, after
quenching the tunneling parameter $J$, the system is described by some density
matrix $\hat{\varrho}(t)$ and the expectation values $\langle\cdot\rangle$ in
Eqs. (2-4) refer to this density matrix. When the system involves many degrees
of freedom that interact with each other, then its state is quite complex and
it is not practical to inquire about its entire quantum state
$\hat{\varrho}(t)$. Instead, one should consider correlation functions and
indeed these can be obtained experimentally and be used to characterize the
system at different times $t$ [22, 23]. Accordingly, we will indicate
correlation functions obtained at different times by writing explicitly their
dependence on first the spatial variables and then the time variable, as
$\Phi(z_{1},\ldots,z_{n},t)$.

Figure 3: Observation of Gaussification (_left_) and the recurrence of non-
Gaussianity (_right_) in the experiment. Plots of the non-Gaussianity measure
$M^{(4)}$ as a function of time in various experimental scans corresponding to
different initial states, as presented in Ref. [24] (the plots are based on
the data published in Ref. [36]). The time window ranges from the initial time
until 16 ms _(left)_ or 32.5 ms _(right)_. Notice that the decay of non-
Gaussianity follows a relatively sharp decrease to a small saturation value.
At the recurrence time [30] the non-Gaussianity returns to quite a high value
relative to the initial one. The green shaded area indicates the bias values
of $M^{(4)}$ corresponding to Gaussian states with finite statistics.

We are now in a position to discuss the measure for the non-Gaussianity of the
phase fluctuations based on correlation functions that has been presented in
Refs. [22, 24]. It is given by

$\displaystyle
M^{(4)}(t)=\frac{S_{\text{con}}^{(4)}(t)}{S_{\text{full}}^{(4)}(t)}=\frac{\sum_{z_{1},\ldots,z_{4}}|\Phi_{\text{\rm
con}}(z_{1},\ldots z_{4},t)|}{\sum_{z_{1},\ldots,z_{4}}|\Phi(z_{1},\ldots
z_{4},t)|}\ .$ (5)

This quantity vanishes for a Gaussian state and meaningfully quantifies its
non-Gaussian character. Note that given that $M^{(4)}$ as defined above is a
non-negative quantity, in an actual experiment it can only go as low as a
small non-vanishing value, reflecting the experimental statistical
fluctuations (finite size of the statistical sample). The summation window is
typically taken over a region around the middle of the system where the atomic
gas is to a good approximation homogeneous and the data is more reliable than
closer to the edges. In all scans of the present study the trap used to
confine the atoms is box-like, more specifically, a superposition of a
harmonic trap with a box trap of much smaller extent so that the density is
approximately homogeneous over the entire system. Therefore, unlike for the
more commonly used harmonic traps, there is no significant inhomogeneity in
our experimental system, except for boundary effects that are still present in
the close vicinity of the edges of the box. As a technical remark, in Fig. 3
for scans 1-4 the box-like trap is $75\text{\,}\mathrm{\SIUnitSymbolMicro m}$
long and we analyze using Eq. (5) the central region of length
$50\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and for scans 5-7 we analyze the
central $38\text{\,}\mathrm{\SIUnitSymbolMicro m}$ region of a
$50\text{\,}\mathrm{\SIUnitSymbolMicro m}$ long trap.

The experiment has observed the decay of this measure for all initial
conditions (and trap geometry), as shown in Fig. 3 presenting the dynamics of
$M^{(4)}$ for various initial states. The quantity $M^{(4)}$ decreases rapidly
to a small value that is indistinguishable from that of a Gaussian state. Note
that, because $M^{(4)}$ is by definition a non-negative quantity and because
of the statistical fluctuations in a finite sample of measurements, the mean
value of $M^{(4)}$ corresponding to measurements on a Gaussian state is not
zero but positive: this is called ‘finite statistics bias’. Interestingly, for
a box-like trapping geometry, at a particular later time, $M^{(4)}$ becomes
large again. This revival of non-Gaussianity has been discussed above and can
be fully accounted for by the mechanism of Gaussification by canonical
transmutation which will be introduced in more detail in the following
section.

Figure 4: Characterization of the initial states of the experimental scans of
Ref. [24]. Plot of the observed value of $M^{(4)}$ in the initial equilibrium
state as a function of the observed coherence factor $\langle\cos\phi\rangle$
as first explored in Ref. [22]. The initial state of scan 1 and 7 corresponds
to large $J$, i.e., the Klein-Gordon (KG) regime and is Gaussian, while those
of scans 2 to 6 are non-Gaussian, with the highest value of non-Gaussianity
reached for scan 2. For small $J$ the system is in the Tomonaga-Luttinger
liquid regime (TLL). The green curves correspond to theoretical simulations of
the sine-Gordon model at thermal equilibrium ($\lambda_{\mathrm{T}}$
represents the thermal coherence length) based on the classical field
approximation [27].

A characterization of the initial states based on the strength of the
interactions is presented in Fig. 4 where the observed value of $M^{(4)}$ in
the initial equilibrium state is plotted against the observed coherence factor
$\langle\cos\phi\rangle$. The latter is controlled by the interaction coupling
$J$ and unlike that it is directly measurable in the experiment, so that it
can play the role of the control parameter for the strength of interactions.
For $J$ increasing from zero to a relatively large value the equilibrium
coherence factor changes from zero to one. It should be mentioned that in the
experiment $J$ cannot be chosen to be larger than $\hbar\omega_{\text{tr}}$
where $\omega_{\text{tr}}$ is the transverse trap frequency, otherwise
transverse mode excitations are present. In this case the system can no longer
be considered as effectively one-dimensional. The initial state of scan 1
corresponds to the heavy mass regime of the sine-Gordon model (KG regime)
where $\langle\cos\phi\rangle\approx 1$ and is therefore Gaussian. Quenching
to the TLL regime results in a highly squeezed initial state due to the fast
change of the effective mass parameter from a large value to zero. Scans 2, 3
and 4, on the other hand, correspond to non-Gaussian states at decreasing
values of the interaction strength $J$. The initial non-Gaussianity reaches
its highest observed value in scan 2 and progressively decreases in scans 3
and 4. A theoretical explanation of the observed $M^{(4)}$ versus
$\langle\cos\phi\rangle$ curve is given by the classical field approximation
of the sine-Gordon model [22] using a stochastic method developed in Ref. [27]
whose results are consistent with numerical simulations of the actual quantum
model [37].

##  2 Theoretical discussion of two readings of Gaussification

In what follows, we will consistently refer to _Gaussification_ when we speak
about the memory loss effect wherein non-Gaussianity decays in an isolated
system following an interaction quench in time. It is a main theoretical
insight of this work that there are basically two readings of this effect,
which we will study and complement with each other in detail in this section.
The actual experimental findings from cold atomic quantum field systems will
later be put into the context of these Gaussification mechanisms, and we will
develop and elaborate on what the actually dominant effect in the experiment
is. In what follows, we do not aim at providing a mathematically rigorous
framework of Gaussification as presented in Refs. [5, 14, 7, 19], but instead
keep the discussion on an intuitive, physically minded level.

###  2.1 Gaussification by spatial scrambling

The first mechanism, and this is what this section will focus on, can be seen
as being Gaussification by spatial scrambling. The basic idea to understand is
that: Correlation functions of fields delocalized over regions much larger
than the correlation length are effectively Gaussian. Expanding on this
statement is the goal of this subsection. We keep the discussion close to the
experimental setting at hand. In this context, let us consider the field in
the above statement to be the particle velocity field, defined as (see
Appendix A)

$\displaystyle{\hat{u}}(z)=\partial_{z}\hat{\varphi}(z)\ .$ (6)

This is a local field, hence one would expect that its correlations cluster in
typical initial conditions. This is the first requirement of Gaussification by
spatial scrambling constituting one of the pillars in Fig. 1.

######  Condition 1 (Clustering of correlations).

We assume the initial correlation functions to be exponentially decaying in
the distance as

$\displaystyle\left|\left\langle\prod_{i=1}^{n}{\hat{u}}(z_{i})\right\rangle_{\text{con}}\right|\leq\text{const}\times\mathrm{e}^{-\tfrac{1}{2}\sum_{i,j=1}^{n}|z_{i}-z_{j}|/\xi}\quad(\text{for
}z_{i}\neq z_{j})\ .$ (7)

States that satisfy this bound will be referred to as exhibiting exponentially
clustering correlations.

This nomenclature is motivated by the fact that connected correlation
functions are substantial only for $z_{i}$’s that cluster together within a
range of the order of the correlation length $\xi>0$. This property is
ubiquitously valid in many systems, specifically in ones prepared in ground
states of gapped quantum systems and for those close to thermal equilibrium at
sufficiently high temperatures. Whenever all the positions $z_{i}$ are
sufficiently far from each other

$\displaystyle\left\langle\prod_{i=1}^{n}{\hat{u}}(z_{i})\right\rangle\approx\left\langle\prod_{i=1}^{n}{\hat{u}}(z_{i})\right\rangle_{\text{Wick}}\
.$ (8)

For general lattice models with finite-dimensional constituents, it has indeed
been proven that not only ground states [38], but also high temperature states
exhibit a clustering of correlations [39, 40]. Thus, initial states exhibiting
exponentially clustering correlations can be strongly correlated states within
the range of the correlation length, while exhibiting approximately Gaussian
correlations at long distances. The correlations for points separated beyond
the correlation length may intuitively be thought of as playing the role of a
Gaussian ‘bath’ in the sense that Gaussification occurs as the result of the
dynamical mixing of the initially non-Gaussian component of correlations into
the much larger Gaussian component.

The delocalization of the field can result in Gaussianity both at long and
short distances as compared to the correlation length. Let us see how this
works by means of a simple example. Consider the field integrated over a
region of the system $R$

$\displaystyle\overline{\Delta\hat{\varphi}}=\frac{1}{\sqrt{|R|}}\int_{R}\mathrm{d}z~{}{\hat{u}}(z)\
.$ (9)

Here we weigh this expression by the size $|R|$ of the region. This is similar
to considering independent identically distributed variables $X_{i}$ and
forming the central limit variable $\overline{X}=n^{-1/2}\sum_{i=1}^{n}X_{i}$.
Classically, the inverse-square-root normalization of $\overline{X}$ is
crucial and then the distribution becomes Gaussian in the limit
$n\rightarrow\infty$. We can estimate the connected part of this observable as

$\displaystyle\left\langle\overline{\Delta\hat{\varphi}}^{\,4}\right\rangle_{\text{con}}\approx\frac{1}{|R|^{2}}\int_{R^{4}}\mathrm{d}^{4}z\left\langle\prod_{i=1}^{4}{\hat{u}}(z_{i})\right\rangle_{con}\approx\text{const}\times\frac{\xi^{3}}{|R|}$
(10)

where the constant is bounded by the maximum of the correlations at short
distances determined by the UV cut-off. We hence see that if we take a field
whose correlations exponentially cluster, then the average over a large region
results in connected functions scaling inversely proportional to the size of
that region. Thus in the limit of $|R|\rightarrow\infty$ the connected part
vanishes.

A type of such delocalization can be implemented by unitary Gaussian dynamics
of an isolated quantum system [16, 17, 18, 5, 7, 19, 20]. This is the second
pillar of spatial scrambling in Fig. 1.

######  Condition 2 (Delocalizing dynamics).

The dynamics generated by the non-interacting Hamiltonian is expected to be
delocalizing. Said in different words, the propagator of an initially local
field decays with time in the entire space.

Let us expand on this condition in more detail. For a Hamiltonian that is non-
interacting i.e. quadratic in terms of a set of fields, the Heisenberg
equations of motion that determine the unitary time evolution of these fields
are linear differential equations with respect to time. As such they can be
solved for general initial conditions and the time evolved fields are
expressed as linear combinations of the initial fields. Conversely, if the
unitary time evolution is of this linear form then the governing Hamiltonian
must be non-interacting. In the present context, for example, the Heisenberg
equation of motion satisfied by
$\hat{\varphi}(x,t)=\hat{U}^{\dagger}(t)\hat{\varphi}(x,0)\hat{U}(t)$ with
$\hat{U}(t)=\mathrm{e}^{-\mathrm{i}Ht}$ is a linear partial differential
equation of second order with respect to time, or equivalently a linear system
of two first order equations for the relative phase field $\hat{\varphi}(x,t)$
and its canonically conjugate field $\delta\hat{\varrho}(x,t)$ which expresses
the relative density fluctuations of the two-component atomic gas (see
Appendix B for details). The solution to the associated initial value problem
can be expressed in terms of $\hat{\varphi}(y,0)$ and
$\partial_{t}\hat{\varphi}(y,0)\propto\delta\hat{\varrho}(y,0)$ by means of
the Green’s function method, which finally gives the field ${\hat{u}}(x,t)$
expressed as

$\displaystyle{\hat{u}}(x,t)$
$\displaystyle=\int\mathrm{d}{y}~{}G_{u}(x,y,t){\hat{u}}(y,0)+\int\mathrm{d}{y}~{}G_{\rho}(x,y,t)\delta\hat{\varrho}(y,0)\
.$ (11)

Here $G_{u}(x,y,t)$ and $G_{\rho}(x,y,t)$ are the retarded Green’s functions,
commonly called propagators, corresponding to the equations of motion for
${\hat{u}}(x,t)$ under the initial conditions
${\hat{u}}(y,0)=\delta(x-y),\delta\hat{\varrho}(y,0)=0$ and
${\hat{u}}(y,0)=0,\delta\hat{\varrho}(y,0)=\delta(x-y)$, respectively. They
encode the dynamics triggered by the initial conditions and express the
impulse response of the field ${\hat{u}}(x,t)$ coming from a localized
disturbance of ${\hat{u}}(y,0)$ or $\delta\hat{\varrho}(y,0)$, respectively.
In practice, they can be calculated by expanding in the eigen-modes of the
Hamiltonian (Fourier modes in the translationally invariant case).

If the dynamics admits a bound of the form

$\displaystyle|G_{u/\rho}(x,y,t)|\leq\text{const}\times t^{-\alpha},$ (12)

for some $\alpha>0$, then we say that it features delocalization, as
summarized in Fig. 1. For example, in a typical lattice system with local
dynamics the propagators are upper bounded as above with an exponent
$\alpha=1/2$. The causality of the time evolution (Lieb-Robinson bound) is
implemented in $G$ as exponential decay with the distance outside the light-
cone region $|x-y|<ct$, where $c$ is the speed of sound. The light-cone
propagation bound is also a characteristic of the low-energy effective theory
describing the present continuous system. These properties imply then that
${\hat{u}}(x,t)$ equals to an effective average (with an oscillatory weight
function) of ${\hat{u}}(y,0)$ within the light-cone. This is similar to the
above discussed case with $|R|=ct$ and so we expect a similar bound on the
connected part of the correlations

$\displaystyle|\langle{\hat{u}}(x,t)^{4}\rangle_{\text{con}}|\approx\text{const}\times\frac{1}{ct}\
.$ (13)

We hence encounter the following scenario: If a system exhibits dispersive
dynamics, then an initially local field spreads in space under the time
evolution and its connected correlation functions should then decay at least
as a power-law in time. If instead the dispersion relation is linear, then
information about initial correlations travels coherently and non-Gaussianity
may be preserved.

It should be remarked that the dynamics of lattice models does not implement
precisely a square-root decay of the Green’s function in the entire space [5,
7, 19, 20]. Rather, as can be seen, e.g., by the stationary phase
approximation, there is a wave-front at the edge of the effective causal cone
which decays more slowly than is true for the interior of the cone. This
complication is one of the reasons why a complete derivation of Gaussification
by spatial scrambling goes beyond the intuitive picture sketched here.

More specifically, in the above we have considered for illustration purposes a
typical case of dynamics characterized by an effective light-cone spreading.
In non-relativistic continuous systems, however, correlations can in principle
spread quite fast with no effective maximum velocity bound. Even in this case
where field spreading is not constrained the Gaussification by spatial
scrambling arguments may still be applicable [17]. In the illustrative
discussion of the mechanism we assumed the scaling exponent to be everywhere
in space equal to $1/2$ which implies uniform decrease of the Green’s function
inside the light-cone. This is instructive for a general discussion but again,
the delocalization mechanism is valid more generally. In particular, it holds
for dynamics generated by local Hamiltonians where the Green’s function scales
according to the $1/2$ exponent in the bulk of the effective light-cone but at
the wave-fronts at its edge the asymptotic scaling in time is $1/3$ [19, 17,
20]. In any case the essential arguments available in the literature in all
cases can be argued to rely closely on the delocalization of the fields as
explained here (see Ref. [20] for an accessible discussion of these
complications along the lines that we presented here and Ref. [19] for a
general derivation of such types of scaling for local lattice models).

###  2.2 Gaussification by canonical transmutation

Gaussian correlations can also emerge based on a mechanism which does not
involve spatial scrambling. Instead of that, we want to consider the
possibility of the decay of non-Gaussian correlations that results from the
internal dynamics in each of the independent harmonic modes.

To be specific, let us consider the dynamics generated by a general non-
interacting Hamiltonian

$\displaystyle\hat{H}=\tfrac{1}{2}\sum_{k=1}^{\infty}\hbar({\delta\hat{\rho}}_{k}^{2}+\omega_{k}^{2}{\hat{\phi}}_{k}^{2})$
(14)

with canonical commutation relations
$[{\delta\hat{\rho}}_{k},{\hat{\phi}}_{k^{\prime}}]=\mathrm{i}\delta_{k,k^{\prime}}$
for $k,k^{\prime}=1,2,\ldots$ and $\delta_{k,k^{\prime}}$ is the Kronecker
symbol. The time evolution of each harmonic mode under this non-interacting
Hamiltonian is

$\displaystyle{\hat{\phi}}_{k}(t)=\cos(\omega_{k}t){\hat{\phi}}_{k}(0)-\frac{\sin(\omega_{k}t)}{\omega_{k}}{\delta\hat{\rho}}_{k}(0)\
.$ (15)

Here, we find that the density sector is being rotated into the phase sector
and vice versa, which is what we mean by canonical transmutation: In
particular, for $t=\pi/2\omega_{k}$ we have that
${\hat{\phi}}_{k}(t)\propto{\delta\hat{\rho}}_{k}(0)$ with a complete
transmutation of the role of the operators in the canonically conjugate pair.
This is the crucial qualitative insight and captures the essential physical
process occurring after the quench, which allows us to account qualitatively
for the resulting Gaussification dynamics and at the same time is largely
independent of the model as even perturbed Hamiltonians would feature similar
rotation dynamics. The discussion of this quadrature rotation makes specific
the general idea of the necessary dynamical pillar as listed in Fig. 1.

######  Condition 1 (Quadrature rotation).

The dynamics generated by the non-interacting Hamiltonian is expected to be
implemented by a quadrature rotation to generate transmutation.

This relation is enough to argue that marginals of non-Gaussian states with a
certain structure of correlations are going to become Gaussian over time. To
give the basic idea, consider the case that initially only the phase sector is
non-Gaussian, e.g.,

$\displaystyle\langle{\hat{\phi}}_{k}(0)^{4}\rangle_{\text{con}}>0\ .$ (16)

At time $t=\pi/2\omega_{k}$ we have

$\displaystyle\langle{\hat{\phi}}_{k}(t)^{4}\rangle_{\text{con}}=\langle{\delta\hat{\rho}}_{k}(0)^{4}\rangle_{\text{con}}=0\
.$ (17)

We hence see that the connected correlation function of the phase mode $k$
will become fully Gaussian under this transmutation. Here, we see that if
there is such a structure of Gaussian correlations being contained in one
sector of the canonically conjugate pairs, then we can identify them again as
constituting a Gaussian bath. As suggested by this illustrative case, we wish
to promote this property (one canonical sector being Gaussian and the other
non-Gaussian) to a feature that can also be present in other systems with
possibly different degrees of freedom. For this reason we single out this
characteristic of the initial state as a distinct pillar of the mechanism of
Gaussification by canonical transmutation as listed in Fig. 1.

######  Condition 2 (Existence of a Gaussian canonical sector).

The initial state should have a canonical sector that is Gaussian and
decoupled from the non-Gaussian sector.

This basic observation so far phrased for essentially the non-local eigen-
modes of the system is at play in more complicated situations. First, what
happens at intermediate times? Assuming the Gaussianity of one canonical
sector, we see that the non-Gaussianity of the other sector of each eigen-mode
will oscillate between the initial value and zero. Secondly, such
transmutation dynamics within the eigen-modes will lead to more complicated
dynamics in real-space. For this let us take the instructive example of the
mode decomposition of the fields in a homogeneous system of length $L$

$\displaystyle\hat{\varphi}(z)\sim\sum_{k}\cos(\pi zk/L){\hat{\phi}}_{k}\ .$
(18)

Using Eq. (15), we see that the dynamics of the real-space field
$\hat{\varphi}(z,t)$ will lead to a transmutation into the density sector in
each mode at its own frequency. For this reason the transmutation dynamics
will turn out to be a mixture of some eigen-modes which happened to have
rotated significantly and some which did not, resulting overall in a reduced
value of non-Gaussianity for almost all times. The analysis of mode mixing is
more generally relevant in the study of signatures of space-time propagation
of phase correlations and has been studied in Refs. [41].

######  Condition 3 (Dominance of the Gaussian canonical sector).

The initially Gaussian canonical sector should be much larger than the non-
Gaussian one.

As we will discuss later, the decay of Gaussianity can be substantial if the
Gaussian bath is at high temperature and there is an imbalanced energetic
penalty which leads to strong squeezing, i.e.,
$\langle{\delta\hat{\rho}}_{k}^{2}\rangle\gg\langle{\hat{\phi}}_{k}^{2}\rangle$.
In this case even if the dynamics is weighting both eigen-mode operators
similarly, the Gaussian correlations will dominate, because there will be more
contribution from fluctuations of the Gaussian type at any moment. When this
condition holds, after the quench one would see that fluctuations overall
increase with an increasing portion of these correlations being Gaussian. This
feature is an important aspect of the effect in the experiment in Ref. [24]
and will be dissected in a number of sections below. However, the discussion
geared towards the experiment can be viewed as being a feature that should be
recognized in its own right as potentially playing a role also in other
systems. For this reason, we include this feature as one of the pillars of
Gaussification by canonical transformation and refer to it as ‘dominance of
the Gaussian canonical sector’ (Fig. 1.b). Let us also remark that, even
though the dominance requirement means that the ratio of initial fluctuations
of the Gaussian over the non-Gaussian sector should be much larger than unit,
as we will see below, we actually obtain a large decrease of non-Gaussianity
even for relatively moderate ratios.

###  2.3 Similarities and differences of the two mechanisms

The two mechanisms discussed in this work have in common the setting in which
they can come into play and the type of effect they cause. In both cases, we
are interested in a non-Gaussian state influenced by some interactions in the
initial Hamiltonian. After interactions have been quenched, the connected
correlations become less prominent over time. So both mechanisms can lead to
the emergence of Gaussian correlations. One of the most intriguing
similarities is that effectively in both cases ‘Gaussianity’ is not being
dynamically produced, but rather dynamically ‘uncovered’.

In the case motivated by the experiment, the sector of phase observables has
been non-Gaussian while the canonically conjugate operators have been
Gaussian. Then, any Hamiltonian that mixes the two sectors sufficiently
strongly (quadrature rotation pillar) can lead to Gaussification. Here,
however, the extent of Gaussification is given by the extent in which the
density correlations dominate over time the phase sector. If this is very
substantial then connected correlations will become less important and yield a
negligible relative contribution to the full correlation function. From Eq.
(15) it seems that the ultimate extent of Gaussification is given by the
average between the two sectors which is the result one obtains by a time
average over a period $T\sim\omega_{1}^{-1}$.

One can think of an analogy with two water tanks, one with a high and the
other with a low level of water so that after connecting them the levels
oscillate and even out. Here, the phase sector has a large degree of non-
Gaussianity while the density sector has a low level and both start mixing
after the quench such that the full correlation functions of phases become
mostly Gaussian due to the Gaussian density fluctuations that are dynamically
rotated in after the interaction quench. So here the density correlations play
the role of a Gaussian ‘reservoir’.

One can identify a type of such reservoir that is dynamically mixed in also in
the Gaussification by spatial scrambling. In this case, however, a field
becomes spread over a large portion of the system and as we saw earlier higher
order correlation functions become largely made of initial correlation
functions from positions separated by large distances. These are initially
Gaussian as the correlation function factorizes due to the assumption of
clustering of correlations. So far away correlation functions can be
identified to play the role of the Gaussian ‘reservoir’ and the spreading of
correlations in form of homogeneous spreading of wave-packets ‘couples’ the
local correlations to that reservoir such that local expectation values are
described by a Gaussian state to an increasingly better approximation.

Having said that, there are also crucial differences. In Gaussification by
canonical transmutation, the non-Gaussianity need not necessarily disappear
fully, but rather become only overwhelmed by rotating in other Gaussian
correlations. Here, the degree of Gaussification by canonical transmutation is
limited by how much densities dominate typically over time, and how Gaussian
they are in the first place. Via wave-packet spreading, if the system is
asymptotically large, the local observables can be captured arbitrarily well
by a Gaussian state, so in contrast to Gaussification by canonical
transmutation there is no limit to the convergence to Gaussianity.

Gaussification by canonical transmutation can occur when spatial scrambling
will fail to be present, as exemplified by the experiment in Ref. [24]. The
most important example where this happens is in systems characterized by a
linear dispersion relation, which will not exhibit spreading of wave-packets,
so the only possible type of Gaussification is the relative one. As described
below (Subsec. 3.2), the required condition that an entire sector of initial
correlations is decoupled and Gaussian is possible to hold in a broad class of
systems prepared in a high temperature state which is described by the
classical field approximation. Therefore, the effect is not limited to the
sine-Gordon model as e.g. the $\phi^{4}$ theory should also exhibit similar
equilibrium correlations. In the case of linear dispersion systems discussed
here, the different nature of the canonical transmutation compared to the
spatial scrambling mechanism translates into a different scaling of the
approach to equilibrium in a large system, as explained in Subsec. 5.2.
Another feature of the Gaussification by canonical transmutation as opposed to
spatial scrambling is that the overall amount of non-Gaussianity of connected
correlations should be preserved. Lastly, Gaussification by canonical
transmutation can occur without the condition of initial clustering which is
necessary for the Gaussification by spatial scrambling – as long as there is a
Gaussian ‘reservoir’ sector that is rotated in.

##  3 Characteristics of the initial correlations in the experiment

One overarching observation concerning the two Gaussification mechanisms
discussed in this work is the presence of Gaussian correlations in the initial
state in some form such that they have initially been inaccessible (since we
can only measure directly phase and not density correlations) but then come to
view via the dynamics. We will now focus on the different possibilities for
some type of correlations to be initially Gaussian. We will begin by analyzing
one of the two pillars (as summarized in Fig. 1.a) of Gaussification by
spatial scrambling, namely the question of clustering of correlations. The
next subsection will discuss a microscopic argument that pertains to the
experiment.

###  3.1 Clustering and scaling in initial states

In this section, we will be discussing various characteristics that pertain to
the scaling of correlations in the experimental system at thermal equilibrium.
The relative phase defined in Eq. (1), whose statistics is the direct
experimental observable, is not a local field in a strict sense. This is
because it is only differences between the phases at two different points (one
of which is the reference point) that have a direct physical interpretation.
As a result, the correlations of the field that corresponds to the phase
difference between two points do not decrease with the distance but instead
increase or saturate, due to the cumulative effect of fluctuations in the
intermediate spatial interval. This is in contrast to the typical behavior of
the correlations of local fields which decrease with the distance when we
consider equilibrium states of local systems. The long-range character of
phase correlations has been discussed already in earlier works, see, e.g.,
Refs. [42, 43, 44, 41, 22]. Here we explain this behavior in an alternative
way (by deriving the scaling of correlations of the phase field from that of
its spatial derivative) and argue that in the non-Gaussian states considered
here the growth of phase correlations with the distance is related to the
presence of topological excitations in the system.

The space and time derivatives ${\hat{u}}=\partial_{z}\hat{\varphi}$ and
$\partial_{t}\hat{\varphi}\propto\delta\hat{\varrho}$ of the quantum phase are
expected to be local fields. Connected correlations of such local fields in
local quantum states are expected to decay with the distance: In ground states
of gapped Hamiltonians the decay of local field correlations is exponential
with a correlation length that is controlled by the mass of the lightest
quasi-particle excitation [45, 46]. In particular, this behavior holds for the
sine-Gordon ground and thermal states in the gapped phase [47].

Let us begin the discussion by noting that the phase is related to
fluctuations of the atomic gas. In its microscopic description, the atoms are
described by quantum fields with creation and annihilation operators denoted
as $\Psi^{\dagger}(z),\Psi(z)$, and the correlations of such atomic fields
should decay in space reflecting the locality of the interactions between the
atoms and the fact that finite temperature prevents long-range correlations
that would signify order in the system. The fluctuations of the phase field
are related to the correlations between the atoms such that the former
increase when the latter decrease. For example, in the quadratic harmonic
fluid approximation

$g_{1}(z,z^{\prime})=\langle\Psi^{\dagger}(z)\Psi(z^{\prime})\rangle\propto
e^{-\tfrac{1}{2}\langle(\hat{\varphi}(z)-\hat{\varphi}(z^{\prime}))^{2}\rangle},$
(19)

i.e., the phase difference field is related to the logarithm of the one-
particle density matrix of the atoms in the gas [47, 48]. The decay of off-
diagonal fluctuations of the atoms for increasing distance corresponds to
increasing correlations of the phase difference. In the gapless phase of the
gas, at finite temperature the relation of fluctuations to correlations is
implemented by exponential decrease of $g_{1}$, i.e., linear increase of phase
correlations, while in the ground state the scaling would be inverse algebraic
for $g_{1}$, i.e., logarithmic for phase correlations [48].

Figure 5: Scaling of phase auto-correlations $C^{\phi\phi}(r)$ as a function
of the distance $r$ from the reference point (middle) in the initial state.
Scan 1 corresponds to a massive Gaussian state, scans 2, 3 and 4 correspond to
non-Gaussian states at different interaction strength. The bands indicate the
size of the estimated errors.

In the gapped phase corresponding to the initial states of the experiment,
even though the quadratic approximation is not always valid, the relation
between the two correlation functions is qualitatively similar. Fig. 5
presents the auto-correlations

$\displaystyle
C^{\phi\phi}(r)=\langle(\hat{\varphi}(z_{0}+r)-\hat{\varphi}(z_{0}))^{2}\rangle$
(20)

measured in the experiment for initial state preparations with differing
barrier heights and hence different strength of the tunnel coupling $J$. We
find that over short distances the fluctuations increase seemingly
quadratically, eventually switching to a slower increase, corresponding to
linear growth or saturation. The rate of the increase depends on the strength
of the tunnel coupling.

The cross-over from quadratic to approximately linear scaling can be
interpreted by relating to atom fields as alluded to above. On this level we
find that the short-range correlation decay of $g_{1}$ is rapid and governed
by a Gaussian function while at large separations we find exponential decay of
correlations. The former may be non-universal, affected by the effective cut-
off of the field theory implemented by finite measurement resolution.
Additionally the short-range correlations may depend on _renormalization group
(RG)_ irrelevant terms in the Hamiltonian. On the other hand the long-range
scaling are expected to be robust and depend primarily on RG relevant terms.
This is what we find as the tunnel coupling is effectively described by the
sine-Gordon interaction (discussed in detail below) which is RG relevant.

In order to further interpret the scaling of the initial correlations, we can
think of the $C^{\phi\phi}$ function based on the local velocity field
correlations

$\displaystyle
C^{\phi\phi}(r)=\int_{z_{0}}^{z_{0}+r}\mathrm{d}{x_{1}}\int_{z_{0}}^{z_{0}+r}\mathrm{d}{x_{2}}\langle{\hat{u}}(x_{1}){\hat{u}}(x_{2})\rangle\
.$ (21)

This integral representation can be further simplified by the following
physical considerations. Given that the state is characterized by the presence
of massive excitations of the sine-Gordon model in combination with the fact
that the velocity field is local, its correlations are expected to decay
exponentially with some correlation length $\xi>0$. Under this assumption, as
$r$ becomes larger than $\xi$ we reach a linear scaling of $C^{\phi\phi}(r)$.
To see this, it is instructive to switch to the coordinates $x_{\pm}=x_{1}\pm
x_{2}$. The integral of the off-diagonal direction $x_{-}$ should give a
constant roughly proportional to the correlation length. The remaining
integral over the diagonal direction along $x_{+}$ should yield a scaling
proportional to $r$ because we are integrating a constant function. Together,
this results in a linear scaling

$\displaystyle C^{\phi\phi}(r)\approx\text{const.}\times r$ (22)

for $r\gg\xi$. Inspecting Fig. 5 we indeed find that for
$r>$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ a linear behavior is valid,
namely for all prepared initial conditions we see either a non-trivial linear
increase (scans 3 and 4 corresponding to relatively small $J$) or a leveling-
off (scan 1 corresponding to the largest $J$ and to a lesser extent scan 2 at
an intermediate value of $J$) compatible with a small or essentially
negligible slope constant. In the next paragraph, we will further elaborate on
this argument and explain how one can understand the saturation of the
correlations $C^{\phi\phi}(.)$ for large coupling using the phenomenology
derived from the Gaussian KG field theory.

####  3.1.1 Quantum field simulation of the relativistic field theory models

The state preparation in the experiment involves open system dynamics due to
cooling by evaporation or atom losses [49] which yields initial conditions
closely matching thermal theory [44, 25], consistent with predictions that can
be derived within _Tomonaga-Luttinger liquid (TLL)_ and _Klein-Gordon (KG)
model_ as special limits of the _sine-Gordon (SG) model_. We will now give
more details about this description which allows us to capture the system for
the limits of a strong and weak coupling of the adjacent one-dimensional gases
as depicted. For a high double well barrier, the coupling between the two
wells vanishes and the effective sine-Gordon description reduces to the
Gaussian TLL model. For a low double well barrier, on the other hand, we are
in the limit of large coupling and the description becomes again effectively
Gaussian. This can be seen in terms of a semi-classical description: the
system lies at the bottom of a very steep cosine potential which can therefore
be approximated by a parabola. In this case the effective description is given
by the KG theory (see the experimental study [22] and the numerical
theoretical analysis [37] for detailed discussions on the crossover from the
TLL to the KG regime of the SG model).

Following Refs. [12, 22, 33], we consider the effective field theory model
describing the fast decoupling of two adjacent one-dimensional gases of
neutral atoms at ultra-cold temperatures (see Fig. 2 for a schematic of the
experimental setup). At strong tunnel coupling between the two wells, the
effective model is given by the Hamiltonian

$\begin{split}\hat{H}_{\text{KG}}(J)=\int\mathrm{d}z\biggl{[}&\frac{\hbar^{2}n_{\mathrm{GP}}(z)}{4m}\left(\partial_{z}\hat{\varphi}(z)\right)^{2}+g(z)\delta\hat{\varrho}(z)^{2}+Jn_{\mathrm{GP}}(z)\hat{\varphi}(z)^{2}\biggr{]}\
,\end{split}$ (23)

involving relative fluctuations in phase $\hat{\varphi}(.)$ and density
$\delta\hat{\varrho}(.)$ [50, 30]. These low-energy degrees of freedom
represent the phononic excitations of a one-dimensional Bose gas and satisfy
bosonic commutation relations
$[\delta\hat{\varrho}(z),\hat{\varphi}(z^{\prime})]=\mathrm{i}\delta(z-z^{\prime})$.
Here $m$ is the atomic mass and $n_{\mathrm{GP}}(.)\approx\texttt{const}$ is
the mean density profile. The phononic operators are defined within the atomic
cloud whose spatial extension is given by the support of $n_{\mathrm{GP}}$.
Lastly, $g(.)$ is the density-broadened interaction strength [51, 30] and the
parameter $J$ is the tunnel coupling which is tuned by the double-well barrier
height.

The last term of the Hamiltonian that involves $J$ can be viewed as a mass
term. Ramping up the barrier height in the experiment, effectively quenches
the mass term from a nonzero value of $J$ to zero, so that only the first two
terms that make up the TLL Hamiltonian remain, i.e.,

$\displaystyle\hat{H}_{\text{KG}}(J=0)=\int\mathrm{d}z\biggl{[}$
$\displaystyle\frac{\hbar^{2}n_{\mathrm{GP}}(z)}{4m}\left(\partial_{z}\hat{\varphi}(z)\right)^{2}+g(z)\delta\hat{\varrho}(z)^{2}\biggl{]}=\hat{H}_{\text{TLL}}.$
(24)

Note that the above effective description in terms of a quantum quench from
the massive KG to the massless TLL model is valid only under the condition
that the density $n_{GP}$ is constant in both space and time. This holds to a
very good approximation in the experimental system when a box-shaped external
trap is used, and it is also a sufficiently good approximation in the middle
region of the system when a parabolic trap is used instead.

Such a KG to TLL quench is performed in experimental scan 1 and the
corresponding measurements of the dynamics will be shown below. Here we will
focus on the properties of the initial state preparation. The experiment tends
towards the KG regime for very strong tunnel couplings as evidenced here as
well. In the KG limit there is strong pinning of the compactified phase field
to the minimum of the cosine potential which results in the absence of phase
winding, that is, absence of soliton excitations. This means that the phase
difference between the edges of the system is very small

$\displaystyle\int\mathrm{d}{z}~{}\partial_{z}\hat{\varphi}(z,t=0)=0\ .$ (25)

Accordingly, in scan 1 for which the initial state is in the KG regime we
observe that the slope in (22) is close to zero and the large distance
asymptotics is saturated to a constant (Fig. 5), in contrast to other scans,
especially 3 and 4 corresponding to relatively small $J$, which are consistent
with an approximately linear increase with the distance.

Summarizing, the correlations of the phase should be thought of as
fluctuations as they govern the melting of the ordering of the atomic system
and tend to increase with the distance. On the other hand, the correlations of
the derivative field ${\hat{u}}$ in equilibrium states corresponding to the KG
limit should be decaying with the distance as a result of pinning of the phase
due to heavy phononic excitations. The scaling observed in the experiment
indeed reproduces the phenomenology that correlations of the derivative field
${\hat{u}}$ decay with the distance.

####  3.1.2 Velocity field correlations in the initial state

In this subsection, we are interested in the scaling, in particular, the range
of correlations of the local velocity field in the initial equilibrium states,
in order to assess the presence of one of the ‘pillars’ of the spatial
scrambling mechanism. Extracting the velocity field two point correlations

$\displaystyle
C^{uu}(z_{1},z_{2})=\langle\partial\hat{\varphi}(z_{1})\partial\hat{\varphi}(z_{2})\rangle$
(26)

directly is impossible in the experiment. Instead we look at the approximation
obtained by looking at the correlations of phase differences between adjacent
pixels

$\displaystyle C^{uu}(z_{1},z_{2})\approx\delta
z^{(-2)}\langle(\hat{\varphi}(z_{1})-\hat{\varphi}(z_{1}+\delta
z))(\hat{\varphi}(z_{2})-\hat{\varphi}(z_{2}+\delta z))\rangle.$ (27)

In practice, we have $\delta z\approx$2\text{\,}\mathrm{\SIUnitSymbolMicro
m}$$ which in typical configurations amounts to $2$ to $4\%$ of the system
length.

Figure 6: _Top_ : Longitudinally averaged velocity field correlations
$\overline{C^{uu}}(r)$ for state preparations used in Ref. [24]. From left to
right the tunneling strength is reduced. For all scans we see an auto-
correlation peak for $r\approx 0$. We see that in the KG regime (scan 1) there
is a rapid decay of the auto-correlation switching over to negative values,
i.e., an anti-correlation of the velocity field at a distance of about
$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$ which is about $20\%$ of the system
length. As $J$ decreases the anti-correlation dip becomes more shallow,
eventually vanishing almost entirely. For all scans the velocity field 2-point
functions are consistent with exponential fall-off envelopes with a length-
scale which is a small portion of the system size. _Bottom_ : Theoretical
plots of the correlation function of the velocity field $C^{uu}(r)$ as a
function of the distance $r$ in the Tomonaga-Luttinger liquid and Klein-Gordon
model (large coupling limit), at one choice of parameter values $c=1$, $M=1$,
$T=0$ and various values of the high-momentum cutoff $\Lambda$. We observe the
characteristic negative correlation dips on the two sides of the central peak.

Fig. 6 shows the extracted profiles where we average over all pairs of
positions $z_{1}$ and $z_{2}$ with a fixed distance $r=z_{1}-z_{2}$. For large
$J$, i.e., in the KG regime (scan 1), we find once again a profile consistent
with a vanishing integral in the distance $r$, as discussed above, implemented
by a strong anti-correlation on the two sides of the central peak that cancels
the auto-correlation at zero distance $r\approx 0$. The profile matches
qualitatively with the theoretical prediction for equilibrium correlations in
the KG model, also plotted in Fig. 6 for comparison with the experimental plot
for scan 1. The theoretical formula for the distance dependence of the
correlations in the KG model in a large homogeneous system is (see for example
[52])

$C^{uu}_{KG}(r)=-\partial_{r}^{2}C^{\phi\phi}_{KG}(r),\qquad
C^{\phi\phi}_{KG}(r)\propto\int_{-\infty}^{+\infty}\mathrm{d}k\;\frac{e^{ikr-\tfrac{1}{2}(k/\Lambda)^{2}}}{E(k)}\left(\frac{1}{2}+n_{\text{BE}}(E(k),T)\right),$
(28)

where $E(k)=\sqrt{c^{2}k^{2}+M^{2}c^{4}}$ and $n_{\text{BE}}(E,T)$ is the
Bose-Einstein distribution at temperature $T$. The parameter $M$ is the KG
particle mass, which is $M\sim\sqrt{J}$ here, and $\Lambda$ is a high-momentum
cutoff. At large distances the above correlation function decays exponentially
with a correlation length determined by the mass $M$ and temperature $T$,
while at short distances the theoretical scaling is the same as in the TLL
ground state, $M=0$, $T=0$, which is $\sim-1/r^{2}$. However, the contribution
of high momentum modes is suppressed by the cutoff $\Lambda$, which in the
experiment is dominated by the finite imaging resolution which has the form of
a Gaussian weight function [33]. As a result the profile of the correlation
function switches from a singular function to the shape shown in Fig. 6. The
precise height and width of this profile is controlled by all parameters $M,T$
and $\Lambda$. For this reason and given that a precise estimate of $M$ and
$T$ is not available, a quantitative comparison of the experimental and
theoretical plots is not possible, however, the qualitative behavior is
similar.

Decreasing $J$ to intermediate values in the sine-Gordon regime, the anti-
correlation dips are weakened and the integral over $r$ becomes nonzero in
agreement with the linear increase of ${C^{\phi\phi}}(r)$ at large distances
observed in the previous section. In the context of the sine-Gordon model the
physical meaning of the velocity field, i.e., the phase derivative, is the
density of solitons in the system. A non-vanishing value for the correlations
of the velocity field integrated over a spatial region signifies fluctuations
of the number of solitons in that region. This is precisely what one would
expect for equilibrium states of the sine-Gordon model in the strongly
correlated regime. Therefore the above scaling analysis provides an indirect
signature of the presence of solitons and of their thermodynamics in the
experiment, already observed in the full counting statistics of the phase
field in Ref. [22].

###  3.2 Canonical decoupling into Gaussian and non-Gaussian sectors

It is often the case that the internal dynamics of a physical system can be
very well described by a Hamiltonian consisting of two parts where
individually each part involves only a commuting set of operators and the non-
commutative, i. e., quantum, character of the model comes from the fact that
the two parts are non-commuting. This is very well illustrated by the phononic
degrees of freedom that we have in mind where typically, e.g., for equilibrium
conditions at low temperatures we have the generic form

$\displaystyle\hat{H}=\hat{H}_{\phi}+\hat{H}_{\rho}\ .$ (29)

Individually, the thermal state of only one part of the canonical pairs, so
$\hat{H}_{\phi}$ or $\hat{H}_{\rho}$ would agree with a classical probability
distribution, but in general for $\hat{H}$ this need not be true.
Interestingly, at sufficiently large temperatures the bosonic statistics of
the degrees of freedom is expected to become less prominent. Whenever this is
true, such an effect makes sure that the thermal state of the entire
Hamiltonian

$\displaystyle\sigma_{\beta}=e^{-\beta\hat{H}}/Z_{\beta},\,\,Z_{\beta}=\mathrm{tr}[e^{-\beta\hat{H}}]$
(30)

agrees closely with a classical probability distribution. To illustrate this
point, in the most extreme case of very high temperatures, we have

$\displaystyle\sigma_{\beta}\approx(1-\beta\hat{H})/Z_{\beta}\ .$ (31)

Here, we see that the correlation function of either the phase or density
operator depends on its respective part of the Hamiltonian and so the non-
commuting aspect of the fields does not even have a chance to play a role.

This observation coming from the simple high-temperature expansion may be
valid only at temperatures much higher than in the experiment. However, to
make the extrapolation to lower temperatures plausible, one should notice that
as the temperature increases, the energy of the system must also increase, so
we must have that $\langle{\hat{\phi}}_{k}^{2}\rangle$ and
$\langle{\delta\hat{\rho}}_{k}^{2}\rangle$ must also grow with temperature.
This means that the occupation numbers of the modes involved will become
increasingly larger than the vacuum level. However, only for states close to
the ground state the canonical commutation relations play a prominent role
because the expectation value of the commutator operator, which is independent
of the state and proportional to $\hbar$, is comparable to the mode occupation
numbers.

This leads to the _classical field approximation (CFA)_ being a good
approximation, where we assume that the phase and density operators are
effectively commuting with each other. This suggests that we can effectively
‘trace-out’ one of the sectors of the operators without quantum corrections,
so by treating the phononic fields as independent classical fields not coupled
with each other. For the lack of a better term, we will refer to such an
independence of the correlation functions of fields in each canonical sector
from the canonically conjugate sector as canonical decoupling. Whenever such a
feature is at play we can still make use of further observations to plausibly
infer characteristics of the correlations of the quantum many-body system
considered.

Specifically, let us assume that one of the canonical parts of the Hamiltonian
is a quadratic form of the degrees of freedom. Anticipating the relation to
the experiment, let us take the density part to be quadratic while the phase
part can include a many-body potential going beyond the quadratic term, as in
the sine-Gordon model. Under the assumption that the two canonically conjugate
fields are effectively decoupled for the temperature in question, we obtain a
very crucial prediction, namely that the density-density fluctuations can be
approximated as

$\displaystyle\langle\delta\hat{\varrho}(z_{1})\ldots\delta\hat{\varrho}(z_{n})\rangle\approx\mathrm{tr}\left[\delta\hat{\varrho}(z_{1})\ldots\delta\hat{\varrho}(z_{n})e^{-\beta\hat{H}_{\rho}}\right]/Z_{\beta}^{(\rho)}$
(32)

where again $Z_{\beta}^{(\rho)}=\mathrm{tr}[e^{-\beta\hat{H}_{\rho}}]$ for
$\beta>0$. This correlation can then be obtained using Wick expansion of the
only non-trivial correlation function, i. e., the second moments

$\displaystyle\mathcal{Q}(z,z)=\langle\delta\hat{\varrho}(z)\delta\hat{\varrho}(z^{\prime})\rangle.$
(33)

This expression makes specific what we mean by tracing out: The correlations
in the density sector are taken to be described by the thermal state of the
density Hamiltonian $\hat{H}_{\rho}$. We take the same temperature as we would
take for the full thermal state of the full Hamiltonian. The former
observation should be stressed as the independence of tunnel coupling means
that at all its values $\hat{H}_{\rho}$ is Gaussian so actually the
correlations in the density sectors should be Gaussian too. This means that
the density fluctuations should be _independent_ of the tunnel coupling and
approximately thermal. In summary, the argument suggests that, crucially, CFA
implies that at sufficiently high temperatures _all_ higher-order moments in
the density sector should be _Gaussian_. Thus the Gaussian bath in the
canonical transmutation mechanism is a result of canonical decoupling and at
least one of the canonical sectors being Gaussian.

By allowing for tunneling of atoms between the adjacent gases an effective
interaction between the phonons becomes relevant which can give rise to kink
excitations according to the _sine-Gordon_ (SG) model whose Hamiltonian is
given by

$\displaystyle\hat{H}_{\text{sG}}=\hat{H}+J\int_{-L/2}^{L/2}\mathrm{d}z\
n_{\mathrm{GP}}(z)\cos(\hat{\varphi}(z))\ .$ (34)

Here $J/(2\pi\hbar)$ is the single particle tunneling rate, which can be tuned
by the barrier height.

Having this specific model in mind, lets us state what the implications of
canonical decoupling would be in our system. First, the factorization of the
thermal state density matrix at high temperature means that the density
correlations would be given by expectation values computed in a state
$\sim\mathrm{e}^{-\beta\hat{H}_{\rho}}$ which is a diagonal quadratic form of
the density field. This means that all higher-order functions of the density
would be Wick contractions based on the two-point function which at high-
temperatures reads

$\displaystyle\mathcal{Q}(z,z^{\prime})\approx(\beta
g)^{-1}\delta(z-z^{\prime})\ .$ (35)

The value of phase correlation functions can in turn be computed as explained
in Ref. [27] by replacing the quantum phase operators by classical phase
variables since the canonically conjugate operators should not play a role at
high temperatures, so only commuting operators would be involved. Put
differently, we would take the phase correlation functions to be the
correlation functions in the classical sine-Gordon model.

In the large temperature limit the correlation length of the density
fluctuations is very small. This in particular should imply that to a good
approximation we can write for eigen-mode operators

$\langle\delta\hat{\varrho}_{k}^{2}\rangle=\text{Const}.$ (36)

We can recover this qualitative feature by starting from the formula for
thermal correlations of a set of eigen-modes

$\displaystyle\left\langle\delta\rho(x_{1})\delta\rho(x_{2})\right\rangle_{0}$
$\displaystyle\propto\sum_{n=1}^{\infty}\left[1+2n_{\text{BE}}(E;\beta)\right]E({k_{n}})\mathrm{e}^{\mathrm{i}k_{n}(x_{1}-x_{2})}$
(37)

where $n_{\text{BE}}(E;\beta)=({\mathrm{e}^{\beta E}-1})^{-1}$ is the Bose-
Einstein distribution at inverse temperature $\beta>0$. In the above the
eigen-mode energy is denoted by $E({k_{n}})$, while $c$ and $K$ are the sound
velocity and Luttinger parameter, respectively (see Appendix A for details).
We now see that for $\beta^{-1}\to\infty$ we have
$n_{\text{BE}}(E;\beta)\approx 1/(\beta E)$ from which we find

$\displaystyle\left\langle\delta\rho(x_{1})\delta\rho(x_{2})\right\rangle$
$\displaystyle\approx\text{Const}.\times\delta(x_{1}-x_{2}).$ (38)

By inspecting formula (37), we see that in the Gaussian case we will see
quantum corrections, related to vacuum fluctuations, only once the first order
expansion of the exponential in the Bose-Einstein distribution becomes an
inadequate approximation. This sets a scale for the temperature for Gaussian
CFA to be closely linked to the eigen-mode frequency of the low-lying modes
observable in the experiment. To account for such short-range correlations we
would indeed see that

$\langle\delta\hat{\varrho}_{k}^{2}\rangle\approx\text{Const}.$ (39)

Let us remark that in the case of an interacting system, an approximate
Gaussian description at large temperature would include effects of the
interaction through the self-energy corrections [52]. This asymptotic behavior
is generally valid independently of the precise form of $E({k})$ and therefore
independent from the interaction which enters only through the self-energy
insertions in $E({k})$. In the present case, the above result justifies (36)
and explains why that constant is independent of the sine-Gordon tunnelling
$J$.

Summarizing the entire section, here we have argued that phase and density
decouple at high temperatures if there are no terms in the Hamiltonian
directly involving both types of operators. In particular, this would imply
that density fluctuations should be Gaussian in the SG model which seems to be
the case for temperatures relevant in the experiment. This is the canonical
decoupling pillar summarized in Fig. 1. Of course, this property alone does
not imply that the Gaussian density sector dominates the phase sector during
the dynamics. The next section explains that this is actually expected to be
the case for the experimental system.

###  3.3 Dominance of density over phase initial correlations in the
experimental states

####  3.3.1 Intuitive understanding of the experimental findings

We are now turning to discussing the third ‘pillar’ of Gaussification by
canonical transmutation, which is the condition where the initial correlations
of one of the two canonical sectors dominate over those of the other. The
correlations of each sector at intermediate times $t>0$ that are not integer
multiples of a half recurrence period will be a mixture of the initial
correlations of both sectors. Therefore, if one of the sectors has larger
correlations than the other initially, it will dominate in the time evolved
correlations of both sectors.

Figure 7: Semi-classical explanation of the dominance of density over phase
fluctuations in the sine-Gordon model. Ignoring the spatial variation of the
phase, the sine-Gordon model can be approximated as the quantum analogue of
the classical pendulum [53]. _Top:_ At a fixed temperature (corresponding to
mean total energy denoted by the dashed black line), increasing $J$ results in
a steeper and narrower potential well that can be approximated by a parabola
(dashed colored lines). _Bottom:_ Seeing from a phase-space point of view,
increasing $J$ results in stronger ‘squeezing’ of the wave-function,
illustrated through the elongation of the dashed line contours that represent
the energy levels (or the shape of a phase space distribution of a classical
thermal ensemble). The values of $J$ from purple to red are 0.01, 0.025, 0.05,
0.1, 0.25, 0.5, 1., 2.5 and the mean total energy level is $E=0.1$.

Let us now explain why this can well be the case in the experiment. As already
mentioned, the initial states are effectively equilibrium states of the SG
model or, in case of large coupling, of the KG model. For clarity, consider
the phase and density fluctuations in thermal states of the two limiting cases
of the SG model, i.e., the KG and TLL models, (23) and (24) respectively,
which are quadratic. In the TLL model expressed as a sum of decoupled modes
(14), we know from the equipartition theorem that the energy contributions of
each of the two quadratures to the total energy are equal

$\displaystyle\langle\delta\hat{\varrho}_{k}^{2}\rangle=\omega_{k}^{2}\langle{\hat{\phi}}_{k}^{2}\rangle\qquad\text{(TLL
equilibrium states)}.$ (40)

Similarly, in the KG model equipartition means that

$\displaystyle\langle\delta\hat{\varrho}_{k}^{2}\rangle=\omega_{0k}^{2}\langle{\hat{\phi}}_{k}^{2}\rangle\qquad\text{(KG
equilibrium states)}$ (41)

where $\omega_{k}=ck,\omega_{0k}=\sqrt{c^{2}k^{2}+M^{2}c^{4}}$ are the mode
frequencies of the TLL and KG model respectively with $M$ the effective mass
of the KG excitations, which increases with the coupling as $\sqrt{J}$. From
this relation, we find that because the coupling $J$ is large, the density
fluctuations $\langle\delta\hat{\varrho}_{k}^{2}\rangle$ are much larger than
the phase fluctuations $\langle{\hat{\phi}}_{k}^{2}\rangle$. In physical
terms, this is due to the energetic penalty imposed by the steep parabolic KG
potential on the phase fluctuations. The same argument applies on equilibrium
states of the SG model for any nonzero value of the coupling $J$ since the
energy is related to the correlations of the density and phase field and there
is an energetic penalty on phase fluctuations only. The underlying semi-
classical argument is illustrated in Fig. 7. More generally, for any model in
which there is an energetic penalty, induced by a potential, which applies on
_only_ one of the two canonically conjugate fields (here, the phase), the
corresponding fluctuations will be suppressed compared to those of the other.

In the experiment, there is an initial penalty on phase fluctuations due to
the tunnel-coupling $J$, while the density fluctuations are free to fluctuate
according to the given temperature. Hence, their amount will be larger than
the fluctuations in the phase sector. By this argument, we find that the phase
correlation function (initially non-Gaussian) should increase in magnitude
after the quench in the form of an increase of the Gaussian component.

####  3.3.2 Estimation of density fluctuations via phonon tomography

Figure 8: _Top:_ Reconstructed initial state correlations in the eigen-mode
basis. _Bottom:_ The tomographic recovery is performed by matching the
measured phase second moments at times
$t={0,1,\ldots,}$12\text{\,}\mathrm{ms}$$. Here, we show the comparison for
the initial time. The cost function $\Delta\Phi$ is the relative deviation of
the reconstructed second moments $\Phi_{\text{Rec}}$ from the data
$\Phi_{\text{Data}}$, weighted by the estimate of the root-variance at the
given point.

Figure 9: The tomographically reconstructed diagonal quadrature moments
confirm the dominance of density over phase initial correlations in the
experiment: The plots show the correlations
$V^{\phi\phi}_{k,k}=\omega_{k}\langle{\hat{\phi}}_{k}{\hat{\phi}}_{k}\rangle$,
$V^{\rho\rho}_{k,k}=\omega_{k}^{-1}\langle{\delta\hat{\rho}}_{k}{\delta\hat{\rho}}_{k}\rangle$
and
$V^{\phi\rho}_{k,k}=\tfrac{1}{2}\langle{\hat{\phi}}_{k}{\delta\hat{\rho}}_{k}\rangle$
as functions of the mode number $k$ in the initial state of each of the four
scans, as reconstructed by means of the tomographic recovery of second moments
of quadratures. We observe that in all scans,
$V^{\rho\rho}_{k,k}>V^{\phi\phi}_{k,k}$ at least by a factor of 2 in the
lowest modes, and also
$V^{\phi\rho}_{k,k}<V^{\phi\phi}_{k,k},V^{\rho\rho}_{k,k}$ by even an order of
magnitude in most cases.

In the experimental setting considered, only one quantum field quadrature is
directly experimentally accessible (phase), but not the canonically conjugate
quadrature (density). Using the quantum read-out method of Ref. [25], however,
we can estimate the content of density fluctuations as suggested above. In
this procedure, we reconstruct the full covariance matrix of the initial
unknown state by relating the non-equilibrium second moments of the phase
$\Phi(z,z^{\prime},t)$ at points $z,z^{\prime}$ and for the time $t\geq 0$
through the known TLL evolution equations, collecting all suitable second
moments. The results of this most simple theoretical model agree very well
with the experimental observations.

In our case, the phase correlations $\Phi(z,z^{\prime},t)$ defined in Eq. (2)
can be measured at a discrete set of points $z,z^{\prime}$ (with a spacing
given by the pixel size of the camera $\delta
z\approx$2\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ [33]) at time various
$t\geq 0$. Apart from the pixel size, other effects, including diffraction,
limit the spatial resolution. The measured values can be related to
theoretical continuum predictions by implementing a real-space cut-off via a
Gaussian convolution with standard-deviation
$\sigma\approx$3.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$$. From Eq. (15), we
see that for a sufficiently large number of pixels and time snapshots, the
eigen-mode correlations of the phase can be extracted first and from them the
corresponding eigen-mode correlations of the density through fitting the
dynamics to those predicted by phase space rotation, making a quantitative
reconstruction feasible. The implementation of the tomographic reconstruction
[25] extends and optimizes this intuitive idea using _convex optimization
techniques_ to find the full covariance matrix (including density
correlations) $\Gamma$ at the initial time $t=0$ such that the corresponding
forwards propagation $\Gamma(t)$ exhibits phase correlations matching the
observed data, under the constraints arising from the Heisenberg uncertainty
principle. This approach makes a lot of sense in the light of the fact that
the Heisenberg uncertainty principle for covariance matrices sets a lower
bound on quadrature fluctuations, giving rise to a semi-definite and hence
convex constraint. In Ref. [25], the functioning of this method has been laid
out in detail and convincingly demonstrated in a study of effectively Gaussian
state dynamics with large initial tunneling and quenching to zero tunneling
[30]. Here, we apply it to derive the second moments of the (generally, non-
Gaussian) quantum initial states corresponding to the present quench
parameters.

Fig. 8 shows the application of the tomographic recovery of second moments. We
find that the eigen-mode quadratures, defined through the rescaled phase and
density mode variables

$\displaystyle\tilde{{\hat{\phi}}}_{k}={{\hat{\phi}}}_{k}\sqrt{\omega_{k}}\quad,\qquad\tilde{{\delta\hat{\rho}}}_{k}={{\delta\hat{\rho}}}_{k}/\sqrt{\omega_{k}},$
(42)

are squeezed in the sense that the diagonal correlations in the density sector
$V^{\rho\rho}_{k,k}$ are larger than those in the phase sector
$V^{\phi\phi}_{k,k}$. This statement pertains directly only to second moments.
The above is further verified at quantitative level in Fig. 9 where the
diagonals of the eigen-mode quadrature correlations are plotted for each of
the four main scans of the experiment. In all scans, $V^{\rho\rho}_{k,k}$ is
larger than $V^{\phi\phi}_{k,k}$ for almost all modes, at least by a factor of
2 in the lowest modes which are the ones that weigh in most in the calculation
of spatial averages of correlations. Also note that $V^{\phi\rho}_{k,k}$ is
smaller than $V^{\rho\rho}_{k,k}$ and $V^{\phi\phi}_{k,k}$ for almost all
modes, as expected for a rapid quench. This is because

$V^{\phi\rho}_{k,k}\propto-\frac{\rm{d}}{\rm{d}t}\langle{\hat{\phi}}_{k}{\hat{\phi}}_{k}\rangle|_{t=0},$
(43)

which vanishes for instantaneous quenches since there is no dynamics prior to
$t=0$. For quenches of a nonzero ramp duration $\Delta t$ that is short
compared to the characteristic time scale of the dynamics, i.e., $\Delta t\ll
1/\omega_{k}$, the corresponding $V^{\phi\rho}_{k,k}$ is nonzero, given that
the time evolution starts already before $t=0$ at the beginning of the quench
ramp, but still small compared to $V^{\phi\phi}_{k,k},V^{\rho\rho}_{k,k}$. The
negative sign of $V^{\phi\rho}_{k,k}$ in most of the scans is due to the fact
that the phase correlations increase as a result of the mixing with the
density correlations.

These observations verify quantitatively the accuracy of the fitted model for
the dynamics. Moreover, they partially verify the validity of the condition of
dominance of density over phase fluctuations in the initial state, as required
for the canonical transmutation mechanism. We find that the former are indeed
larger than the latter, even though not sufficiently larger to claim that they
are dominant. We rather find that their ratio is of order $\sim 2$.
Nevertheless, as we will later see in Fig. 15, this magnitude is in agreement
with estimates based on the classical field approximation and even though only
moderately large it is sufficient to explain the observed decay in the
experiment. Importantly, the tomographic reconstruction provides an
independent verification that the canonical sectors are sufficiently squeezed
as required for the canonical transmutation mechanism, based on directly the
experimental data rather than on the classical field approximation.

Let us remark that in principle the tomographic recovery method employed here
can be generalized to higher order correlations, but there are two
obstructions that make this step largely impractical. Firstly, if the
reconstructed correlations are corrupted by noise, they would not necessarily
precisely reflect those of a physical quantum state, as is well-known. The
constraints that stem from quantum uncertainty bounds on the correlations can
only be easily implemented on the level of second moments as a _semi-definite_
, that is to say a convex, constraint. This alone does not ensure, however,
that the full quantum state is a positive semi-definite operator. To ask
whether such a consistent quantum state can be found is related to the quantum
marginal problem, a computational problem that has no efficient solution. This
is an obstruction commonly faced when generalizing Gaussian tomography to full
quantum state tomography [54]. The second issue is that the available number
of experimental runs is limited, in other words, there is a limitation on the
maximum size of the statistical sample. As discussed and exemplified in the
supplementary material of Ref. [25], even just the second moment
reconstruction would profit from larger sample sizes as this reduces the
statistical errors on the input second moments and hence improves the
precision of the reconstruction. Increasing the number of experimental runs is
challenging, however, as it is hard to ensure the stability of the calibration
of the experimental setup over extended periods of time. This leads to a
significant obstruction for reconstructing higher-order correlation functions:
The amount of data to be reconstructed grows exponentially with the order but
the sample size and hence precision of the input is constant in practice. This
is the reason why we have not attempted the reconstruction of, e.g., the
4-point correlation tensors to assess the question whether the density sector
has negligible cumulants.

Fig. 10 displays the reconstructed full covariance matrix to show the second
moments of density fluctuations in real space. As shown in Fig. 10, we find
that the reconstructed density correlations match closely the density
fluctuations obtained from a Gaussian thermal state in the KG model. The
thermal theory deviates from the exact computation in continuum (35) due to
the application of an appropriate smearing to account for the finite
measurement resolution which affects the input to the tomography. More
specifically, the tomography is performed over a constant number of the
lowest-lying modes, which effectively implements a hard cut-off on the higher
energy modes. The eigen-mode wave functions are oscillatory, even when
computing them numerically in the presence of inhomogeneities. As is known
from quantum field theory – and more generally from harmonic analysis – a
Fourier transform of a discontinuous signal has a long-range support. We see
this effect also in our case: If we set the occupation of the KG thermal modes
to zero above the chosen cut-off, we obtain long-range off-diagonal artifacts.
This would not be the case for a smooth cut-off, say an exponential or a
Gaussian decay. Nonetheless, we see that the local information, centered
around the diagonal, is already converged in shape and magnitude signifying
the matching in real space. In the appendix, we present the results for a
reconstruction over a doubled-up number of available eigen-modes and we find
that the cut-off effect is reduced, with the concentration around the diagonal
being increased.

Through this, we have verified the presence of the density fluctuations in the
initial state of the system as suggested above. The most significant finding
is the verification based on quench data of the presence of a large amount of
density fluctuations in the system. This is depicted in Fig. 8 and 9, where we
encounter the situation that the values of $V^{\rho\rho}$ are significantly
larger than $V^{\phi\phi}$, similar to the results obtained earlier in Ref.
[25]. We then go beyond previous works and show that the off-diagonal
correlations in the reconstructed eigen-mode covariance matrix do not lead to
significant deviations from a diagonal real-space covariance matrix of the
diagonal correlations expected from thermal theory. This is shown in Fig. 10,
where by comparison to a theoretical computation with a hard cut-off we see a
very close matching.

Additional plots on tomographic reconstructions are presented in Appendix D.

Figure 10: Comparison of the reconstructed initial second moments of the
density with the theoretical predictions for a thermal state (at estimated
parameter values $J=$0.5\text{\,}\mathrm{Hz}$$ and
$T=$60\text{\,}\mathrm{nK}$$). The density plots show the correlations in
coordinate space in the entire system (length
$L=$60\text{\,}\mathrm{\SIUnitSymbolMicro m}$$). The reconstructed
correlations agree quite well with the theoretical ones. The first two plots
show an approximation of the continuum limit correlations corresponding to a
very fine discretization and computing the respective eigen-modes of the TLL
model, while the second two show a discrete grid of points corresponding to
the experimentally measured data, i.e. a coarse-grained representation of the
correlations with pixel size $\sim$2\text{\,}\mathrm{\SIUnitSymbolMicro m}$$.
We find that the transformation of the reconstructed eigen-mode correlations
to real space yields a reliable comparison around the diagonal. The off-
diagonal features can be directly linked to the presence of a maximum mode
cut-off and hence do not represent true correlations in the system.

##  4 Characteristics of the dynamics in the experiment

While we have so far focused on discussing properties of the initial quantum
states, we now turn to elaborating on the actual dynamics in the experiment.
The ultra-cold gas in the experiment can be thought of as a fluid which
responds to local external disturbances forming waves. These wave-packets are
phonons and their effective models, TLL and SG, have been discussed above. The
experiment reproduces the phenomenology of the equilibrium conditions captured
by these models making it an excellent candidate for a quantum simulation
platform of continuous quantum fields. Out of equilibrium the situation can be
markedly different, however, because the universality of thermal equilibrium
properties rooted in the renormalization group theory can be distorted by
previously irrelevant terms that may affect significantly the dynamics. This
section will discuss various qualitative aspects that are important features
for developing a complete phenomenological understanding of the non-
equilibrium dynamics in the system.

###  4.1 Signatures of light-cone dynamics

Figure 11: Phase difference correlations
$\overline{C^{\Delta\phi\Delta\phi}}(z_{1},z_{2},t)$ averaged over the mean
position $z=(z_{1}+z_{2})/2$ and plotted as a function of the distance
$r=|z_{1}-z_{2}|$ and time $t$. The presence of a light-cone shape which
separates the outside region where correlations are largely flat as initially
from the inside region where they almost instantly switch to their equilibrium
values is clearly visible (see also Ref. [41] for a detailed discussion). The
dashed green lines correspond to the theoretical estimate for the sound
velocity as captured by Eq. (58). The averaging, which is done over all pairs
of points at the same distance $r$ in the middle region of the system, i.e.,
along lines parallel to the diagonal of the phase covariance matrix, serves to
reduce the statistical fluctuations. At a given time, the middle region agrees
with the dynamics in a thermodynamically large system. After time
$t\sim$15\text{\,}\mathrm{ms}$$, phase correlations in the middle region of
the system start being causally connected to the edges of the system and the
dynamics is affected by finite size and boundary effects.

_Light-cone dynamics_ constitutes a crucial feature of the phenomenology of
phononic dynamics. The effective Hamiltonians often possess a relativistic
form giving rise to an effective Lorentz symmetry parametrized by the speed of
sound. Whenever the dynamics is governed by a gapless Hamiltonian which is to
a good approximation translation invariant and quadratic, then the low energy
spectrum is linear, i.e., wave-packets propagate practically without
dispersion. This is the characteristic property of TLLs and it is responsible
for the coherent propagation of particle and hole excitations at low energies
in the interacting one-dimensional Bose gas. This dispersion-less propagation
of wave-packets is generally sensitive to deviations from homogeneity,
phononic interactions and other perturbations. Nevertheless, the effects of
such perturbations in the experiment do not become important before a
relatively long time passes, which means that the dynamics is practically
dispersion-less for the entire duration of the experiment. This is what allows
us to observe a clear light-cone wave-front in the propagation of phase
correlations.

Light-cone dynamics have been first studied experimentally on the atom chip in
Ref. [41] on the level of the phase field. The propagation velocity (_speed of
sound_) can also be extracted by studying the appearance of recurrences for
different system sizes [30], the very appearance of which has been enabled by
making the trap homogeneous [55]. This can be even further improved by making
use of the innovation in the atom control brought about by the implementation
of the _digital-mirror device_ [56, 57] that allows for programmable
potentials in space and time.

A density plot demonstrating the light-cone propagation of phase correlations
in the experimental system is shown in Fig. 11, which presents the dynamics
that ensues following an interaction quench from a very strongly tunnel-
coupled state of two adjacent gases to an effectively decoupled system of two
independently evolving systems. As discussed above, the initial Hamiltonian
and the one after the quench involve two extreme regimes of tunnel-coupling
values and for this reason both systems are effectively Gaussian [22]. The
parameters of this experimental scan are similar to those of scan 1, with the
difference that the box trap size is larger
($100\text{\,}\mathrm{\SIUnitSymbolMicro m}$) and the time snapshots are
equally spaced at steps of $2\text{\,}\mathrm{ms}$. The observation of the
shape of the correlations wave-fronts and estimation of their width as a
function of time is rather challenging due to the non-trivial profile of the
initial correlations of both the phase and density fluctuations, which are
mixed together and parallelly transported through the system. Here, the fact
that the initial correlations of the local derivative fields are of short
range especially in the limit of large initial coupling $J$ allows us to
observe the shape of the light-cone edges which follow quite straight lines as
expected for a homogeneous system. Even though it is hard to tell if the wave-
fronts spread with time or not, which is what is crucial for the spatial
mixing mechanism, we can verify that their spreading is negligible at least
for times up to half a recurrence period. This is exactly the time window that
is relevant for the equilibration stage of the dynamics.

The data presented in the plots carry two principal sources of distortion. The
first is due to statistical uncertainties arising from finite sample
statistics. Secondly, the measured phase profiles are affected by limitations
in the spatial resolution of the experimental read-out [33]. The comparison to
a theoretical thermal state becomes possible after applying a Gaussian
convolution with variance $\sigma=$3.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$$
(this should be compared with the system size which is
$L\approx$70\text{\,}\mathrm{\SIUnitSymbolMicro m}$$). The observation of
light-cone propagation is an indication that the linear phonon dispersion, an
otherwise idealized property, is actually an accurate approximation of the
dynamics of the experimental system.

###  4.2 Recurrences of the initial state

Another characteristic of the dispersion-less nature of the TLL dynamics is
the emergence of recurrences. This is because recurrences fade easily when
energy levels deviate from being commensurate. Recurrences in an isolated
quantum many-body system have first been observed in Ref. [30] on the level of
the coherence phase correlations. Moreover, in Ref. [24] (and as shown in Fig.
3 based on the results of that experiment) it has been demonstrated that the
same recurrences can be observed even on the level of higher order
correlations, as a more stringent test of returning back to the initial state.

Let us try to better understand why recurrences are linked to the linear
dispersion of the dynamical model. First, recall that the energy levels of a
TLL confined in a box trap with hard wall boundary conditions (of either
Neumann or Dirichlet type) are integer multiplies of a fundamental frequency
$\hbar\omega_{0}=\hbar c\pi/L$ (see Appendix A), therefore the initial state
is expected to be fully recovered after time equal to $2L/c$. Physically, this
is the time needed for the lowest-energy phononic excitation to travel around
the entire length of the system twice and return back to its initial position
in the same direction of motion. At the same time all higher energy
excitations travel an integer number of times around the system. At half of
this time the lowest excitations have travelled once around the system and the
state is a mirror reflection of the initial one, which also results in a
recurrence of phase coherence, since given the reflection symmetry of the
initial state phase correlations are identical to the initial ones. Therefore
the first recurrence is expected at time $T_{rec}=L/c$. In Ref. [30], not only
one but two recurrences have been clearly observed, yet exhibiting a
decreasing amplitude indicating the eventual loss of coherence in the
experiment after a sufficiently long time.

The emergence of recurrences is a property very sensitive to the
commensurability of the energy spectrum. More than that, in a TLL the energy
levels are equally spaced, a characteristic consequence of the linear
dispersion, which makes the observation of recurrences easier. If the above
requirement is not met, approximate recurrences are expected, however, the
fidelity between the recurrence state and the initial one would be strongly
suppressed in the experimental settings. For this reason, it is not by
accident that the observation of recurrences was only achieved when it became
possible to construct external trapping potentials of a hard-wall box shape
[30]. Satisfying this condition required a new technique for the
implementation of a blue detuned optical dipole potential. Until then the
trapping potential used in the experiment was of the standard parabolic shape,
resulting in an incommensurate energy spectrum $E_{n}=\hbar c\sqrt{n(n+1)}/R$
where $R$ is the Thomas–Fermi radius of the trapped gas. Even though
approximate recurrences of coherence are expected for this type of spectrum
[58], these are strongly suppressed and so never observed in the experiment.
From the above, it should be now clear how fragile recurrences are as an
aspect of the dynamics and, for that reason, how accurate witness of the
linearity of the phononic dispersion relation their experimental observation
is.

###  4.3 Deviations from the non-interacting phonon dynamics in the experiment

The effective field theory model describing the dynamics in the system on the
level of non-interacting models constitutes a good approximation, but is
nonetheless an approximation after all. Here, we will discuss it in general,
naming the possible deviations that may be relevant in the experiment. The
model at hand can actually be derived by a low energy approximation of the
_Lieb-Liniger model_ by looking at the fluctuations around the classical
density and phase profiles obtained from the _Gross-Pitaevskii equation_.
Hence, one obtains a hierarchy of effective models each containing terms
including only up to a fixed number of phase and density fluctuation
operators. The first non-trivial order is quadratic and yields an effective
description of the system in terms of non-interacting phonons. As we have
shown, this first order approximation matches quantitatively the experimental
observations up to intermediate time scales, but for long evolution times
various deviations start becoming visible in the experiment. One aspect of the
dynamics that is generally sensitive to deviations is the recurrences that
occur because the system is isolated and the phonons have a sufficiently
linear spectrum [30]. The occurrence of recurrences quantitatively matches
with the predictions of the quadratic phononic model, in particular in terms
of the timing (recurrence period) which is dictated by the speed of sound and
length of the system. The former can be controlled by changing the mean atom
density and the latter by choosing the width of the box-like potential. For
the experimentally accessible values the quadratic model accounts for the
observed recurrence times. Furthermore, recurrences are not visible for non-
homogeneous systems which arise, e.g., from harmonic trapping, again in
accordance with this model [55].

While the overall timing of the recurrences seems to be robust, the quality of
the recurrences deteriorates over time [30]. This damping may be accounted for
by fluctuations in the atom number in each experimental realization, which
result in fluctuations of the speed of sound, by deviations from homogeneity
of the experimental mean density profile, or by higher-order interaction
effects present in a stochastic Gross-Pitaevskii model [30]. On the
experimental level, a tomographic analysis has shown that the damping of the
revivals is closely connected with a steady and irreversible reduction of the
thermal squeezing in the eigen-mode quadratures. During this squeezing damping
process there are only small indications of the overall growth of phonon
occupation numbers. This suggests that the higher-order terms responsible for
this effect are irrelevant under the renormalization group as the equilibrium
features remain intact, while out-of-equilibrium ones are not. A small drift
can also arise from a linear coupling between the relative and the symmetric
phase-density sectors of the two coupled condensates, which, however, is
expected to be less substantial than the non-linearity effects.

There are a number of additional effects that seem to play a negligible role
and are discussed in detail in Appendix C. These include in particular the
possibility of wave-packet spreading due to phononic dispersion non-
linearities. One source of such effects which is significant at large
wavelengths occurs when the trap is not homogeneous. For harmonic traps, in
particular, the energy dispersion is $\omega_{n}\propto\sqrt{n(n+1)}$, i.e.,
the spectrum is non-linear for the lowest eigen-modes, which correspond to
large wave-lengths. As the spectrum converges for large $n$ towards a linear
function, like in the homogeneous trap case, sufficiently compact wave-packets
are unaffected by this non-linearity and do not disperse substantially enough
to play a key role in the Gaussification dynamics. Hence, we conclude that
even in the case of dynamics in a parabolic trap the effect is not strong
enough in the experiment to facilitate the delocalizing dynamics pillar of the
Gaussification by spatial scrambling mechanism.

The second possibility for a non-linearity that could potentially lead to
dispersion of wave-packets comes from the finite healing length of the atomic
gas. The effect of this on the spectrum can be accounted for by the Bogoliubov
dispersion relation, which gradually deviates from linear at small momenta to
parabolic at high momenta. This deviation enters through the Galilean kinetic
energy part of the Bose gas energy. However, as discussed in Appendix C, the
$k^{3}$ correction to the linear dispersion relation $\omega_{k}=ck$ is
negligible compared to the momentum cut-off dictated by the read-out spatial
resolution and its spreading effect is also negligible for the duration of the
experiment. Again, while cubic perturbations to the spectrum could potentially
induce a dynamical delocalization effect as needed for the spatial scrambling
mechanism, we find that this does not play any role in the experiment.

Finally, the effective boundary conditions at the edge of the system can
affect the dynamics. While Neumann conditions seem to be the physically most
soundly motivated, since they correspond to a vanishing particle current at
the edges, the fact that the density profile falls off smoothly rather than
sharply means that the assignment of an effective length for the experimental
system cannot be unambiguously done at arbitrary precision. For this reason,
it is reasonable that the right effective boundary condition of the system is
of mixed type with a small admixture of the Dirichlet type. However, for all
practical purposes Neumann boundary conditions seem to allow accounting for
all observations made.

##  5 Characteristics of Gaussification in the experiment

Before we analyze the above experimental observations into the context of the
earlier presented Gaussification mechanisms, we discuss some further
characteristics of Gaussification in the experiment and different ways of how
one can observe and study it.

###  5.1 Gaussification in the dynamics of full counting statistics

The picture of Gaussification put forth here complements the previous
discussion, in that Gaussification is now being considered from the
perspective of full counting statistics and its dynamics. Each interferometric
measurement in the experiment yields a phase profile which consists of a few
dozen data points for typical system sizes. Treating the measured phase
profiles as samples of a random variable (classical, because the quantum phase
operators at equal times are mutually commuting for any two positions) we can
study the time evolution of its distribution by plotting the corresponding
histograms at various times.

Figure 12: Histograms of the phase difference distributions and their dynamics
for each of the four main experimental scans. _Top_ : Histograms of the state
just before the quench ramp ($t=-2$ ms). Scan 1 corresponds to a narrow bell-
shaped distribution that is close to the Gaussian shape. Scans 2, 3 and 4,
while still characterized by bell-shaped distributions in the center,
progressively deviate from the Gaussian for larger $\varphi$, exhibiting long
tails, as one would expect from thermal states of an interacting model with a
cosine potential. _Bottom_ : Histograms at different times in logarithmic
scale: just before the quench ramp, at an intermediate time and at the last
measured time. The plots are presented in logarithmic scale to facilitate the
identification of the Gaussian bell shape (inverse parabola). The time $t=0$
corresponds to the end of the quench ramp.

Fig. 12 shows experimental results for the phase difference distribution
function, initially and during the dynamics. More specifically, the histograms
correspond to phase differences $\Delta\varphi$ between all pairs of points in
the middle part of the system at the same distance $r=18\delta z$ where
$\delta z=1.95\mu m$ is the pixel size, at the same time $t$. We observe that
the initial distributions are visibly non-Gaussian in scans 3 and 4: There is
a central peak around $\varphi\approx 0$, but also long tails extending far
from the peak whose presence renders the distribution non-Gaussian. This is
also clear in the log-scale histograms of Fig. 12 which display significant
deviations from the inverse parabola of a Gaussian distribution.

Figure 13: Time evolution of the variance (_left_) and kurtosis (_right_) of
the phase difference distribution as a function of time $t$, in the four
different scans. The variance increases, while the kurtosis tends to the
Gaussian value (i.e., to 3).

The central peak can be intuitively understood by the presence of a strong
energetic penalty on phase fluctuations. Because the sine-Gordon potential is
of the form $\propto\cos(\hat{\varphi}(z))$, the weight of the probability
distribution of the phase values is expected to be concentrated around integer
multiples of $2\pi$, reflecting the presence of soliton configurations. The
presence of smooth jumps by $2\pi$ in some of the observed phase profiles is a
clear demonstration of soliton configurations in the physical system which are
a typical characteristic of the sine-Gordon model. These manifest themselves
as satellite peaks in the phase distribution located at $\Delta\phi\approx\pm
2\pi$. Such peaks have indeed been observed when fast cooling is used for the
preparation of the quantum states [22]. In the present case, however, where
the initial states are prepared by slow cooling, to a good approximation they
correspond to thermal equilibrium states of the sine-Gordon model in which
soliton configurations turn out to be strongly suppressed for the temperatures
of the experiment. Nevertheless, the non-parabolic shape of the cosine
interaction is reflected in the presence of the non-Gaussian long tails in the
phase histograms of the initial states.

Over time, the deviations from the Gaussian shape relax: As we see in Fig. 12,
the probability distributions remain centered around $\varphi\approx 0$ (which
is a gauge choice implemented by referencing the profiles) but becomes
substantially wider, gradually converging towards a profile closely resembling
a Gaussian one. This is quantitatively demonstrated by the time evolution of
the moments (the variance and kurtosis) of the phase distribution, as shown in
Fig. 13.

The dynamics of the phase distributions can be explained in the same way as
the dynamics of the $M^{(4)}$ measure. In the TLL model, which effectively
describes the dynamics after the interaction quench, the non-Gaussian features
of the initial phase distributions are diluted in the dominant Gaussian bath
of the density fluctuations. Additionally, after the quench the phase locking
due to the cosine potential is removed resulting in the broadening of the
phase distribution. In more detail, because of the phase space rotation the
initial phase field is mixed with the initial density field and since the