latter shows purely Gaussian fluctuations, is independent from the former and
exhibits a much larger variance, the resulting distribution broadens with
time. This can be seen from (15) expressed in coordinate space

$\displaystyle\hat{\varphi}(x,t)$
$\displaystyle=\frac{1}{2}\left(\hat{\varphi}(x+ct)+\hat{\varphi}(x-ct)\right)+\frac{\pi}{2K}\int_{x-ct}^{x+ct}\delta\hat{\varrho}(x^{\prime})\,\mathrm{d}x^{\prime}.$
(44)

As we see, the fluctuations of the phase field at any time and any position
can be traced back to initial fluctuations of the phase and density fields. In
a large system, these initial fluctuations are independent from each other
because, on the one hand, the phase fluctuation contributions originate from
distant points, on the other hand, the density fluctuations are independent
from the phase ones. In this case, the distribution of $\hat{\varphi}(.,.)$ is
the convolution of those of the independent constituent initial fields. As a
result, given the above characteristics of $\delta\hat{\varrho}$, the field
$\hat{\varphi}(.,.)$ approaches a Gaussian distribution and broadens with
time. The same arguments apply for the statistics of phase differences
$\Delta\varphi$.

###  5.2 Temporal properties of the decaying connected correlations

We can understand the temporal onset of Gaussification in the experiment by
assuming decaying initial correlations of the velocity field in combination
with the CFA properties of the Gaussianity of the density fluctuations and
independence from the phase fluctuations. In Appendix B, we show that under
these assumptions the magnitude of an $n$-point correlation function will be
reduced by factor $2^{-n+1}$. In the case of the $4$-point correlation
functions which are in focus here, this amounts to almost an order of
magnitude reduction of the extent of connected functions in the phase sector.

This reduction of the connected functions is derived from a similar starting
point as the considerations leading up to Gaussification by spatial scrambling
but using only the first pillar, i.e., correlation clustering for the density
and velocity fields, and without assuming validity of the second pillar, i.e.,
dynamical delocalization. The presence of this second pillar which pertains to
dynamics wherein wave-packets delocalize over time and would essentially imply
Gaussification by spatial scrambling is not verified in the experiment. In
contrast, we will now assume non-dispersive TLL dynamics and hence no
delocalization of the correlation wave-fronts which remain concentrated at the
edges of the effective light-cone. Intuitively, this consideration shows how
the absence of the delocalization pillar affects the physics that would be
seen if it was present.

In the appendix, we make the following intuition more precise. While in TLL
dynamics there is no delocalization, any local operator will split into two
parts which stay local but propagate in opposite directions to the left and
right of the initial position. If we then consider a 4-point function of the
velocity fields

$\displaystyle C^{u}(z,t=0)=\langle{\hat{u}}(z,0)^{4}\rangle_{\text{con}},$
(45)

then after time $t$ we will find that only $1/8$ of the possible correlation
functions resulting from this simple splitting have significant correlations.
These are those that precisely overlap in time

$\displaystyle C^{u}(z,t)=\frac{1}{8}\langle{\hat{u}}(z\pm
ct,0)^{4}\rangle_{\text{con}}+O(e^{-ct/\xi})$ (46)

where the second term comprises by all the other contributing terms which are
tightly upper bounded in time and become suppressed as long as the size of the
light-cone becomes larger than the correlation length $ct\gg\xi$.

Figure 14: Schematic explanation of the equilibration time scaling in the
canonical transmutation mechanism. Due to the linear dispersion of the TLL
liquid, the phase field propagates according to the wave equation, meaning
that it depends on initial data in the interior of the past light-cone (a).
However, the phase field is not a physical observable itself: It is only phase
differences with respect to some reference point
$\Delta\phi(x)=\phi(x)-\phi(x_{0})$ (b) or derivatives of the phase field (c)
(corresponding to the hydrodynamic fields like the velocity
${\hat{u}}\sim\partial_{z}\phi$ and density field
$\delta\hat{\varrho}\sim\partial_{t}\phi$) that are well-defined observables.
The propagation of these fields as derived from that of the phase field is
constrained to remain localized in the vicinity of the light-cone edges. The
time evolved correlations between two (or more) such fields relax shortly
after their initial-time spatial supports stop overlapping, which happens when
the light-cone edges distance from each other, thus resulting in a linear
scaling with time for phase differences and exponential scaling for derivative
fields (c.f. Appendix B).

We see that the above arguments result in a prediction for the dynamics of
connected correlation functions of the velocity fields. Namely, we expect a
linear in time decrease to a plateau value that is reduced compared to the
initial correlations. By integrating the velocity correlations this
suppression pre-factor will be translated to a decrease of the connected part
of phase correlations. This explanation is schematically illustrated in Fig.
14 and in more detail in Appendix B.

The temporal profile of linear scaling followed by saturation is a general
feature of the dynamics of correlations in the TLL model. An intuitive
explanation is given by the quasi-particle description of quantum quenches. In
the particular case of TLL dynamics this interpretation says that the quench
creates quasi-particle excitations that travel with the speed of sound $c$ and
spread correlations from their initial positions throughout the system [59,
6]. If the initial correlations are of short range then this mechanism gives
rise to the above described profile and the change from linear scaling to
saturation is sharp, otherwise the temporal profile gets smeared. If the
dispersion relation is not linear as in the TLL, then the quasi-particles
travel with a range of different velocities, which also results in a smearing
of the temporal profile. Moreover, it results in a slower equilibration
process typically characterized by power-law tails, since the slowest quasi-
particles arrive later than the fastest and the transient dynamics lasts
longer, consistently with the discussion in Sec. 2.1. These arguments apply
also in the case of non-Gaussian initial states and higher-order connected
correlation functions [16]. Even though it is harder to estimate the scaling
of $M^{(4)}$ given that it is a spatially integrated measure, the same
qualitative behavior is valid also here.

We hence obtain an interesting answer to the question what happens if the
delocalization pillar of Gaussification by spatial scrambling is missing: We
find that non-Gaussianity can still reduce over time, but it does not tend
towards zero and is rather reduced by a certain factor. It should be pointed
out that this decrease can be significant. Nevertheless, as shown in the
Supplementary Information of Ref. [24] based on numerical simulations in the
classical fields approximation, this decrease is still insufficient to explain
the experimentally observed decrease. Indeed the experimental decrease stops
at a much smaller value that was only possible to explain through the
canonical transmutation mechanism.

Figure 15: _Left:_ Time evolution of $M^{(4)}$ in a theoretical simulation of
the experiment for typical parameters. The simulations are based on a
stochastic method for the construction of thermal states of the classical
sine-Gordon model [27], followed by TLL dynamics. The parameter values of the
theoretical model used here have been estimated from the experimental data as
discussed in detail in [24]. The red line corresponds to one of the non-
Gaussian initial states and its dynamics agrees very well with those observed
in the experiment. The blue line corresponds to the artificial dynamics of the
same initial state if the contribution of density correlations is ignored.
Equivalently, it corresponds to an initial state with the same phase
fluctuations but with the initial density fluctuations set to zero. The dashed
lines correspond to Gaussian fluctuations with finite statistics and set the
bias level. _Right:_ CFA estimate of the ratio of phase over density initial
fluctuations of the lowest cosine mode, calculated in thermal states of the
sine-Gordon model as a function of the coherence factor
$\langle\cos\phi\rangle$, which is monotonically increasing for increasing
coupling $J$. The gray vertical lines indicate the values corresponding to the
experimental initial states (c.f. Fig. 3). We observe that in this parameter
range the estimated ratio varies between 1/2 and close to zero, values
consistent with the tomographic estimates of Fig. 9. For higher modes the
corresponding ratio is always larger than that of the first mode. This
suggests that even for moderately small values of the ratio, the corresponding
time evolution shown on the left exhibits a significant decrease of non-
Gaussianity comparable to the experimental observations.

In Fig. 3, we see that the decay of non-Gaussianity corresponds to a rapid
decrease to a level that is practically indistinguishable from mere Gaussian
bias due to finite sample size. We find that the profile of the decay of the
non-Gaussianity measure $M^{(4)}$ appears to match with the description of a
fast linear decrease to a value much lower than the initial one. The decay
profile can be contrasted with the power-law decay of the upper bound
envelopes obtained via spatial scrambling. Therefore we may conclude that the
decay profile is fully consistent with that of the canonical transmutation
mechanism rather than the slower decay of the spatial scrambling mechanism.
This is further corroborated by theoretical simulations based on the classical
fields approximation as shown in Fig. 15. Using parameter values consistent
with the experimental data to model the initial states and subsequently
applying TLL dynamics on them, we find that the time evolution of non-
Gaussianity, in particular the rapid and almost linear decrease to the bias
value, agrees very well with the experimental observations. Moreover, if we
ignore the contribution of the initial density fluctuations field to the
dynamics, we clearly see that, although the non-Gaussianity decreases from the
initial value also in this case, it does not reach the bias value but stops
decreasing far above it.

What we should stress is that the non-Gaussianity measure
$M^{(4)}={S_{\text{con}}^{(4)}(t)}/{S_{\text{full}}^{(4)}(t)}$ shown in Fig. 3
is a relative measure of non-Gaussianity. The absolute value of the integrated
four-point connected correlation function $S_{\text{con}}^{(4)}(t)$ in the
experiment decreases but does not decay to almost zero. It is partly this
decrease and partly the increase in ${S_{\text{full}}^{(4)}(t)}$ that results
in the drastic decrease of $M^{(4)}$. This is indeed observed and has been
pointed out in Ref. [24]. Eq. (46) indicates a linear in time decrease of the
non-Gaussianity of connected correlation functions in absolute, not relative
size. Note that extracting estimates of higher-order correlations of the
derivative fields from experimental data seems futile due to the stronger
effects from noise and error sources. For this reason, it has not been
possible to use non-Gaussianity measures based on derivative fields in the
experiment.

The analysis of Appendix B can explain these qualitative characteristics of
the temporal profile of the non-Gaussianity measure $M^{(4)}$ as observed in
the experiment and predicted by the numerical simulation of Fig. 15. In
particular, it explains the reduced but nonzero value of $M^{(4)}$ at large
times when the density correlations of the initial state are artificially set
to zero. It also explains why it is only the relative measure $M^{(4)}$ that
decays and not the absolute value $S_{\text{con}}^{(4)}(t)$ of the integrated
four-point connected correlation function. According to the above, the latter
is expected to decrease to a relatively small but considerable fraction of the
initial value. The lesson drawn is that, while TLL dynamics do lead to a
decrease of non-Gaussianity, which is due to the partial mixing of initial
left and right moving field correlations induced by the dynamics, that alone
is not sufficient to explain the experimentally observed decrease, which
instead can be easily explained by the mixing of phase and density
fluctuations according to the canonical transmutation mechanism.

###  5.3 Interpretation of the experimental data

We now turn to a discussion of the experimental findings of Gaussification and
how they relate to the two mechanisms laid out above. Let us first summarize
our findings regarding the occurrence of the conditions of the two mechanisms
in the experiment and then reach what could be called a verdict on which of
the two is mainly responsible for the observed Gaussification. While this
analysis will not come to a fully consistent conclusion that settles the
question in all detail beyond any doubt, the comprehensive analysis put forth
here as well as the substantial data taken will paint a pretty clear picture
on the mechanisms that can be held responsible for Gaussification in the
experiment.

We again would like to start by evaluating the applicability of the spatial
scrambling mechanism for the data at hand. We have argued that the initial
state exhibits clustering (referred to as mechanism 1, pillar 1, in the above
description) only for a class of the physical fields, specifically those that
are space or time derivatives of the phase field $\hat{\varphi}$, which
correspond to the fluid velocity ${\hat{u}}$ and density displacement
$\delta\hat{\varrho}$. These are the fundamental (local) fields of the
hydrodynamic description. A demonstration of the decay of ${\hat{u}}$ initial
correlations with the distance has been shown for the two point function (Fig.
6). The same can be inferred for the $\hat{\varphi}$ initial correlations
based on the tomographic reconstruction, as shown in Fig. 10. The phase field
$\hat{\varphi}$ itself can only be measured through the difference between two
points and therefore its correlations increase with the distance from the
reference point. It is no surprise that this field does not exhibit
clustering. Moreover, from the observation of relatively sharp light-cone
fronts and – more importantly – of strong recurrences of the initial state, we
have to conclude that the experimental dynamics do not exhibit delocalization
(referred to as mechanism 1, pillar 2 in the above description). This is
consistent with the theoretical dynamics of a TLL, which is characterized by a
linear phononic dispersion relation. In the absence of this second pillar, the
first mechanism results in only a mild decrease of non-Gaussianity that can be
attributed to the mixing of left and right moving components of the
hydrodynamic fields ${\hat{u}}$ and $\delta\hat{\varrho}$, as shown in Fig. 15
and discussed in Subsection 5.2.

Having said that, we now turn to elaborating on the canonical transmutation
mechanism. The requirement of quadrature rotation (so mechanism 2 and pillar
1) is automatically satisfied on the basis of the above verification that the
system is subjected to TLL dynamics. Based on the classical fields
approximation, we have argued that the initial state, modelled as a thermal
state of the sine-Gordon model, is expected to exhibit two crucial properties
with respect to the correlations of the two canonical fields.

Firstly, the two canonical sectors should decouple, i.e., the correlations of
the density and the phase should be independent, which at the same time means
that only the latter are non-Gaussian (mechanism 2, pillar 2). Secondly, the
density correlations should be larger than those of the phase (mechanism 2,
pillar 3) for any value of the coupling $J$ and temperature $\beta^{-1}$. Both
of these statements have been verified in the experiment, at least on the
level of second moments, using the tomographic reconstruction of the initial
state: Density-density mode correlations are indeed generally larger than
phase-phase ones, and both are typically much larger than the density-phase
correlations (Fig. 8 and 9). Even though the ratio of density over phase
quadrature correlations is not excessively larger than one in the initial
state, it is still of sufficient magnitude to justify the validity of pillars
2 and 3 and explain the dominance of the Gaussian density correlations in the
subsequent dynamics. This is explicitly verified in the simulation plots of
Fig. 15. We would like to stress that these simulations are based on the
classical field approximation using estimates of the (not directly measurable)
experimental parameters $J$ and $\beta^{-1}$. In addition to the above, the
temporal profile of the observed decay of non-Gaussianity has been shown to be
consistent with simple geometric arguments following from the light-cone
propagation of hydrodynamic fluctuations, which are more rigorously backed by
analytical calculations presented in Appendix B.

Overall, it seems fair to say that the evidence laid out before suggests that
the canonical transmutation mechanism is primarily at play in the experiment
and constitutes the underlying explanation of the decay of non-Gaussian
correlations. The spatial scrambling mechanism also contributes a partial
decrease of non-Gaussianity, despite the absence of dynamical delocalization,
but that alone is not sufficient to explain the observed decrease
quantitatively, without the additional properties giving rise to the canonical
transmutation mechanism.

##  6 Conclusions

To conclude, in this work, we have observed the dynamical emergence of
Gaussianity from an initially non-Gaussian state and comprehensively discussed
the mechanism responsible for it. We explained the experimental observations
and findings with a simple model of canonical rotation and mixing of phononic
modes and properties of the initial state. Our analysis stresses the
importance of the, often overlooked, higher-order correlations in
characterizing and identifying not only equilibrium states but also dynamical
mechanisms in quantum systems.

Future experiments will aim at further exploiting the potential of
controllable space and time dependent cold-atom traps based on digital
micromirror devices [56, 57], which is what made the implementation of box
shaped traps and the observation of recurrences of non-Gaussianity possible.
Given the high levels of control that are technologically reached in the
experimental platform considered here, in conjunction with the perspective of
reaching regimes beyond the effectively classical description, this set-up can
be regarded as a promising _dynamical quantum simulation_ [60, 61]. Such
dynamical quantum simulations provide new insights into the dynamics of
interacting quantum systems, in instances beyond the scope of known classical
simulation methods. The understanding of the precise dynamical mechanisms at
work in this context is also important for realizing _quantum field machines_
[57], an idea that involves ultra-cold continuous atomic systems to perform
quantum thermodynamic tasks and to treat them as instances of thermodynamic
machines in situations in which quantum effects are expected to play a role.

It would be interesting to investigate the effect of particle statistics on
the two mechanisms studied here. Due to the boson-fermion correspondence in
one spatial dimension the effective description of the experimental system can
be expressed equivalently in terms of either bosonic or fermionic degrees of
freedom and in fact the quasi-fermionic nature of excitations has been
demonstrated in a recent experiment [62]. This suggests that, at least the
canonical transmutation mechanism, where locality plays no role, should be
relevant also in fermionic systems. Even though this mechanism is new and has
not been studied before under more general settings, we indeed expect it to be
applicable to certain cases of interaction quenches in superconducting matter,
spinful fermions, spin ladders, or other systems involving two or more types
of fermions that get mixed by the dynamics, but are initially decoupled and
not all are interacting. In such cases the non-Gaussianity of initial
correlations would oscillate between the different types, similarly to the
present problem. This mechanism is expected to be more clearly distinguishable
from the spatial scrambling mechanism in the case of critical dynamics and for
smeared observables, which depend more on the long-wavelength excitations that
evolve following a linear dispersion relation in this case. The spatial
scrambling mechanism, on the other hand, has been shown to work equally well
for bosonic or fermionic degrees of freedom and spin chains that can be mapped
to non-interacting lattice fermions [7], even though the non-local nature of
this mapping (Jordan–Wigner transformation) clearly plays a non-trivial role
in the validity of the relevant conditions.

Another question for future study is the possible effects of topological
excitations (sine-Gordon solitons) on the two Gaussification mechanisms. As
discussed earlier, solitons are rare in the initial state ensembles of scans
of the present analysis, but earlier experiments [22] have demonstrated
significant presence of solitons in states prepared through a fast non-
equilibrium process (for example, fast evaporative cooling) within the coupled
double well regime. Given that the validity of some of the pillars of the two
mechanisms is shown only in equilibrium and can be affected by the presence of
initial solitons, it is unclear what to expect in such a case.

On a more conceptual level, the work presented here can be seen as a
comprehensive discussion of the mechanisms of the emergence of Gaussian
correlations in physical systems, a type of correlations that is ubiquitous in
physics, to say the least. Our theoretical study is matched and underpinned by
a body of fresh experimental data that exemplify Gaussification in time in a
clear-cut fashion. The diagnostic tools that a tomographic recovery offers
provide novel insight into the precise mechanism that is at work here. It is
our hope that the present work inspires further studies on the emergence of
apparent equilibrium in non-equilibrium quantum dynamics, a field of research
at the interface of strongly correlated quantum systems and quantum field
theory, of quantum information theory, and of statistical physics.

## Acknowledgements

Joint work at FU Berlin and TU Wien has been supported by the DFG (FOR 2724 on
‘Thermal machines in the quantum world’ and CRC 183) and the FQXi on ‘Fueling
quantum field machines with information’, for which it constitutes important
inter-node preparatory theoretical-experimental work in the development of
quantum field machines. FUB has also received funding from the European
Union’s Horizon2020 research and innovation programme under grant agreement
No. 817482 (PASQuanS) on programmable quantum simulators. This work touches
also upon identifying platforms of programmable cold atomic quantum
simulators, as being funded by BMBF (FermiQP), and on methods of verification
of quantum simulators, as funded by the Munich Quantum Valley (K8). It has
also been supported by the DFG/FWF Collaborative Research Centre ‘SFB 1225
(ISOQUANT)’ and the ESQ Discovery Grant ‘Emergence of physical laws: From
mathematical foundations to applications in many body physics’ of the Austrian
Academy of Sciences (ÖAW). F. C., F. S. M., B. R., J. Sabino and T. S.
acknowledge support by the Austrian Science Fund (FWF) in the framework of the
Doctoral School on Complex Quantum Systems (CoQuS). T. S. acknowledges support
by the Max Kade Foundation through a postdoctoral fellowship. J. Sabino
acknowledges support from Fundação para a Ciência e a Tecnologia (Portugal)
through Project No. UIDB/EEA/50008/2020 and from the DP-PMI and FCT
(Portugal). S. S. has been also supported by the Slovenian Research Agency
(ARRS) under grant QTE (N1-0109) and by the ERC Advanced Grant OMNES (694544).
J. E., M. G., J. Schmiedmayer and S. S. thank the Erwin Schrödinger Institute
for its hospitality and support under the programme ‘Quantum Simulation–from
Theory to Application’ (LCW 2019).

##  Appendix A Tomonaga-Luttinger liquid description of the experimental
system

The theoretical description of the experimental system is based on
bosonization or Tomonaga-Luttinger liquid (TLL) theory. The system consists of
a cold atomic gas confined in two parallel one-dimensional traps at short
distance. The two components of the gas are coupled to each other with their
coupling being controlled by the height of the barrier between the two traps.
The system is therefore described by the Hamiltonian

$\displaystyle H$
$\displaystyle=\sum_{i=1,2}\biggl{[}\frac{\hbar^{2}}{2m}\int\mathrm{d}x\,\partial_{x}\Psi_{i}^{\dagger}(x)\partial_{x}\Psi_{i}(x)$
$\displaystyle+\int\mathrm{d}x\mathrm{d}x^{\prime}\,V(x-x^{\prime})\Psi_{i}^{\dagger}(x)\Psi_{i}(x)\Psi_{i}^{\dagger}(x^{\prime})\Psi_{i}(x^{\prime})$
$\displaystyle+\int\mathrm{d}x\,(V_{\text{ext}}(x)-\mu)\Psi_{i}^{\dagger}(x)\Psi_{i}(x)\biggl{]}$
$\displaystyle-\hbar
J\int\mathrm{d}x\,\big{(}\Psi_{1}^{\dagger}(x)\Psi_{2}(x)+\Psi_{2}^{\dagger}(x)\Psi_{1}(x)\big{)}$
(47)

where $\Psi_{i}(.)$ is a two component boson field, $V(.)$ the inter-particle
interaction, $V_{\text{ext}}(.)$ the external trap potential, $\mu$ the
chemical potential and $J$ the tunnelling coupling between the two components.
The inter-particle interaction is practically point-like
$V(x-x^{\prime})=\tfrac{1}{2}g\delta(x-x^{\prime})$ for points $x$ and
$x^{\prime}$. The confining trap is in general inhomogeneous in the
longitudinal direction (parabolic $V_{\text{ext}}(x)$), although homogeneous
box traps have also been used in the experiment. In both cases the edges of
the system $x=\pm L/2$ are characterized by vanishing particle current at all
times. Let us assume for the moment that the trap is homogeneous.

In the _bosonization_ or TLL description [63], each of the bosonic fields is
expressed in terms of density and phase

$\displaystyle\hat{\Psi}^{\dagger}(x)=\sqrt{\hat{\rho}(x)}\mathrm{e}^{\mathrm{i}\hat{\phi}(x)}$
(48)

and the density $\rho_{i}(.)$ is represented as

$\displaystyle\hat{\rho}(x)=\left(n(x)-\frac{1}{\pi}\partial_{x}\hat{\theta}\right)\sum_{\ell=-\infty}^{+\infty}\exp[2\ell\pi\mathrm{i}(\textstyle\int^{x}n(x^{\prime})\mathrm{d}x^{\prime}-\frac{1}{\pi}\hat{\theta}(x))],$

where $n(.)$ is the average density. The auxiliary field $\hat{\theta}$
expresses local deviations of the density from the average value, with
$\partial_{x}\hat{\theta}$ corresponding to long-wavelength density
fluctuations and
$\exp[2\ell\pi\mathrm{i}(\int^{x}n(x^{\prime})\mathrm{d}x^{\prime}-\frac{1}{\pi}\hat{\theta}(x))]$
corresponding to short-wavelength kinks at the positions of the particles [63,
47]. As long as we are interested in long-wavelength density fluctuations
only, we can write

$\displaystyle\delta\hat{\rho}(x)\coloneqq\hat{\rho}(x)-n(x)\approx-\frac{1}{\pi}\partial_{x}\hat{\theta},$
(49)

However, it should be kept in mind that the short-wavelength kinks are also
low-energy excitations like the long-wavelength fluctuations and can play a
significant role in the dynamics. From the bosonic commutation relations of
$\hat{\Psi}^{\dagger}_{i}(.)$, it can be shown that $\phi,\theta$ obey the
canonical commutation relations

$\displaystyle[\delta\hat{\rho}_{i}(x),\hat{\phi}_{j}(y)]=\mathrm{i}\delta_{i,j}\delta(x-y).$
(50)

In the same long-wavelength approximation, the particle current $\hat{j}(.)$
is given by

$\displaystyle\hat{j}(x)\coloneqq-\mathrm{i}\hbar\left(\hat{\Psi}^{\dagger}(x)\partial_{x}\hat{\Psi}(x)-\partial_{x}\hat{\Psi}^{\dagger}(x)\hat{\Psi}(x)\right)\approx\hbar
n(x)\partial_{x}\hat{\phi}(x),$ (51)

from which we can recognize $\hbar\partial_{x}\hat{\phi}(.)$ as playing the
role of the local velocity field. The long-wavelength fields
$\delta\hat{\rho}\sim\partial_{x}\hat{\theta}$ and
$\hat{j}\sim\partial_{x}\hat{\phi}$ representing the particle density and
current respectively are the two fundamental local fields in the hydrodynamic
description of the quantum liquid.

From now on we focus on the homogeneous case where the mean density $n$ is
assumed to be constant in space (and also in time). Replacing the bosonic
field ${\Psi}$ in the Hamiltonian (47) using the above representation,
expanding in powers of $\phi,\delta{\rho}$ and keeping only quadratic and
lowest gradient terms, we obtain the Hamiltonian

$\displaystyle H$ $\displaystyle=H_{1}+H_{2}-\hbar
J\sqrt{n_{1}n_{2}}\int\mathrm{d}x\cos({\phi_{2}(x)-\phi_{1}(x)}),$ (52)

where $H_{1},H_{2}$ are of the form of the TLL Hamiltonian

$\displaystyle H_{LL}$
$\displaystyle=\int\mathrm{d}x\biggl{[}\frac{\hbar^{2}}{2m}n\left(\partial_{x}\hat{\phi}_{i}(x)\right)^{2}+\frac{1}{2}g\delta\hat{\rho}_{i}(x)^{2}\biggr{]}.$
(53)

The Hamiltonian (52) provides a low-energy description of the system. Since
the system is symmetric under interchange of the two components, introducing
the symmetric
$\phi_{s}=\phi_{1}+\phi_{2},\delta\rho_{s}=\tfrac{1}{2}(\delta\rho_{1}+\delta\rho_{2})$
and anti-symmetric fields
$\varphi=\phi_{2}-\phi_{1},\delta\varrho=\tfrac{1}{2}(\delta\rho_{2}-\delta\rho_{1})$,
the Hamiltonian decouples into independent symmetric and anti-symmetric parts
$H=H_{s}+H_{a}$. The symmetric part $H_{s}$ is a TLL Hamiltonian for the
symmetric fields $\hat{\phi}_{s}$ and $\delta\hat{\rho}_{s}$

$\displaystyle H_{LL}$
$\displaystyle=\int\mathrm{d}x\biggl{[}\frac{\hbar^{2}}{4m}n\left(\partial_{x}\hat{\phi}_{s}(x)\right)^{2}+g\delta\hat{\rho}_{s}(x)^{2}\biggr{]}$
(54)

while the anti-symmetric part is

$\displaystyle
H_{sG}=\int\mathrm{d}x\biggl{[}\frac{\hbar^{2}}{4m}n\left(\partial_{x}\hat{\varphi}(x)\right)^{2}+g\delta\hat{\varrho}(x)^{2}\biggr{]}-2\hbar
Jn\int\mathrm{d}x\cos\hat{\varphi}(x),$ (55)

where $\hat{\varphi}$ and $\delta\hat{\varrho}$ are the antisymmetric or
simply the relative phase and density fluctuation fields. The expression in
the last equation is the sine-Gordon Hamiltonian. We see that the coupling
between the two components of the gas plays the role of a Josephson junction
corresponding to a cosine self-interaction of the relative phase field.

From now on, our focus will be on the relative density and phase fields
described by Eq. (55). Due to the presence of interaction, the ground and
thermal states of this model are non-Gaussian in terms of the canonical fields
$\hat{\varphi}(.),\delta\hat{\varrho}(.)$. By setting the barrier height to a
large value, the two components are decoupled, i.e., $J\to 0$ and the sine-
Gordon reduces to

$\displaystyle H_{LL}$
$\displaystyle=\int\mathrm{d}x\biggl{[}\frac{\hbar^{2}}{4m}n\left(\partial_{x}\hat{\varphi}(x)\right)^{2}+g\delta\hat{\varrho}(x)^{2}\biggr{]}$
(56)

which is of the standard TLL form

$\displaystyle\hat{H}_{\text{LL}}=\frac{\hbar c}{2}\int
dx\,\left(\frac{\pi}{K}\left(\delta\hat{\varrho}(x)\right)^{2}+\frac{K}{\pi}(\partial_{x}\hat{\varphi}(x))^{2}\right),$
(57)

where

$c:=\sqrt{\frac{gn}{m}}$ (58)

is the _speed of sound_ and

$K:=\frac{\hbar\pi}{2}\sqrt{\frac{n}{mg}}$ (59)

is the _Luttinger parameter_. In the opposite limit of small barrier height,
the two components are strongly coupled, i.e., $J$ is large, and the
equilibrium properties are determined by the parabolic approximation of the
cosine interaction, that is, the KG Hamiltonian

$\displaystyle
H_{KG}=\int\mathrm{d}x\biggl{[}\frac{\hbar^{2}}{4m}n\left(\partial_{x}\hat{\varphi}(x)\right)^{2}+g\delta\hat{\varrho}(x)^{2}\biggr{]}+\hbar
Jn\int\mathrm{d}x\,\hat{\varphi}^{2}(x).$ (60)

The phonon excitations are now massive with a mass proportional to $\sqrt{J}$.

##  Appendix B Tomonaga-Luttinger liquid dynamics of correlations in a
homogeneous infinite system

In this section, we are going to calculate the large time asymptotics of phase
correlations in a thermodynamically large system based on the assumptions of
i) Tomonaga-Luttinger liquid dynamics, ii) initial clustering of correlations,
iii) homogeneity of the system and iv) the validity of classical field
approximation (CFA) for the initial state which as discussed in subsections
3.2 and 3.2 means that the density fluctuation operator has vanishing
connected correlation functions and vanishing correlations with the phase
operator.

Let us begin by expressing the time evolved phase correlations in terms of
initial correlations using the Heisenberg picture and the fact that TLL
dynamics is Gaussian. Specifically, the Hamiltonian describing the dynamics is
assumed to be (56), therefore the Heisenberg equation of motion for the phase
field is the wave equation

$\displaystyle\partial_{t}^{2}\hat{\varphi}(x,t)-c^{2}\partial_{x}^{2}\hat{\varphi}(x,t)=0$
(61)

with the general solution in infinite space and for arbitrary initial
conditions $\hat{\varphi}(x)=\hat{\varphi}(x,0)$,
$\delta\hat{\varrho}(x)=\delta\hat{\varrho}(x,0)$, given by the D’Alembert
solution

$\displaystyle\hat{\varphi}(x,t)$
$\displaystyle=\frac{1}{2}\left(\hat{\varphi}(x+ct)+\hat{\varphi}(x-ct)\right)+\frac{\pi}{2K}\int_{x-ct}^{x+ct}\delta\hat{\varrho}(x^{\prime})\,\mathrm{d}x^{\prime},$
(62)

where we have used
$\partial_{t}\hat{\varphi}|_{t=0}=(c\pi/K)\,\delta\hat{\varrho}$. As already
mentioned, given that the phase field is only measurable in the form of
differences between two points, we will need to consider either the non-local
phase difference field with respect to some reference point or the local phase
derivative field. Let us focus first on the latter from which we can derive
the former by integration. From the last relation, we see that the time
evolution of ${\hat{u}}=\partial_{x}\hat{\varphi}$ is given by

$\displaystyle{\hat{u}}(x,t)=\partial_{x}\hat{\varphi}(x,t)$
$\displaystyle=\frac{1}{2}\left({\hat{u}}(x+ct)+{\hat{u}}(x-ct)\right)+\frac{\pi}{2K}\left(\delta\hat{\varrho}(x+ct)-\delta\hat{\varrho}(x-ct)\right).$
(63)

The time evolution of ${\hat{u}}=\partial_{x}\hat{\varphi}$ can be
characteristically expressed as the sum of _left-_ and _right-moving_ fields
defined as

$\displaystyle{\hat{u}}(x,t)$
$\displaystyle=\hat{\psi}_{+}(x,t)+\hat{\psi}_{-}(x,t)$ (64)

where

$\displaystyle\hat{\psi}_{\pm}(x,t)\coloneqq\frac{1}{2}\left({\hat{u}}(x,t)\mp\frac{\pi}{K}\delta\hat{\varrho}(x,t)\right)\
.$ (65)

Their name originates from the observation that they follow the simple time
evolution $\hat{\psi}_{\pm}(x,t)=\hat{\psi}_{\pm}(x\pm ct,0)$. From Eq. (63),
it is easy to derive time evolved phase correlation functions of any order
$\left\langle\prod_{i=1}^{n}{\hat{u}}(x_{i},t)\right\rangle$ from initial
ones. All we have to do is substitute (63) and expand to express the result as
a sum of initial velocity and density correlations. The same applies to
connected correlation functions, since they are multi-linear with respect to
the fields. In this way we see that the time evolution of correlations between
a set of points is nothing but the result of mixing initial correlations
originating from the past light-cone projections of these points.

This expansion simplifies significantly when we take into account the
classical field approximation and the property of clustering of correlations
which constrain the initial state. Below we state more formally the precise
definition that we need for our argument.

######  Definition 1 (Classical field approximation (CFA)).

We say that a state satisfies the CFA property, if

$\displaystyle\left\langle\prod_{j=1}^{n}\delta\hat{\varrho}(x_{j})\right\rangle_{\mathrm{con}}=0,\quad\mathrm{for}\;n>2\
$ (66)

and

$\displaystyle\left\langle\prod_{i=1}^{n_{1}}{\hat{u}}(x_{i})\prod_{j=1}^{n_{2}}\delta\hat{\varrho}(x_{j})\right\rangle_{\text{con}}=0,\quad\mathrm{for\,all}\;n_{1},n_{2}>0\
.$ (67)

######  Definition 2 (Clustering of correlations).

We say that a state has exponentially clustering correlations if

$\displaystyle\left\langle\prod_{i=1}^{n_{1}}{\hat{u}}(x_{i})\prod_{j=1}^{n_{2}}\delta\hat{\varrho}(y_{j})\right\rangle_{\mathrm{con}}=O(e^{-\sum_{i<j}|x_{i}-x_{j}|/\xi-\sum_{i<j}|y_{i}-y_{j}|/\xi}),\quad\mathrm{for}\;n_{1},n_{2}>0,$
(68)

where $\xi>0$ is a correlation length. If the bound is a polynomially decaying
function then we speak of polynomially clustering correlations.

Let us consider a higher-order connected correlation function with $n>2$ and
unequal positions $x_{i}\neq x_{j}$ unless $i=j$

$\displaystyle\left\langle\prod_{i=1}^{n}{\hat{u}}(x_{i},t)\right\rangle_{\mathrm{con}}$
$\displaystyle=\frac{1}{2^{n}}\sum_{\begin{subarray}{c}\sigma\in\\{-1,1\\}^{\times
n}\\\ \mu\in\\{0,1\\}^{\times
n}\end{subarray}}\left\langle[{\hat{u}}(x_{1}+\sigma_{1}ct)]^{\mu_{1}}[\frac{\pi\sigma_{1}}{K}\delta\hat{\varrho}(x_{1}+\sigma_{1}ct)]^{1-\mu_{1}}\ldots\right.$
(69)
$\displaystyle\left.\times[{\hat{u}}(x_{n}+\sigma_{n}ct)]^{\mu_{n}}[\frac{\pi\sigma_{n}}{K}\delta\hat{\varrho}(x_{n}+\sigma_{n}ct)]^{1-\mu_{n}}\right\rangle_{\mathrm{con}}.$

Notice that this formula is already simpler than the general case of Gaussian
dynamics (11) where the Green’s functions characterizing the propagation are
not simply localized ($\delta$-like) functions. This property is directly
related to the fact that under TLL dynamics wave-packets do not spread. For
$n=2$, we have

$\displaystyle\left\langle{\hat{u}}(x_{1},t){\hat{u}}(x_{2},t)\right\rangle$
$\displaystyle=\frac{1}{4}\sum_{\sigma_{1},\sigma_{2}=\pm}\bigg{[}\left\langle{\hat{u}}(x_{1}+\sigma_{1}ct){\hat{u}}(x_{2}+\sigma_{2}ct)\right\rangle$
(70)
$\displaystyle+\left(\frac{\pi}{K}\right)^{2}\left\langle\delta\hat{\varrho}(x_{1}+\sigma_{1}ct)\delta\hat{\varrho}(x_{2}+\sigma_{2}ct)\right\rangle\bigg{]}.$

Here, we see explicitly that whenever $\sigma_{1}=\sigma_{2}=\sigma$ then the
time dependent correlation function is just a rigid translation of the initial
one between the points $x_{1}+\sigma ct$ and $x_{2}+\sigma ct$.

The field operators in the correlation functions can be freely reshuffled for
$n>2$ because any field commutators resulting from this reordering do not
contribute to the connected correlations. After ordering the operators, we use
the CFA property to evaluate the terms involving only $\hat{\varphi}(.)$
fields or both $cf(.)$ and $\hat{\varphi}(.)$ fields, which simply vanish
since $n>2$ and we are left with the terms involving ${\hat{u}}(.)$ fields
only

$\displaystyle\left\langle\prod_{i=1}^{n}{\hat{u}}(x_{i},t)\right\rangle_{\mathrm{con}}$
$\displaystyle=\frac{1}{2^{n}}\sum_{\begin{subarray}{c}\sigma\in\\{-1,1\\}^{\times
n}\\\ \mu\in\\{0,1\\}^{\times
n}\end{subarray}}\left\langle\prod_{i=1}^{n}[{\hat{u}}(x_{i}+\sigma_{i}ct)]^{\mu_{i}}\prod_{i=1}^{n}[\frac{\sigma_{i}\pi}{K}\delta\hat{\varrho}(x_{i}+\sigma_{i}ct)]^{1-\mu_{i}}\right\rangle_{\mathrm{con}}$
$\displaystyle=\frac{1}{2^{n}}\sum_{{\sigma\in\\{-1,1\\}^{\times
n}}}\left\langle{\hat{u}}(x_{1}+\sigma_{1}ct)\ldots{\hat{u}}(x_{n}+\sigma_{n}ct)\right\rangle_{\mathrm{con}}.$
(71)

Further simplification applies when we consider the large time limit of
correlations and exploit the clustering of initial correlations of local
fields. In this limit the initial correlations involved in the above
expansions are non-trivial only if all $\sigma_{i}$ have the same sign. This
is because the fields with $\sigma_{i}>0$ propagate to the right and those
with $\sigma_{i}<0$ to the left and therefore become quickly uncorrelated
since they originate from distant points. For sufficiently long times, the
light-cone separation of these points $ct$ is substantially longer than the
constant diameter of the initial points $\text{diam}\\{x_{i}\\}$ and the
correlation length $\xi>0$. We can hence bound asymptotically the connected
function as

$\displaystyle\lim_{t\to\infty}\left\langle\prod_{i=1}^{n}{\hat{u}}(x_{i}+\sigma_{i}ct)\right\rangle_{\text{con}}$
$\displaystyle=\lim_{t\to\infty}O(e^{-ct/\xi})=0,\quad\text{for signs
}\sigma_{i}\text{ not all equal}.$ (72)

Therefore we can simplify the connected correlation functions accordingly

$\displaystyle\left\langle\prod_{i=1}^{n}{\hat{u}}(x_{i},t)\right\rangle_{\mathrm{con}}$
$\displaystyle=\frac{1}{2^{n}}\left\langle\prod_{i=1}^{n}{\hat{u}}(x_{i}+ct)\right\rangle_{\mathrm{con}}+\frac{1}{2^{n}}\left\langle\prod_{i=1}^{n}{\hat{u}}(x_{i}-ct)\right\rangle_{\mathrm{con}}.$
(73)

This means that for late times, what is left are correlations between points
at fixed distances. Assuming also translational invariance of the initial
state, the two terms above are equal

$\displaystyle\lim_{t\to\infty}\left\langle\prod_{i=1}^{n}{\hat{u}}(x_{i}+ct)\right\rangle_{\text{con}}$
$\displaystyle=\lim_{t\to\infty}\left\langle\prod_{i=1}^{n}{\hat{u}}(x_{i}-ct)\right\rangle_{\text{con}}=\left\langle\prod_{i=1}^{n}{\hat{u}}(x_{i})\right\rangle_{\text{con}}\
.$ (74)

This allows us to drop the dependence on signs and we finally obtain the
following result for the large time correlations under the above assumptions

$\displaystyle\left\langle\prod_{i=1}^{n}{\hat{u}}(x_{i},t)\right\rangle_{\mathrm{con}}$
$\displaystyle=\frac{1}{2^{n-1}}\left\langle\prod_{i=1}^{n}{\hat{u}}(x_{i})\right\rangle_{\mathrm{con}}+O(e^{-ct/\xi})$
(75)

for $n>2$ and for $n=2$

$\displaystyle\lim_{t\to\infty}\left\langle{\hat{u}}(x_{1},t){\hat{u}}(x_{2},t)\right\rangle$
$\displaystyle=\frac{1}{2}\bigg{[}\left\langle{\hat{u}}(x_{1}){\hat{u}}(x_{2})\right\rangle+\left(\frac{\pi}{K}\right)^{2}\left\langle\delta\hat{\varrho}(x_{1})\delta\hat{\varrho}(x_{2})\right\rangle\bigg{]}.$

In particular, for an auto-correlation function, we have

$\displaystyle\left\langle{\hat{u}}(x,t)^{n}\right\rangle_{\mathrm{con}}$
$\displaystyle=\frac{1}{2^{n-1}}\left\langle{\hat{u}}(x,0)^{n}\right\rangle_{\mathrm{con}}+O(e^{-ct/\xi}).$
(76)

Figure 16: Geometric explanation of the temporal profile of cumulants of the
phase difference field
$\Delta\hat{\varphi}=\hat{\varphi}(x,t)-\hat{\varphi}(x_{0},t)$. (a) Seen from
a space-time perspective, the splitting of $\Delta\hat{\varphi}$ into two
components (left and right moving) can be depicted by two stripes of different
color. For times $t<|x-x_{0}|/2c$ the two projections on the $t=0$ line still
overlap with each other, while for $t>|x-x_{0}|/2c$ they do not. (b) To
compute the 2nd order cumulant of
$C_{2}(x,t)=\left\langle\Delta\hat{\varphi}(x,t)^{2}\right\rangle_{\mathrm{con}}$
we integrate
$\left\langle{\hat{u}}(y_{1},t){\hat{u}}(y_{2},t)\right\rangle_{\mathrm{con}}$
over the variables $y_{1}$ and $y_{2}$ in the overlapping region of the
corresponding stripes which are shifted from the original interval $[x_{0},x]$
by $\pm ct$. Given that this correlation function is of short range, i.e.,
decays quickly away from the diagonal $y_{1}=y_{2}$, the integration gives
simply a value proportional to the length of the intersections of these
regions and the diagonal (red lines). The stripes move apart with speed $2c$,
so the overall length of the intersection intervals at first decreases
linearly with time and eventually stops changing. (c) The resulting temporal
profile of $C_{2}(x,t)$ exhibits a linear decrease to a plateau value, which
is half of the initial one. The exponential tails of the correlation function
away from the diagonal have the effect of smearing the otherwise sharp
temporal profile [59]. Following the same arguments applied to higher order
cumulants $C_{n}(x,t)$ we get a decrease by a factor $1/2^{n-1}$.

In a similar way, we can derive the scaling of phase difference correlations
which turn out to decrease linearly with time until they reach a saturation
value. This is due to a purely geometric reason, as we can see by integrating
the phase derivative correlations. From (73) the integration interval
$[x_{0},x]$ is split into two intervals $[x_{0}-ct,x-ct]$ and
$[x_{0}+ct,x+ct]$ corresponding to the light-cone projections when we trace
the time-evolved correlations to initial correlations. The two intervals
overlap in the sub-interval $[x_{0}+ct,x-ct]$ whose length is linearly
decreasing with time, up to the time $t_{*}=|x_{0}-x|/(2c)$ beyond which they
do not overlap anymore. Since the initial correlations of the phase derivative
are short-range, the integral is simply proportional to the length of the
overlap, as illustrated in Fig. 14 and in more detail in Fig. 16, resulting in
the above described behavior. The large time asymptotic value can be easily
related to the initial one

$\displaystyle\left\langle\Delta\hat{\varphi}(x,t)^{n}\right\rangle_{\mathrm{con}}$
$\displaystyle=\int_{[x_{0},x]^{n}}\prod_{i}^{n}dy_{i}\left\langle\prod_{i}^{n}{\hat{u}}(y_{i},t)\right\rangle_{\mathrm{con}}$
(77)
$\displaystyle=\frac{1}{2^{n}}\int_{([x_{0}-ct,x-ct]\cup[x_{0}+ct,x+ct])^{n}}\prod_{i}^{n}dy_{i}\left\langle\prod_{i}^{n}{\hat{u}}(y_{i},0)\right\rangle_{\mathrm{con}}$
$\displaystyle\to\frac{1}{2^{n-1}}\left\langle\Delta\hat{\varphi}(x,0)^{n}\right\rangle_{\mathrm{con}}\
.$

For the second cumulant the decrease is by a factor 1/2, for the fourth
cumulant by 1/8 and so on. As discussed in Subsection 5.2, the above analysis
can explain the qualitative characteristics of the temporal profile of the
non-Gaussianity measure $M^{(4)}$ in the experiment.

##  Appendix C Deviations from the Tomonaga-Luttinger liquid dynamics

The TLL Hamiltonian provides a very good description of equilibrium properties
of cold atom systems as can be shown by means of _renormalization group_
theory [47]. However, it is less clear how accurate this description is in far
from equilibrium dynamics, like after a quantum quench. As discussed in the
main text, deviations from the TLL model that are irrelevant at equilibrium in
the renormalization group sense may be important out of equilibrium. One class
of deviations comes from the fact that, in passing from the original
Hamiltonian (47) to (54), we ignored powers and gradients of the density-phase
fields of order higher than two. Higher power terms in particular induce an
effective self-interaction of the phonons resulting in non-Gaussian dynamics.
Moreover, some of these corrections correspond to coupling between the
symmetric and anti-symmetric modes, meaning that the dynamics of the anti-
symmetric modes is not completely closed. In Sec. 4.3 of the main text we have
discussed the dynamical effects of deviations from the TLL description and
presented evidence that such deviations are negligible in the experiment. As
we argued, despite its simplicity the Hamiltonian (54) actually captures the
dynamics of the system sufficiently well within the time scales of the
experiment.

Here, we will analyze in more detail two of the main deviations from this
model that are potentially relevant for the Gaussification mechanism. The
first one is related to the presence of inhomogeneity in the system. In
general cold atom gases are inhomogeneous due to the longitudinal trapping
potential, which is typically parabolic. However, in the present experiment a
homogeneous box trap with hard walls was used, so that the atom density was
practically homogeneous. It is still interesting to observe that even in the
parabolic trap case, the time evolution does not induce delocalization as
required for the spatial scrambling mechanism. The second type of deviation
refers to non-linear corrections to the phonon dispersion relation, which are
present due to the higher gradient terms that were ignored in the derivation
of (54). These are quadratic in the density and phase fields and therefore
preserve the Gaussianity of dynamics but induce dispersive spreading of the
phononic excitations. Since these physical effects are relevant in
Gaussification, it is important to estimate their role in the experiment.

###  C.1 Effects of a non-uniform trap potential

The trap inhomogeneity can be taken into account by allowing the density
profile $n$ in (54) to be a function of the spatial coordinate $x$

$\displaystyle\hat{H}_{\text{inhom}}=\int\mathrm{d}x\biggl{[}$
$\displaystyle\frac{\hbar^{2}n(x)}{4m}\left(\partial_{x}\hat{\varphi}(x)\right)^{2}+g\delta\hat{\varrho}(x)^{2}\biggr{]}\
.$ (78)

In the experimental settings, the Luttinger parameter $K$ is sufficiently
large in which case the _Thomas-Fermi approximation_ is applicable. In this
approximation and for a parabolic trap
$V_{\text{ext}}(x)=\tfrac{1}{2}\omega^{2}x^{2}$, the density profile $n(x)$
turns out to be

$n(x)=n_{0}\bigg{(}1-\frac{x^{2}}{R^{2}}\bigg{)}\Theta(R^{2}-x^{2}),$ (79)

where $n_{0}=\tfrac{1}{2}\omega^{2}R^{2}$ is the density at the middle and
$R=\sqrt{2\mu}/\omega$ is the semi-classical (or Thomas-Fermi) radius of the
system, which is half of the system size.

The Hamiltonian (78) corresponds to an inhomogeneous TLL [64]. In the present
special case (weak boson interaction and parabolic trap) an exact solution is
possible by expanding the fields $\hat{\varphi}$ and $\delta\hat{\varrho}$ in
Legendre polynomials instead of the cosine plane waves of the homogeneous
case. The energies of these modes are $E_{n}=\sqrt{n(n+1)}/R$ (instead of
integer multiples of $2\pi/L$). Approximate recurrences are expected to occur
at integer multiples of $T=2\pi R$ (instead of $T=2L$), since at higher mode
numbers all energies are close to integer multiples of $1/R$. The
inhomogeneity of the trap affects significantly the dynamics. A wave packet
initially localized at the centre propagates to the edges following curved
light-cones at a position-dependent speed, as expected for a curved
background. In addition, however, edge effects are present, which grow and
spread inwards to the centre of the trap. Instead of (62), the time evolution
of the phase field is now given by

$\displaystyle\hat{\varphi}(x,t)=\int\mathrm{d}y\;G_{\phi\phi}(x,y,t)\hat{\varphi}(y,0)+\int\mathrm{d}y\;G_{\phi\rho}(x,y,t)\delta\hat{\varrho}(y,0)$
(80)

where the qualitative form of the phase-phase and density-phase propagators
$G_{\phi\phi}(x,y,t)$ and $G_{\phi\rho}(x,y,t)$ in the parabolic trap is
presented in Fig. 19, which should be compared with Figs. 17 and 18 for the
hard-wall box trap with Neumann or Dirichlet boundary conditions,
respectively.

Figure 17: Phase-phase $G_{\phi\phi}(x,y,t)$ (_top_) and density-phase
$G_{\phi\rho}(x,y,t)$ propagators (_bottom_) for a homogeneous box with
Neumann boundary conditions (real space density plots at various times from
zero to half of the recurrence time).  Figure 18: Same as Fig. 17 but for a
homogeneous box with Dirichlet boundary conditions.  Figure 19: Same as Fig.
17 but for a parabolic trap as described by the Thomas-Fermi approximation.

From our analysis we conclude that despite its special characteristics, the
dynamics does not induce delocalization that could justify Gaussification via
wave packet spreading. In contrast we find that for all relevant aspects, the
dynamics up to the maximum time studied in the experiment are unaffected by
edge effects and practically the same as that of a thermodynamically large
homogeneous system.

###  C.2 Effects of a non-linear phonon dispersion

The linear dispersion of the TLL in its standard form (54) is a rather
idealistic property that is far from true in realistic systems. The most
important source of nonlinear dispersion comes from the Galilean form of the
kinetic energy in the original Hamiltonian (47) [65]. In the density-phase
representation this reads

$\displaystyle\hat{H}_{\text{kin}}$
$\displaystyle=\int\mathrm{d}x\biggl{[}\frac{\hbar^{2}n}{4m}\left(\partial_{x}\hat{\varphi}(x)\right)^{2}+\frac{\hbar^{2}}{4mn}\left(\partial_{x}\delta\hat{\varrho}(x)\right)^{2}\biggr{]}\
,$ (81)

meaning that the extra term

$\displaystyle\hat{H}_{\text{disp}}$
$\displaystyle=\frac{\hbar^{2}}{4mn}\int\mathrm{d}x\left(\partial_{x}\delta\hat{\varrho}(x)\right)^{2}\
,$ (82)

should be added to the TLL Hamiltonian (54). The new Hamiltonian is Gaussian
with a different dispersion relation

$E(k)=\hbar|k|\sqrt{c^{2}+\left(\frac{\hbar k}{2m}\right)^{2}}=\hbar
c|k|\sqrt{1+\left(\frac{\xi_{h}k}{2}\right)^{2}}$ (83)

with $c=\sqrt{gn/m}$ as in Eq. (58) and $\xi_{h}=\hbar/\sqrt{gnm}$ is the
_healing length_. This is the Bogoliubov dispersion relation, which is linear
at small $k$ ($E(k)\sim\hbar c|k|$ ) consistently with TLL theory, and
quadratic at large $k$ ($E(k)\sim\hbar^{2}k^{2}/(2m)$). The linear dispersion
of Eq. (54) is an excellent approximation at equilibrium as long as we are
interested in length scales much larger than $\xi_{h}$. In the experiment
phase measurements refer to smeared local fields due to the finite imaging
(spatial) resolution which is of the order of $3\mu m$. On the other hand, the
healing length is estimated to be $\xi_{h}\approx 0.35\mu m$, much smaller
than the smearing length. Therefore for equilibrium states the short distance
effects of the nonlinear dispersion (83) are unobservable and negligible. This
is not necessarily true, however, for the dynamics: As the non-linearity of
the dispersion relation induces spreading of local fields following an
algebraic scaling with time, this spreading effect will eventually become
significant and manifest itself in the measurements, even if these are
restricted to smeared local fields. Nevertheless, by evaluating the spreading
effect in the time scale of the experiment, i.e., up to the recurrence time,
we find that it is still negligible and unimportant.

Figure 20: A reconstruction for the same experimental scan as in Fig. 8 but
with a restriction of the eigen-mode correlations to the diagonal, i.e.,
instead of a reconstruction over unrestricted Gaussian states, here, we
restrict the variational states with an additional product-state constraint in
eigen-mode space. We notice a remarkable stability of the values of
correlations on the diagonal. The absence of off-diagonal eigen-mode
correlations due to the restriction leads to a reduced fidelity of the fit.
While the precision is a little worse, the accuracy of extrapolation of the
time evolution is improved as the restriction forces proximity to steady
states which behave more stably upon extrapolation outside of the input time
window.

##  Appendix D Further tomographic reconstruction plots

In this section we include additional plots of tomographic reconstructions of
one of the experimental scans (Figs. 20, 21 and 22) that provide more details
on the quality of the reconstruction for different choices of the variational
states or of the number of modes used in the tomographic method.

Figure 21: A reconstruction for the same experimental scan as in Fig. 8 but
now allowing for an larger number of eigen-modes $\mathcal{N}=20$. Comparing
with the earlier reconstruction using $\mathcal{N}=10$ eigen-modes we find a
rather good convergence in the sense that the correlations which are well-
resolved in space and time (modes with $k\leq 7$) are accurately reproduced in
both reconstructions. It should be stressed that the time step of the
measurements and the spatial resolution of the read-out camera do not allow
the accurate resolution of higher momentum modes, and the reconstruction
favors diagonal non-squeezed correlations in these modes. The inclusion of
additional modes softens the hard cut-off at high momentum modes leading to
less off-diagonal artifacts.  Figure 22: The real space representation of the
reconstructed density fluctuation covariance matrix for the same experimental
scan as in Fig. 8 but based on $\mathcal{N}=20$ eigen-modes. We see a much
more pronounced localization in coordinate space since the inclusion of higher
momentum modes allows us to observe more short-range structures. Most notably
in comparison to $\mathcal{N}=10$ we observe the typical anti-correlation
stripes close to the diagonal.

## References

  * [1] A. Polkovnikov, K. Sengupta, A. Silva and M. Vengalattore,  _Nonequilibrium dynamics of closed interacting quantum systems_ ,  Rev. Mod. Phys. 83, 863 (2011),  10.1103/RevModPhys.83.863.
  * [2] C. Gogolin and J. Eisert,  _Equilibration, thermalisation, and the emergence of statistical mechanics in closed quantum systems_ ,  Rep. Prog. Phys. 79, 56001 (2016),  10.1088/0034-4885/79/5/056001.
  * [3] R. Vasseur and J. E. Moore,  _Non-equilibrium quantum dynamics and transport: From integrability to many-body localization_ ,  J. Stat. Mech. p. 064010 (2016),  10.1088/1742-5468/2016/06/064010.
  * [4] J. Eisert, M. Friesdorf and C. Gogolin,  _Quantum many-body systems out of equilibrium_ ,  Nature Phys. 11, 124 (2015),  doi:10.1038/nphys3215.
  * [5] M. Cramer, C. M. Dawson, J. Eisert and T. J. Osborne,  _Exact relaxation in a class of nonequilibrium quantum lattice systems_ ,  Phys. Rev. Lett. 100, 030602 (2008),  10.1103/PhysRevLett.100.030602.
  * [6] P. Calabrese and J. Cardy,  _Entanglement and correlation functions following a local quench: a conformal field theory approach_ ,  J. Stat. Mech. p. P10004 (2007),  10.1088/1742-5468/2007/10/P10004.
  * [7] M. Gluza, C. Krumnow, M. Friesdorf, C. Gogolin and J. Eisert,  _Equilibration via Gaussification in fermionic lattice systems_ ,  Phys. Rev. Lett. 117, 190602 (2016),  10.1103/PhysRevLett.117.190602.
  * [8] M. Rigol, V. Dunjko, V. Yurovsky and M. Olshanii,  _Relaxation in a completely integrable many-body quantum system: An ab initio study of the dynamics of the highly excited states of 1D lattice hard-core bosons_ ,  Phys. Rev. Lett. 98, 050405 (2007),  10.1103/PhysRevLett.98.050405.
  * [9] M. Rigol, V. Dunjko and M. Olshanii,  _Thermalization and its mechanism for generic isolated quantum systems_ ,  Nature 452, 854 (2008),  10.1038/nature06838.
  * [10] P. Calabrese, F. H. L. Essler and M. Fagotti,  _Quantum quench in the transverse-field ising chain_ ,  Phys. Rev. Lett. 106, 227203 (2011),  10.1103/PhysRevLett.106.227203.
  * [11] E. Ilievski, J. De Nardis, B. Wouters, J.-S. Caux, F. H. L. Essler and T. Prosen,  _Complete generalized Gibbs ensembles in an interacting theory_ ,  Phys. Rev. Lett. 115, 157201 (2015),  10.1103/PhysRevLett.115.157201.
  * [12] T. Langen, S. Erne, R. Geiger, B. Rauer, T. Schweigler, M. Kuhnert, W. Rohringer, I. E. Mazets, T. Gasenzer and J. Schmiedmayer,  _Experimental observation of a generalized Gibbs ensemble_ ,  Science 348, 207 (2015),  10.1126/science.1257026.
  * [13] R. Haag, R. V. Kadison and D. Kastler,  _Asymptotic orbits in a free fermi gas_ ,  Commun. Math. Phys. 33, 1 (1973),  10.1007/bf01645603.
  * [14] M. Cramer and J. Eisert,  _A quantum central limit theorem for non-equilibrium systems: Exact local relaxation of correlated states_ ,  New J. Phys. 12, 055020 (2010),  10.1088/1367-2630/12/5/055020.
  * [15] A. Bastianello and S. Sotiriadis,  _Cluster expansion for ground states of local Hamiltonians_ ,  Nucl. Phys. B 909, 1020 (2016),  10.1016/j.nuclphysb.2016.06.021.
  * [16] S. Sotiriadis,  _Memory-preserving equilibration after a quantum quench in a one-dimensional critical model_ ,  Phys. Rev. A 94, 031605 (2016),  10.1103/PhysRevA.94.031605.
  * [17] S. Sotiriadis,  _Equilibration in one-dimensional quantum hydrodynamic systems_ ,  J. Phys. A 50, 424004 (2017),  10.1088/1751-8121/aa8aa5.
  * [18] S. Sotiriadis and P. Calabrese,  _Validity of the GGE for quantum quenches from interacting to noninteracting models_ ,  J. Stat. Mech. 2014, P07024 (2014),  10.1088/1742-5468/2014/07/P07024.
  * [19] M. Gluza, J. Eisert and T. Farrelly,  _Equilibration towards generalized Gibbs ensembles in non-interacting theories_ ,  SciPost Phys. 7, 38 (2019),  10.21468/SciPostPhys.7.3.038.
  * [20] C. Murthy and M. Srednicki,  _Relaxation to Gaussian and generalized Gibbs states in systems of particles with quadratic Hamiltonians_ ,  Phys. Rev. E 100, 012146 (2019),  10.1103/PhysRevE.100.012146.
  * [21] T. Ishii and T. Mori,  _Strong eigenstate thermalization within a generalized shell in noninteracting integrable systems_ ,  Phys. Rev. E 100, 012139 (2019),  10.1103/physreve.100.012139.
  * [22] T. Schweigler, V. Kasper, S. Erne, I. Mazets, B. Rauer, F. Cataldini, T. Langen, T. Gasenzer, J. Berges and J. Schmiedmayer,  _Experimental characterization of a quantum many-body system via higher-order correlations_ ,  Nature 545, 323 (2017),  10.1038/nature22310.
  * [23] T. V. Zache, T. Schweigler, S. Erne, J. Schmiedmayer and J. Berges,  _Extracting the field theory description of a quantum many-body system from experimental data_ ,  Phys. Rev. X 10, 011020 (2020),  10.1103/PhysRevX.10.011020.
  * [24] T. Schweigler, M. Gluza, M. Tajik, S. Sotiriadis, F. Cataldini, S.-C. Ji, F. S. Møller, J. Sabino, B. Rauer, J. Eisert and J. Schmiedmayer,  _Decay and recurrence of non-Gaussian correlations in a quantum many-body system_ ,  Nature Phys. 17, 559 (2021),  10.1038/s41567-020-01139-2.
  * [25] M. Gluza, T. Schweigler, B. Rauer, C. Krumnow, J. Schmiedmayer and J. Eisert,  _Quantum read-out for cold atomic quantum simulators_ ,  Comm. Phys. 3, 12 (2020),  10.1038/s42005-019-0273-y.
  * [26] J. Eisert, D. Hangleiter, N. Walk, I. Roth, D. Markham, R. Parekh, U. Chabaud and E. Kashefi,  _Quantum certification and benchmarking_ ,  Nature Rev. Phys. 2, 382 (2020),  10.1038/s42254-020-0186-4.
  * [27] S. Beck, I. E. Mazets and T. Schweigler,  _Nonperturbative method to compute thermal correlations in one-dimensional systems_ ,  Phys. Rev. A 98, 023613 (2018),  10.1103/PhysRevA.98.023613.
  * [28] R. Folman, P. Krüger, D. Cassettari, B. Hessmo, T. Maier and J. Schmiedmayer,  _Controlling cold atoms using nanofabricated surfaces: Atom chips_ ,  Phys. Rev. Lett. 84, 4749 (2000),  10.1103/PhysRevLett.84.4749.
  * [29] T. Schumm, S. Hofferberth, L. M. Andersson, S. Wildermuth, S. Groth, I. Bar-Joseph, J. Schmiedmayer and P. Krüger,  _Matter-wave interferometry in a double well on an atom chip_ ,  Nature Phys. 1, 57 (2005),  10.1038/nphys125.
  * [30] B. Rauer, S. Erne, T. Schweigler, F. Cataldini, M. Tajik and J. Schmiedmayer,  _Recurrences in an isolated quantum many-body system_ ,  Science 360, 307 (2018),  10.1126/science.aan7938.
  * [31] Y. D. van Nieuwkerk, J. Schmiedmayer and F. Essler,  _Projective phase measurements in one-dimensional Bose gases_ ,  SciPost Physics 5, 046 (2018),  10.21468/scipostphys.5.5.046.
  * [32] P. Hauke, M. Lewenstein and A. Eckardt,  _Tomography of band insulators from quench dynamics_ ,  Phys. Rev. Lett. 113, 045303 (2014),  10.1103/PhysRevLett.113.045303.
  * [33] T. Schweigler,  _Correlations and dynamics of tunnel-coupled one-dimensional Bose gases_ ,  Ph.D. thesis, TU Wien,  arxiv:1908.00422 (2019).
  * [34] V. Gritsev, A. Polkovnikov and E. Demler,  _Linear response theory for a pair of coupled one-dimensional condensates of interacting atoms_ ,  Phys. Rev. B 75, 174511 (2007),  10.1103/PhysRevB.75.174511.
  * [35] V. Gritsev, E. Demler, M. Lukin and A. Polkovnikov,  _Spectroscopy of collective excitations in interacting low-dimensional many-body systems using quench dynamics_ ,  Phys. Rev. Lett. 99, 200404 (2007),  10.1103/PhysRevLett.99.200404.
  * [36] T. Schweigler, M. Tajik, F. Cataldini, S.-C. Ji, F. S. Møller, J. Sabino, B. Rauer and J. Schmiedmayer,  _Raw data for ‘Decay and recurrence of non-Gaussian correlations in a quantum many-body system’_ ,  Zenodo (CERN),  10.5281/zenodo.3753155 (2020).
  * [37] I. Kukuljan, S. Sotiriadis and G. Takacs,  _Correlation functions of the quantum sine-Gordon model in and out of equilibrium_ ,  Phys. Rev. Lett. 121, 110402 (2018),  10.1103/physrevlett.121.110402.
  * [38] M. B. Hastings and T. Koma,  _Spectral gap and exponential decay of correlations_ ,  Commun. Math. Phys. 265, 781 (2006),  10.1007/s00220-006-0030-4.
  * [39] M. Kliesch, C. Gogolin, M. J. Kastoryano, A. Riera and J. Eisert,  _Locality of temperature_ ,  Phys. Rev. X 4, 031019 (2014),  10.1103/PhysRevX.4.031019.
  * [40] A. Molnar, N. Schuch, F. Verstraete and J. I. Cirac,  _Approximating Gibbs states of local Hamiltonians efficiently with projected entangled pair states_ ,  Phys. Rev. B 91, 045138 (2015),  10.1103/PhysRevB.91.045138.
  * [41] T. Langen, R. Geiger, M. Kuhnert, B. Rauer and J. Schmiedmayer,  _Local emergence of thermal correlations in an isolated quantum many-body system_ ,  Nature Phys. 9, 640 (2013),  10.1038/nphys2739.
  * [42] V. Gritsev, E. Altman, E. Demler and A. Polkovnikov,  _Full quantum distribution of contrast in interference experiments between interacting one-dimensional bose liquids_ ,  Nature Phys. 2, 705 (2006),  10.1038/nphys410.
  * [43] A. Polkovnikov, E. Altman and E. Demler,  _Interference between independent fluctuating condensates_ ,  Proc. Natl. Ac. Sc. 103, 6125 (2006),  10.1073/pnas.0510276103.
  * [44] S. Hofferberth, I. Lesanovsky, T. Schumm, A. Imambekov, V. Gritsev, E. Demler and J. Schmiedmayer,  _Probing quantum and thermal noise in an interacting many-body system_ ,  Nature Phys. 4, 489 (2008),  10.1038/nphys941.
  * [45] M. Peskin and D. Schroeder,  _An introduction to quantum field theory_ ,  Frontiers in Physics. Westview Press,  10.1201/9780429503559 (1995).
  * [46] J. Zinn-Justin,  _Quantum field theory and critical phenomena_ ,  Clarendon Press,  10.1093/acprof:oso/9780198509233.001.0001 (1996).
  * [47] T. Giamarchi,  _Quantum physics in one dimension_ ,  Clarendon Press, Oxford,  10.1093/acprof:oso/9780198525004.001.0001 (2004).
  * [48] D. S. Petrov, D. M. Gangardt and G. V. Shlyapnikov,  _Low-dimensional trapped gases_ ,  J. Phys. IV (Proceedings) 116, 5 (2004),  10.1051/jp4:2004116001.
  * [49] B. Rauer, P. Grišins, I. E. Mazets, T. Schweigler, W. Rohringer, R. Geiger, T. Langen and J. Schmiedmayer,  _Cooling of a one-dimensional Bose gas_ ,  Phys. Rev. Lett. 116, 30402 (2016),  10.1103/physrevlett.116.030402.
  * [50] C. Mora and Y. Castin,  _Extension of Bogoliubov theory to quasicondensates_ ,  Phys. Rev. A 67, 053615 (2003),  10.1103/PhysRevA.67.053615.
  * [51] L. Salasnich, A. Parola and L. Reatto,  _Effective wave equations for the dynamics of cigar-shaped and disk-shaped Bose condensates_ ,  Phys. Rev. A 65, 043614 (2002),  10.1103/PhysRevA.65.043614.
  * [52] M. Le Bellac,  _Thermal field theory_ ,  Cambridge Monographs on Mathematical Physics. Cambridge University Press,  10.1017/CBO9780511721700 (1996).
  * [53] M. Pigneur and J. Schmiedmayer,  _Analytical pendulum model for a bosonic josephson junction_ ,  Phys. Rev. A 98, 63632 (2018),  10.1103/physreva.98.063632.
  * [54] A. Steffens, M. Friesdorf, T. Langen, B. Rauer, T. Schweigler, R. Hübener, J. Schmiedmayer, C. A. Riofrio and J. Eisert,  _Towards experimental quantum field tomography with ultracold atoms_ ,  Nature Comm. 6, 7663 (2015),  10.1038/ncomms8663.
  * [55] R. Geiger, T. Langen, I. E. Mazets and J. Schmiedmayer,  _Local relaxation and light-cone-like propagation of correlations in a trapped one-dimensional Bose gas_ ,  New J. Phys. 16, 053034 (2014),  10.1088/1367-2630/16/5/053034.
  * [56] M. Tajik, B. Rauer, T. Schweigler, F. Cataldini, J. Sabino, F. S. Møller, S.-C. Ji, I. E. Mazets and J. Schmiedmayer,  _Designing arbitrary one-dimensional potentials on an atom chip_ ,  Optics Exp. 27, 33474 (2019),  10.1364/OE.27.033474.
  * [57] M. Gluza, J. Sabino, N. H. Y. Ng, G. Vitagliano, M. Pezzutto, Y. Omar, I. Mazets, M. Huber, J. Schmiedmayer and J. Eisert,  _Quantum field thermal machines_ ,  PRX Quantum 2, 030310 (2021),  10.1103/PRXQuantum.2.030310.
  * [58] R. Geiger, T. Langen, I. E. Mazets and J. Schmiedmayer,  _Local relaxation and light-cone-like propagation of correlations in a trapped one-dimensional bose gas_ ,  New J. Phys. 16, 053034 (2014),  10.1088/1367-2630/16/5/053034.
  * [59] P. Calabrese and J. Cardy,  _Time dependence of correlation functions following a quantum quench_ ,  Phys. Rev. Lett. 96, 136801 (2006),  10.1103/PhysRevLett.96.136801.
  * [60] S. Trotzky, Y.-A. Chen, A. Flesch, I. P. McCulloch, U. Schollwöck, J. Eisert and I. Bloch,  _Probing the relaxation towards equilibrium in an isolated strongly correlated one-dimensional Bose gas_ ,  Nature Phys. 8, 325 (2012),  10.1038/nphys2232.
  * [61] J. I. Cirac and P. Zoller,  _Goals and opportunities in quantum simulation_ ,  Nature Phys. 8, 264 (2012),  10.1038/nphys2275.
  * [62] F. Cataldini, F. Møller, M. Tajik, J. Sabino, T. Schweigler, S.-C. Ji, B. Rauer and J. Schmiedmayer,  _Emergent Pauli blocking in a weakly interacting Bose gas_ ,  arXiv preprint, arXiv:2111.13647 (2021).
  * [63] F. D. M. Haldane,  _Effective harmonic-fluid approach to low-energy properties of one-dimensional quantum fluids_ ,  Phys. Rev. Lett. 47, 1840 (1981),  10.1103/PhysRevLett.47.1840.
  * [64] M. Gluza, P. Moosavi and S. Sotiriadis,  _Breaking of Huygens–Fresnel principle in inhomogeneous Tomonaga–Luttinger liquids_ ,  J. Phys. A: Math. Theor. 55(5), 054002 (2022),  10.1088/1751-8121/ac39cc.
  * [65] M. A. Cazalilla, R. Citro, T. Giamarchi, E. Orignac and M. Rigol,  _One dimensional bosons: From condensed matter systems to ultracold gases_ ,  Rev. Mod. Phys. 83, 1405 (2011),  10.1103/RevModPhys.83.1405.