# Condensation and thermalization of an easy-plane ferromagnet in a spinor
Bose gas

Maximilian Prüfer<EMAIL_ADDRESS>Kirchhoff-Institut für Physik,
Universität Heidelberg, Im Neuenheimer Feld 227, 69120 Heidelberg, Germany
Vienna Center for Quantum Science and Technology, Atominstitut, TU Wien,
Stadionallee 2, 1020 Wien, Austria  Daniel Spitz  Institut für Theoretische
Physik, Universität Heidelberg, Philosophenweg 16, 69120 Heidelberg, Germany
Stefan Lannig  Kirchhoff-Institut für Physik, Universität Heidelberg, Im
Neuenheimer Feld 227, 69120 Heidelberg, Germany  Helmut Strobel  Kirchhoff-
Institut für Physik, Universität Heidelberg, Im Neuenheimer Feld 227, 69120
Heidelberg, Germany  Jürgen Berges  Institut für Theoretische Physik,
Universität Heidelberg, Philosophenweg 16, 69120 Heidelberg, Germany  Markus
K. Oberthaler  Kirchhoff-Institut für Physik, Universität Heidelberg, Im
Neuenheimer Feld 227, 69120 Heidelberg, Germany

###### Abstract

The extensive control of spin makes spintronics a promising candidate for
future scalable quantum devices zutic_spintronics_2004 . For the generation of
spin-superfluid systems sonin_spin_2010 , a detailed understanding of the
build-up of coherence and relaxation is necessary. However, to determine the
relevant parameters for robust coherence properties and faithfully witnessing
thermalization, the direct access to space- and time-resolved spin observables
is needed. Here, we study the thermalization of an easy-plane ferromagnet
employing a homogeneous one-dimensional spinor Bose gas stamper-
kurn_spinor_2013 ; sadler_spontaneous_2006 . Building on the pristine control
of preparation and readout kunkel_simultaneous_2019  we demonstrate the
dynamic emergence of long-range coherence glauber_quantum_1963  for the spin
field and verify spin-superfluidity by experimentally testing Landau’s
criterion landau_theory_1941 . We reveal the structure of the emergent quasi-
particles: one ‘massive´ (Higgs) mode, and two ‘massless´ (Goldstone) modes -
a consequence of explicit and spontaneous symmetry breaking, respectively. Our
experiments allow for the first time to observe the thermalization of an easy-
plane ferromagnetic Bose gas; we find agreement for the relevant momentum-
resolved observables with a thermal prediction obtained from an underlying
microscopic model within the Bogoliubov approximation uchino_bogoliubov_2010 ;
phuc_effects_2011 ; symes_static_2014 . Our methods and results pave the way
towards a quantitative understanding of condensation dynamics in large
magnetic spin systems and the study of the role of entanglement and
topological excitations for its thermalization.

Figure 1: Homogeneous spinor Bose gas and easy-plane ferromagnetic properties.
a, We realize a homogeneous spinor BEC of 87Rb in a box-like trapping
potential by a combination of an elongated attractive potential (red) and two
repulsive end caps (green; see Methods for details). The total density (grey
shading) is flat over the extent of the cloud. b, Level structure of the $F=1$
hyperfine manifold. We control the offset energy between the $m$-states by
microwave dressing (blue shading) such that the system features easy-plane
ferromagnetic properties in its ground state. c, The spatial degree of freedom
is continuous, however, in the analysis discretized by the finite pixel size
of the camera and the imaging resolution ($\approx 1.2\,\mu$m). Each imaging
volume (boxes) contains $\approx 500$ atoms which are described by continuous
fields for density and spin. The spins orient themselves in the (easy-)plane
orthogonal to the external magnetic field B. d, The transversal spin features
two different types of excitations: A Goldstone mode and a Higgs mode related
to the excitation of the orientation $\phi_{\text{L}}$ and length
$|F_{\perp}|$, respectively. e, Histogram of the local spin normalized by the
atom number, combining all spatial points and experimental realizations. In
every realization the phase of the central spatial point is subtracted. The
dashed line indicates $|F_{\perp}|=0.75$.

In recent years analog quantum simulators with ultracold atoms allowed for
unprecedented insights by implementing building blocks of complex condensed
matter systems bloch_quantum_2012 ; gross_quantum_2017 . This opens up new
possibilities for studying pressing questions concerning quantum many-body
dynamics and thermalization rigol_thermalization_2008 ; eisert_quantum_2015 ;
reimann_typical_2016 ; gogolin_equilibration_2016 ; mori_thermalization_2018 ;
kaufman_quantum_2016 ; neill_ergodic_2016 ; clos_time-resolved_2016 ;
evrard_many-body_2021 . For probing these phenomena in macroscopic systems
often either the timescales are too short or the control to extract
information is not given such that direct observation of the dynamical
processes is not possible.

Many dynamical phenomena emerging in the many-body limit, such as the build-up
of long-range coherence, superfluidity or spontaneous symmetry breaking, can
be studied in Bose-Einstein condensates (BEC) griffin_bose-einstein_1995 .
Here, the macroscopic occupation of the ground state together with a
spontaneously broken symmetry manifests itself in a globally well-defined
phase of the complex-valued order parameter in each realization. This phase
can be probed experimentally by interferometric measurements, which has been
demonstrated with different platforms bloch_measurement_2000 ;
deng_spatial_2007 ; damm_first-order_2017 . In an easy-plane ferromagnetic
system the order parameter is characterized by a well-defined magnitude in the
transversal plane whereas all orientations in the plane are equally likely.
Theoretically, this is due to a spatial anisotropy, breaking the full
rotational SO(3) symmetric part of the Hamiltonian down to a transversal SO(2)
symmetry. In condensed matter physics prototype models include the XXZ model
proukakis_bose-einstein_2017-2 , which recently has also been realized with
ultracold atoms in lattice systems jepsen_spin_2020  and Rydberg atoms
geier2021floquet ; scholl_microwave-engineering_2022 .

Here, we realize a spinor BEC of 87Rb with easy-plane ferromagnetic properties
stamper-kurn_spinor_2013 ; uchino_spinor_2015  in a quasi-one-dimensional box
trap navon_quantum_2021  (see Fig. 1a). It consists of three internal states,
labelled by their magnetic quantum number $m\in\\{0,\pm 1\\}$. The system
features rotationally invariant ferromagnetic spin-spin interactions described
by $\hat{\mathcal{H}}_{{s}}=c_{1}\int\text{d}V\,\hat{\boldsymbol{F}}^{2}/2$,
where $c_{1}<0$ is the spin-spin interaction constant and
$\hat{\boldsymbol{F}}$ denotes the spin operator (see Methods for details). A
quadratic Zeeman shift $q$ induced by the magnetic field plays the role of the
isotropy-breaking term; it shifts the energy of the $m=\pm 1$ levels (see Fig.
1b) and is explicitly given by
${\hat{\mathcal{H}}_{q}=q\int\text{d}V\,(\hat{N}_{+1}+\hat{N}_{-1})}$. We
adjust $q$ by using off-resonant microwave dressing gerbier_resonant_2006
such that the mean-field ground-state exhibits easy-plane ferromagnetic
properties ($0<q<2n|c_{1}|$; $n$ is the atomic density) and our initial
conditions restrict the dynamics to the spatially averaged longitudinal ($z$-)
spin being zero. In addition to the spin interactions, our system exhibits
SO(3)-invariant density-density interactions described by
$\hat{\mathcal{H}}_{{d}}=c_{0}\int\text{d}V\,\hat{N}^{2}/2$, with interaction
constant $c_{0}$ and $|c_{0}/c_{1}|\approx 200$ stamper-kurn_spinor_2013 .

The capability to extract the relevant order-parameter field
kunkel_simultaneous_2019  allows us to study the build-up of long-range
coherence in a time- and space-resolved fashion; accessing the full structure
factor of the observables defining the Hamiltonian is the handle to faithfully
witness thermalization. We experimentally examine the order parameter, which
is the transversal spin degree of freedom, by acquiring many realizations of
the complex-valued field $F_{\perp}(y)=F_{x}(y)+\text{i}F_{y}(y)$ using
spatially resolved joint measurements based on positive operator valued
measures (POVM) kunkel_simultaneous_2019 ; prufer_experimental_2020 . We
obtain a value for the transversal spin
$F_{\perp}(y)=|F_{\perp}|e^{-i\phi_{L}}$ with length $|F_{\perp}|$ and
orientation in the plane $\phi_{L}$. The position $y$ along the long axis of
the cloud is discretized by our imaging resolution; in each typical imaging
volume we infer the spin from an average over $\sim 500$ atoms which are
described by a spin field, i.e. taking nearly continuous values (see Fig. 1b),
which we identify as the macroscopic order-parameter field describing the spin
condensation.

For studying the condensation dynamics, we initialize the system far from
equilibrium without well-defined spin length, and fluctuations solely in the
plane. We visualize the emergence of a spin ($F_{\perp}$) field by evaluating
the histogram of $F_{\perp}$ taking into account all spatial positions and
realizations (see Fig. 1e). After 5 s, which corresponds to $\approx 10\times$
the typical time scale of the spin interaction energy $t_{s}=h/(n|c_{1}|)$,
the spin is still far from equilibrium and shows large fluctuations in
orientation and length. After 30 s ($\approx 60\times t_{s}$) of evolution
time we find that the fluctuations settle around a well-defined spin length
$|F_{\perp}|$ and the phase $\phi_{\text{L}}$ becomes well-defined over the
whole sample, i.e. long-range order emerges. This is expected for a thermal
state incorporating spontaneous symmetry breaking in the transversal spin
degree of freedom and can be intuitively grasped by looking at the underlying
mexican-hat-like free-energy potential (see Fig. 1d and kawaguchi_spinor_2012
).

Figure 2: Emergence of long-range coherence and superfluidity. a, Absolute
value of first order coherence $|g_{1}(x,y=0)|$ of transversal spin
$F_{\perp}$; reference point $(y=0)$ is chosen at the left edge of the cloud
with system size $L=74\,\mu\text{m}$. We observe a build-up of long-range
order, i.e. for long times the system features non-zero coherence over its
whole size. (Inset) Two-dimensional coherence function $|g_{1}(x,y)|$ after 27
s evolution time. For long times we find the correlations to be translation-
invariant. b, Second order coherence of the transversal spin showing the
evolution and character of spin length fluctuations. c, Superfluid properties
of the spin condensate. Standard deviation along the cloud of spin length
(purple) and density (grey) for different speeds $v$ of the local
perturbation. The rapid increase at finite speed indicates superfluid
properties of spin and density. Insets show representative single realizations
of the spin length and total density in the different regimes.  Figure 3:
Local spin and density control enables probing of quasi-particle properties.
a, Schematics of the three Bogoliubov dispersion relations. b-d, Time
evolution of orientation (b,), total density (c,) and spin length (d,) after
local perturbation of the spin condensate. The upper panels show all evolution
times and lower panels selected 1D cuts. We find a splitting of the imprinted
wavepacket for phase and total density according to the expected linear
dispersion relations (green and blue); solid lines are Gaussian fits.
Strikingly, the speed of sound differs by one order of magnitude reflecting
the energy scales (solid black line corresponds to $v=110\mu\text{m}/\text{s}$
and dashed line to $v=1100\mu\text{m}/\text{s}$). In contrast, the spin length
excitation disperses. e, Results of Gaussian fits to spin length excitation.
$1\sigma$ width (circles), the amplitude $a$ (diamonds) and the integral
$\propto a\sigma$ (squares) are shown. We find an increasing (decreasing)
width (amplitude) while the integral stays nearly constant; this is in
accordance with an underlying gapped quadratic dispersion. f, Oscillation of
the $m=0$ population after perturbing the $k=0$ mode of the spin length. The
oscillation frequency is a measure of the gap $\Delta$ of the quadratic mode
identified in d. All shown error bars are 1 s.d. of the mean. Figure 4:
Structure factors of different observables at late times. We show experimental
structure factors (green diamonds; error bars smaller than plot marker) after
30 s evolution time where the system behaves stationary. Experimental
uncertainties are smaller than marker size. We compare to a thermal prediction
within number-conserving Bogoliubov theory (green line; grey band indicates
68% confidence interval of statistical and systematic uncertainties; for
details on the fit and parameters see Methods). Grey diamonds indicate
reference noise level of a coherent spin state prepared from a single
component gas by a global spin rotation.

To test for eventual spin condensation, we characterize the coherence
properties of the transversal spin by evaluating first and second order
coherence functions glauber_quantum_1963  with
$g_{1}(x,y)\propto{\braket{{\hat{F}}_{\perp}^{\dagger}(x){\hat{F}}_{\perp}(y)}}$
and
${g_{2}(x,y)\propto{\braket{{\hat{F}}_{\perp}^{\dagger}(x){\hat{F}}^{\dagger}_{\perp}(y){\hat{F}}_{\perp}(x){\hat{F}}_{\perp}(y)}}}$,
respectively (see Methods for details). In contrast to earlier experiments
observing the emergence of long-range coherence in one-component BECs
bloch_measurement_2000 ; kohl_growth_2002 ; ritter_observing_2007 , we do not
rely on spatial interference as we access the relevant spin field directly by
joint measurements entailing interference in the internal degrees of freedom
kunkel_simultaneous_2019 ; kunkel_detecting_2022 . We find that coherence is
built up dynamically and the system finally features long-range order, i.e.
non-zero $|g_{1}(x,y)|$ over the whole extent of the atomic cloud (see Fig.
2a). At the same time the spin length fluctuations, quantified by $g_{2}(x,y)$
(see Fig. 2b), settle close to unity at zero distance as expected for a weakly
interacting Bose-Einstein condensate naraschewski_spatial_1999 ;
ottl_correlations_2005 ; hodgman_direct_2011 ; perrin_hanbury_2012 ;
cayla_hanbury-brown_2020 .

To characterize the final state, we first test for superfluidity of the spin
as well as the density; in the spirit of Landau landau_theory_1941 , we drag a
well-localised obstacle (see Methods) coupling to density and spin through the
BEC raman_evidence_1999 ; weimer_critical_2015 ; kim_observation_2020  and
measure the response of the system. We quantify the response by evaluating the
mean standard deviation of the total density and the transversal spin length
along the cloud. The breakdown of superfluidity is signalled by a rapid
increase of the response at a non-zero critical velocity. We find two
different critical velocities for spin and density (see Fig. 2c). While the
spin shows superfluidity up to $v_{c,s}\simeq 3\times 10^{-2}\,$mm/s, the
density tolerates a moving barrier for up to 10 times faster speeds. This is
consistent with the interaction strengths and the corresponding stiffness of
the degrees of freedom.

In the following we address the underlying structure in more detail. With two
spontaneously broken symmetries, the U(1) symmetry of the total density and
the SO(2) symmetry of the spin orientation, we anticipate two Goldstone-like
modes with linear dispersions in the infrared. The different energy scales of
density and spin interactions are reflected in two associated sound speeds;
they are theoretically expected to differ by more than an order of magnitude
which is consistent with the observed critical velocities. Additionally, the
symmetry explicitly broken by $\hat{\mathcal{H}}_{s}+\hat{\mathcal{H}}_{q}$
leads to a Higgs-like gapped mode (see Fig. 3a). Compared to two-component
BECs we find an additional mode due to the increased number of degrees of
freedom recati_breaking_2019 ; cominotti_observation_2021 .

Experimentally, we probe the three different modes by applying local
perturbations (see Fig. 3 and Methods for details). After the perturbation we
observe and analyze the temporal evolution to learn about the underlying
structure of the dispersion relations. First, we probe the linear mode
associated with the spin orientation by imprinting a spatially varying
orientation $\phi_{L}$ pattern onto the thermalized state. The probing scheme
is based on our capabilities to combine global and local radio frequency spin
rotations with fixed relative phases. The initially imprinted Gaussian
wavepacket splits up into two wavepackets travelling with velocities of $\pm
v_{{s}}$ – a clear indication for a linear dispersion relation. To access the
density degree of freedom we imprint a Gaussian-shaped density reduction of
$\sim 5\%$ and observe again a splitting of the initially prepared wave
packet. The difference in energy scales of density and spin is reflected in
the two observed sound speeds $v_{{d}}$ and $v_{{s}}$ which we find to be
$v_{{d}}/v_{{s}}=10\approx\sqrt{|c_{0}/c_{1}|}$. Finally, the gapped mode is
associated with excitations of the spin length which we perturb with a
Gaussian length modulation. Strikingly, we find no splitting but a decaying
amplitude and growing width of the prepared perturbation. Additionally we
measure the gap of this mode; it manifests as a finite oscillation frequency
when exciting the $k=0$ mode (see Fig. 3f). We find a $q$-dependent, non-zero
oscillation frequency consistent with expectations concerning the nature of
the underlying Bogoliubov mode in the easy-plane ferromagnetic phase (see
Ext.Data.Fig.1 and Methods for details).

We now turn to the question if we experimentally realized a thermal ensemble;
to ease its realization we prepare an elongated spin initial condition and
evolve it for 30 s. Experimentally, we extensively characterize the
thermalized state of our system by using different structure factors, such as
those of the spin in transversal as well as in longitudinal ($z$-) direction
and the densities of the three $m$-components. To separate the single
components we employ a Stern-Gerlach magnetic field gradient and a short time-
of-flight (2 ms). We are able to extract momentum-resolved structure factors
by Fourier-transforming the spatial profiles. Special care has to be taken to
measure the total density fluctuations since phase fluctuations are
transformed into density fluctuations during any time-of-flight betz_two-
point_2011 . Therefore we employ an in-situ imaging of the total density to
access the low level of fluctuations which is in our regime of the order of
atom shot-noise (see Methods for details).

The structure factor (see Fig. 4) as well as the local fluctuations (see Ext.
Data Fig. 2) of all observables are consistently described using a thermal
prediction. The latter is obtained for the spinor Bose gas with contact
interactions within the Bogoliubov approximation bogolubov_theory_nodate ;
uchino_bogoliubov_2010  using a single temperature for all three quasiparticle
modes. We thus conclude that the system has evolved to a thermal state within
experimentally accessible timescales. The found temperature is $(57\pm 3)\,$Hz
($\sim 3\,$nK) and thus a factor of $\approx 5\times$ smaller than the
density-density interaction ($nc_{0}=252\pm 54\,$Hz) and more than one order
of magnitude larger than the spin-spin interaction energy scale
($|nc_{1}|=1.17\pm 0.25\,$Hz). The temperature is consistent with theoretical
estimates for the critical temperature for the emergence of the easy-plane
ferromagnetic phase phuc_effects_2011 ; kawaguchi_finite-temperature_2012 .
Interestingly, the single $m$-densities feature high fluctuations also in
comparison with three independent single-component condensates at
corresponding density, interaction and temperature.

We repeat our measurement close to the phase boundary at $q=0$ and find
structures beyond thermal Bogoliubov theory in that case (see Ext.Data Fig.
3). The observed enhanced fluctuations can be associated with long-lived non-
linear excitations of the spin blakie_solitons_2022  which are energetically
less suppressed for lower $q$.

In conclusion, the high degree of control allows us to experimentally observe
the thermalization process of an easy-plane ferromagnet. This sets the
foundations for studies in quantum field settings addressing the microscopic
processes for thermalization as well as its absence due to e.g. long-lived
topological defects. The robust generation of a spin superfluid is a
prerequisite for spin Josephson junctions where finite temperature effects and
spin-density separation can now be studied on a new quantitative level due to
the direct access to the order parameter.

Acknowledgements
We thank Jan Dreher for experimental assistance.
This work is supported by ERC Advanced Grant Horizon 2020 EntangleGen
(Project-ID 694561), the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) under Germany’s Excellence Strategy EXC2181/1-390900948 (the
Heidelberg STRUCTURES Excellence Cluster) and the Collaborative Research
Center - Project-ID 27381115 - SFB 1225 ISOQUANT. M.P. has received funding
from the European Union’s Horizon 2020 research and innovation programme under
the Marie Skłodowska-Curie grant agreement No 101032523.
Author contributions, M.P., S.L. and H.S. took the measurement data. M.P.,
D.S., S.L., H.S. and M.K.O. discussed the measurement results and analysed the
data. D.S., M.P., and J.B. elaborated the theoretical framework. All authors
contributed to the discussion of the results and the writing of the
manuscript.

Competing financial interests
The authors declare no competing financial interests.

Data availability statement
Source data and all other data that support the plots within this paper and
other findings of this study are available from the corresponding author upon
reasonable request.

Methods

##  .1 Experimental details

We prepare a spinor Bose-Einstein condensate of 87Rb in a quasi-one-
dimensional trapping geometry. Details concerning the preparation and readout
of the transversal spin can be found in prufer_observation_2018 ;
prufer_experimental_2020 ; kunkel_simultaneous_2019 ; kunkel_detecting_2022

Here, we employ a box-like trapping potential. We use a weakly focused red-
detuned laser beam creating to a quasi-one-dimensional trapping potential with
$\omega_{l}\approx 2\pi\times 1.7\,$Hz and $\omega_{r}\approx 2\pi\times
170\,$Hz; repulsive potential walls are created by two blue-detuned laser
beams which results in a trapping volume of adjustable size around the centre
of the harmonic trap. The longitudinal harmonic potential is in good
approximation constant over the employed sizes and, thus, effectively leads to
a 1D box-like confinement for the atomic cloud. For the measurements of the
thermalized state shown in Fig. 4 we utilize a box size of $\sim 100\,\mu$m.

Initial conditions: For the detailed observation of the emergence of coherence
we prepare the atoms in the state $\ket{F,m}=\ket{1,0}$, the so-called polar
state. To allow for thermalization for shorter times, we initially prepare a
coherent spin state with maximal length. For this we apply a $\pi/2$-rf
rotation with the atoms initially prepared in the state
$\ket{F,m}=\ket{1,-1}$. As a reference noise level (grey diamonds in Fig. 4
and Ext. Data Fig. 3) for the thermalized structure factor, we prepare this
coherent spin state by performing the rotation after holding the atoms in
$\ket{1,-1}$ for 30 s.

Readout: After the evolution time $t$ we image the atomic densities using
spatially resolved absorption imaging. Employing a Stern-Gerlach magnetic
field gradient followed by a short time-of-flight (TOF; $2\,$ms), we are able
to image the atomic densities of all 8 magnetic sublevels of the $F=1$ and
$F=2$ hyperfine manifolds of the electronic ground state. Additional coherent
microwave and radio-frequency manipulations before the imaging allow us to map
the two spin projections, $F_{x}$ and $F_{y}$, of the transversal spin
prufer_experimental_2020  onto measureable densities. With this we infer the
complex-valued transversal spin
$F_{\perp}(y)=F_{x}(y)+iF_{y}(y)=|F_{\perp}|(y)\text{e}^{-i\phi_{\text{L}}(y)}$
as a function of position $y$. The position $y$ is the centre of a spatial bin
which contains $\sim 500$ atoms and has a spatial extension of $\sim
1.2\,\mu$m along the cloud (we bin three adjacent camera pixels where each
pixel corresponds to $420\,$nm in the atom plane).

For the measurement of the density fluctuations $|N_{\text{tot}}|^{2}(k)$ we
take in-situ images without spin resolution (without Stern-Gerlach
separation). This is important since any free propagation will transform phase
fluctuation to density fluctuations leading to strongly enhanced structure
factor manz_two-point_2010 . It is important to note that the observed
increased fluctuations compared to the spin coherent state by a factor of two
can be a result of only one particle per k-mode. For the spin observables and
the single densities we checked that the enhanced fluctuations due to the TOF
are negligible.

##  .2 Local perturbation

To access the superfluid properties of the spin and density degrees of freedom
we use a localized perturbation (r.m.s width $\sim 5\,\mu$m) that we drag
through the thermalized system. Specifically we use a blue detuned, steerable
laser beam (760 nm) which position is controlled by an acusto-optical
deflector. Using a linear frequency ramp we implement a sweep over the cloud
with fixed velocity which we change over two orders of magnitude. The density
is probed after 35 s and the spin after 20 s evolution time. The ramp duration
for the lowest speed is $\approx 18\,$s.

Local perturbation of the Larmor phase: We use a combination of global and
local rf rotations (see lannig_collisions_2020  for details on local rf
rotations). A first global $\pi/2$-rf rotation around the $x$-axis maps the
$z$-axis onto the $y$-axis. Using a local rotation with a well-defined phase
with respect to the global rotation we perform a rotation with variable angle
around the $y$-axis. Because of the performed mapping this effectively leads
to rotation around the $z$-axis in the original coordinate system. At time
$\Delta\tau=210\,\mu$s after the first global rf pulse we apply a global rf
$\pi$-pulse followed by a second global rf $\pi$/2-pulse after another time
delay of $\Delta\tau$, where all pulses rotate around the same axis. This
constitutes a spin echo sequence which additionally executes a full $2\pi$
spin rotation which ensures that the global rotation pulses do not excite the
system. The last $\pi$/2-pulse maps the local rotation axis back to the z-axis
in the original system. The perturbation has an approximate Gaussian shape
with a r. m. s. width of $\sim 5\,\mu$m according to the shape of the used
laser beam.

Local perturbation of the total density: We reduce the total density locally
by $\approx 5\%$ by shining a blue detuned laser beam ($760\,$nm) onto the
centre of the cloud. We adiabatically ramp up the potential in 100 ms such
that we get no further excitations in the density; after the ramp the
potential is instantaneously switched off to generate the wavepacket.

Local perturbation of the transversal spin length $|F_{\perp}|$: We induce a
local density reduction by applying the same blue-detuned laser beam. During
the evolution time of 30 s we let the system thermalize subject to the local
density reduction. This effectively leads to a spatially dependent mean-field
ground state spin length. We linearly ramp down the potential over 50 ms; this
implements an adiabatic ramp for the total density and a rapid switch off for
the spin.

For experimentally accessing the gap we excite the $k=0$ mode of the spin
length by changing the phase of the $m=0$ component (spinor phase) globally.
For this we use two microwave $\pi$-pulses between $\ket{1,0}$ and
$\ket{2,0}$, where the second pulse is phase shifted by $\Delta\phi$. We
record the subsequent oscillations of the $m=0$ population and fit a
sinusoidal function to extract the frequency. The theoretical prediction for
the gap $\Delta$, deduced from the oscillation, and the $m=0$ ground state
population $n_{0}$ in the easy-plane ferromagnetic phase is
uchino_bogoliubov_2010 :

$\Delta=\sqrt{4n^{2}c_{1}^{2}-q^{2}}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\text{and}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,n_{0}=\frac{1}{2}-\frac{q}{4nc_{1}}\,.$
(1)

Assuming $n_{+1}=n_{-1}$ these formulae hold also true for $0>q>2nc_{1}$
(dashed lines in Ext. Data Fig. 1).

##  .3 Experimental coherence functions and structure factor

In every experimental realization $(i)$ we measure atomic densities from which
we infer single shot realizations $O^{(i)}$ of different observables
$\hat{O}$. The quantum expectation value is approximated by averaging over
many realizations as

$O=\braket{\hat{O}}=\frac{1}{N_{S}}\sum_{i=1}^{N_{S}}O^{(i)}\,,$ (2)

where $N_{S}$ is the number of realizations.

The coherence functions of the transversal spin are explicitly given by:

$g_{1}(x,y)=\frac{\braket{\hat{F}_{\perp}^{\dagger}(x)\hat{F}_{\perp}(y)}}{\sqrt{\braket{{\hat{F}}_{\perp}^{\dagger}(x){\hat{F}}_{\perp}(x)}\braket{{\hat{F}}_{\perp}^{\dagger}(y){\hat{F}}_{\perp}(y)}}}$
(3)

and

$g_{2}(x,y)=\frac{{\braket{{\hat{F}}_{\perp}^{\dagger}(x){\hat{F}}^{\dagger}_{\perp}(y){\hat{F}}_{\perp}(x){\hat{F}}_{\perp}(y)}}}{{\braket{{\hat{F}}_{\perp}^{\dagger}(x){\hat{F}}_{\perp}(x)}\braket{{\hat{F}}_{\perp}^{\dagger}(y){\hat{F}}_{\perp}(y)}}}\,.$
(4)

For the inferred single shot results of the transversal spin $F_{\perp}(x)$
the † is treated as the complex conjugate.

The structure factors as a function of the spatial momentum $k$ are defined as

$\braket{}{\hat{O}(k)}{{}^{2}}=|O(k)|^{2}=\frac{1}{N_{\text{tot}}}\frac{1}{N_{S}}\sum_{i=1}^{N_{S}}|DFT_{x\rightarrow
k}\left(O^{(i)}(x)-{O(x)}\right)|^{2}\,,$ (5)

where $DFT_{x\rightarrow k}$ is the discrete Fourier transform, $k=1/\lambda$
the spatial momentum. All structure factors are normalized by the mean total
atom number $N_{\text{tot}}$ to obtain an atom number independent measure for
the fluctuations and allow comparison between theory and experiment; for the
total density structure factor a value of one corresponds to the atomic shot
noise level.

##  .4 Bogoliubov transformations in the easy-plane ferromagnetic phase

We explicitely derive the Bogoliubov transformations in the easy-plane
ferromagnetic phase ($0<q/(n|c_{1}|)<2$). Here we set $\hbar=1$.

In terms of total density and spin operators

$\hat{N}(\mathbf{x})=\sum_{m=-1}^{1}\hat{N}_{m}(\mathbf{x})=\sum_{m=-1}^{1}\hat{\psi}^{\dagger}_{m}(\mathbf{x})\hat{\psi}_{m}(\mathbf{x}),\qquad\hat{F}_{\nu}(\mathbf{x})=\sum_{m,m^{\prime}=-1}^{1}\hat{\psi}_{m}^{\dagger}(\mathbf{x})(f^{\nu})_{mm^{\prime}}\hat{\psi}_{m^{\prime}}(\mathbf{x}),$
(6)

with the spin-1 matrices

$f^{x}=\frac{1}{\sqrt{2}}\left(\begin{matrix}0&1&0\\\ 1&0&1\\\
0&1&0\end{matrix}\right),\qquad
f^{y}=\frac{i}{\sqrt{2}}\left(\begin{matrix}0&-1&0\\\ 1&0&-1\\\
0&1&0\end{matrix}\right),\qquad f^{z}=\left(\begin{matrix}1&0&0\\\ 0&0&0\\\
0&0&-1\end{matrix}\right),$ (7)

the system Hamiltonian reads

$\hat{\mathcal{H}}=\int
d^{3}\mathbf{x}\bigg{[}\sum_{m=-1}^{1}\hat{\psi}_{m}^{\dagger}(\mathbf{x})\bigg{(}-\frac{\nabla^{2}}{2M}+qm^{2}\bigg{)}\hat{\psi}_{m}(\mathbf{x})+\frac{c_{0}}{2}:\hat{N}^{2}(\mathbf{x}):+\frac{c_{1}}{2}\sum_{\nu=x,y,z}:\hat{F}_{\nu}^{2}(\mathbf{x}):\bigg{]}.$
(8)

With momentum-space creation and annihilation operators

$\hat{a}_{\mathbf{k},m}^{\dagger}=\frac{1}{\sqrt{V}}\int
d^{3}\mathbf{x}\,\hat{\psi}^{\dagger}_{m}(\mathbf{x})e^{i\mathbf{k}\mathbf{x}},\qquad\hat{a}_{\mathbf{k},m}=\frac{1}{\sqrt{V}}\int
d^{3}\mathbf{x}\,\hat{\psi}_{m}(\mathbf{x})e^{-i\mathbf{k}\mathbf{x}},$ (9)

the Hamiltonian in the number-conserving Bogoliubov approximation becomes
kawaguchi_spinor_2012

$\displaystyle\hat{\mathcal{H}}_{B}=$
$\displaystyle\,E_{0}+\sum_{\mathbf{k}\neq\mathbf{0},m}(\epsilon_{\mathbf{k}}+qm^{2}-\mu)\hat{n}_{\mathbf{k},m}$
$\displaystyle\,+\frac{N}{V}\sum_{j,j^{\prime},m,m^{\prime}}\sum_{\mathbf{k}\neq\mathbf{0}}(\Gamma_{jj^{\prime},m^{\prime}m}+\Gamma_{jm,m^{\prime}j^{\prime}})\zeta_{j^{\prime}}\zeta_{m^{\prime}}^{*}\hat{a}^{\dagger}_{\mathbf{k},j}\hat{a}_{\mathbf{k},m}$
$\displaystyle\,+\frac{N}{2V}\sum_{j,j^{\prime},m,m^{\prime}}\sum_{\mathbf{k}\neq\mathbf{0}}\Gamma_{jj^{\prime},mm^{\prime}}\big{(}\zeta^{*}_{j}\zeta^{*}_{m}\hat{a}_{-\mathbf{k},m^{\prime}}\hat{a}_{\mathbf{k},j^{\prime}}+\zeta_{m^{\prime}}\zeta_{j^{\prime}}\hat{a}^{\dagger}_{\mathbf{k},j}\hat{a}^{\dagger}_{-\mathbf{k},m}\big{)}.$
(10)

with $\epsilon_{\mathbf{k}}=\mathbf{k}^{2}/(2M)$, atom mass $M$, total atom
number $N$ and
$\hat{n}_{\mathbf{k},m}=\hat{a}^{\dagger}_{\mathbf{k},m}\hat{a}_{\mathbf{k},m}$.
Theoretically, momentum $|\mathbf{k}|$ corresponds to wavelength
$\lambda=2\pi/|\mathbf{k}|$. The spinor $(\zeta_{m})$ specifies the normalized
condensate configuration and we set
$\hat{a}_{\mathbf{k}=0,m}=\sqrt{N}\zeta_{m}$.
$\Gamma_{jj^{\prime},mm^{\prime}}$ denotes density and spin interactions,

$\Gamma_{jj^{\prime},mm^{\prime}}\equiv
c_{0}\delta_{jj^{\prime}}\delta_{mm^{\prime}}+c_{1}\sum_{\nu=x,y,z}f^{\nu}_{jj^{\prime}}f^{\nu}_{mm^{\prime}}.$
(11)

The ground state energy is given by

$E_{0}\equiv
N\bigg{[}\sum_{m}qm^{2}|\zeta_{m}|^{2}+\frac{N-1}{2V}\sum_{j,j^{\prime},m,m^{\prime}}\Gamma_{jj^{\prime},mm^{\prime}}\zeta_{j}^{*}\zeta_{m}^{*}\zeta_{m^{\prime}}\zeta_{j^{\prime}}\bigg{]},$
(12)

the chemical potential reads

$\mu\equiv\sum_{m}qm^{2}|\zeta_{m}|^{2}+\frac{2N-1}{2V}\sum_{j,j^{\prime},m,m^{\prime}}\Gamma_{jj^{\prime},mm^{\prime}}\zeta^{*}_{j}\zeta^{*}_{m}\zeta_{m^{\prime}}\zeta_{j^{\prime}}.$
(13)

We set $\sin\theta=\sqrt{1/2-q/(4n|c_{1}|)}$, such that the mean-field ground
state reads $\zeta=(\sin\theta/\sqrt{2},\cos\theta,\sin\theta/\sqrt{2})$
uchino_bogoliubov_2010 . The initial orthogonal transformation, given by

$\left(\begin{matrix}\hat{a}_{\mathbf{k},d}\\\ \hat{a}_{\mathbf{k},\theta}\\\
\hat{a}_{\mathbf{k},f_{z}}\end{matrix}\right)=\left(\begin{matrix}\sin\theta/\sqrt{2}&\cos\theta&\sin\theta/\sqrt{2}\\\
\cos\theta/\sqrt{2}&-\sin\theta&\cos\theta/\sqrt{2}\\\
1/\sqrt{2}&0&-1/\sqrt{2}\end{matrix}\right)\left(\begin{matrix}\hat{a}_{\mathbf{k},1}\\\
\hat{a}_{\mathbf{k},0}\\\ \hat{a}_{\mathbf{k},-1}\end{matrix}\right)\equiv
A(\theta)\left(\begin{matrix}\hat{a}_{\mathbf{k},1}\\\
\hat{a}_{\mathbf{k},0}\\\ \hat{a}_{\mathbf{k},-1}\end{matrix}\right),$ (14)

leads to a description of the system in terms of longitudinal and transversal
spin fluctuations.

The longitudinal ($z$-) spin fluctuations can be diagonalized using the
Bogoliubov transformation uchino_bogoliubov_2010

$\hat{b}_{\mathbf{k},f_{z}}=u_{\mathbf{k},f_{z}}\hat{a}_{\mathbf{k},f_{z}}+v_{\mathbf{k},f_{z}}\hat{a}^{\dagger}_{-\mathbf{k},f_{z}},$
(15)

where

$u_{\mathbf{k},f_{z}}\equiv\sqrt{\frac{\epsilon_{\mathbf{k}}+q/2+E_{\mathbf{k},f_{z}}}{2E_{\mathbf{k},f_{z}}}},\qquad
v_{\mathbf{k},f_{z}}\equiv\sqrt{\frac{\epsilon_{\mathbf{k}}+q/2-E_{\mathbf{k},f_{z}}}{2E_{\mathbf{k},f_{z}}}},$
(16)

with the dispersion

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

To diagonalize transversal spin fluctuations we follow the diagonalization
procedure outlined in kawaguchi_spinor_2012 . We obtain mode energies $\pm
E_{\mathbf{k},+}$ and $\pm E_{\mathbf{k},-}$ as in kawaguchi_spinor_2012 ,
explicitely given by

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

with

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

Defining

$h_{00}\equiv n(c_{0}+c_{1}-c_{1}),\quad h_{01}\equiv q\sin(2\theta)/2,\quad
h_{11}\equiv-2nc_{1}+nc_{q},\quad h_{211}\equiv nc_{q}$ (20)

and

$\displaystyle u_{\mathbf{k},\pm,1}=$ $\displaystyle\;-\big{(}h_{01}\big{[}\pm
2E_{1}+4\epsilon_{\mathbf{k}}^{2}+2\epsilon_{\mathbf{k}}(2E_{\mathbf{k},\pm}+h_{00}+2h_{11}-h_{211})$
$\displaystyle\;\qquad+(h_{11}-h_{211})(2E_{\mathbf{k},\pm}+h_{11}+h_{211})\big{]}\big{)}\big{/}\big{(}4\epsilon_{\mathbf{k}}h_{01}^{2}+2\epsilon_{\mathbf{k}}(h_{11}-h_{00})h_{211}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+h_{211}(\pm
2E_{1}+h_{11}^{2}-h_{211}^{2})\big{)},$ (21) $\displaystyle
u_{\mathbf{k},\pm,2}=$
$\displaystyle\;\big{(}4\epsilon_{\mathbf{k}}^{2}(h_{00}-h_{11})-2(E_{\mathbf{k},\pm}+h_{11})(\pm
2E_{1}+h_{11}^{2}-h_{211}^{2})$
$\displaystyle\;\qquad-2\epsilon_{\mathbf{k}}(\pm
2E_{1}+4h_{01}^{2}-2E_{\mathbf{k},\pm}(h_{00}-h_{11})-2h_{00}h_{11}+3h_{11}^{2}$
$\displaystyle\;\qquad-
h_{211}^{2})\big{)}\big{/}\big{(}8\epsilon_{\mathbf{k}}h_{01}^{2}+4\epsilon_{\mathbf{k}}(h_{11}-h_{00})h_{211}+2h_{211}(\pm
2E_{1}+h_{11}^{2}-h_{211}^{2})\big{)},$ (22) $\displaystyle
v_{\mathbf{k},\pm,1}=$ $\displaystyle\;\frac{h_{01}\big{[}\pm
2E_{1}+2\epsilon_{\mathbf{k}}(h_{00}+h_{211})+(h_{11}-h_{211})(2E_{\mathbf{k},\pm}+h_{11}+h_{211})\big{]}}{4\epsilon_{\mathbf{k}}h_{01}^{2}+2\epsilon_{\mathbf{k}}(h_{11}-h_{00})h_{211}+h_{211}(\pm
2E_{1}+h_{11}^{2}-h_{211}^{2})},$ (23) $\displaystyle v_{\mathbf{k},\pm,2}=$
$\displaystyle\;1,$ (24)

together with normalization factors

$\mathcal{N}_{\mathbf{k},\pm}\equiv\sqrt{u_{\mathbf{k},\pm,1}^{2}+u_{\mathbf{k},\pm,2}^{2}-v_{\mathbf{k},\pm,1}^{2}-v_{\mathbf{k},\pm,2}^{2}}$
(25)

we find Bogoliubov transformation matrices (in the parametrization of
uchino_bogoliubov_2010 )

$U_{\mathbf{k},d\theta}=\left(\begin{matrix}u_{\mathbf{k},+,1}/\mathcal{N}_{\mathbf{k},+}&u_{\mathbf{k},+,2}/\mathcal{N}_{\mathbf{k},+}\\\
u_{\mathbf{k},-,1}/\mathcal{N}_{\mathbf{k},-}&u_{\mathbf{k},-,2}/\mathcal{N}_{\mathbf{k},-}\end{matrix}\right),\qquad
V_{\mathbf{k},d\theta}=\left(\begin{matrix}-v_{\mathbf{k},+,1}/\mathcal{N}_{\mathbf{k},+}&-v_{\mathbf{k},+,2}/\mathcal{N}_{\mathbf{k},+}\\\
-v_{\mathbf{k},-,1}/\mathcal{N}_{\mathbf{k},-}&-v_{\mathbf{k},-,2}/\mathcal{N}_{\mathbf{k},-}\end{matrix}\right).$
(26)

These fulfil the identities

$U_{\mathbf{k},d\theta}U_{\mathbf{k},d\theta}^{\dagger}-V_{\mathbf{k},d\theta}V_{\mathbf{k},d\theta}^{\dagger}=1,\qquad
U_{\mathbf{k},d\theta}^{*}V_{\mathbf{k},d\theta}^{\dagger}-V_{\mathbf{k},d\theta}^{*}U_{\mathbf{k},d\theta}^{\dagger}=0,$
(27)

as required for the transformations to preserve canonical commutation
relations. This requirement is not fulfilled for the transformations given in
uchino_bogoliubov_2010 .

The complete transformation matrices diagonalizing the Bogoliubov Hamiltonian
(10) read

$U_{\mathbf{k}}=\left(\begin{matrix}U_{\mathbf{k},d\theta}&0\\\
0&u_{\mathbf{k},f_{z}}\end{matrix}\right)A(\theta),\qquad
V_{\mathbf{k}}=\left(\begin{matrix}V_{\mathbf{k},d\theta}&0\\\
0&v_{\mathbf{k},f_{z}}\end{matrix}\right)A(\theta).$ (28)

##  .5 Thermal structure factors from Bogoliubov theory

We provide analytical computations of thermal structure factors in the spinor
Bose gas from Bogoliubov theory. We are interested in correlators of the form

$\langle\hat{C}^{\dagger}(\mathbf{x})\hat{C}(\mathbf{y})\rangle_{\beta,s}$
(29)

for a composite field $\hat{C}(\mathbf{x})$ given by

$\hat{C}(\mathbf{x})=\sum_{m,m^{\prime}=-1}^{+1}\hat{\psi}_{m}^{\dagger}(\mathbf{x})c_{mm^{\prime}}\hat{\psi}_{m^{\prime}}(\mathbf{x})$
(30)

with $c_{mm^{\prime}}$ a $3\times 3$ matrix corresponding to the type of
spectrum under investigation; $c=f^{x}+if^{y}$ leads to the transversal
magnetization spectrum, $c=f^{z}$ describes the spectrum of magnetization in
$z$-direction, $c=\textrm{diag}(1,1,1)$ describes the total density spectrum.
In Eq. (29) $\langle\cdot\rangle_{\beta,s}$ indicates the thermal expectation
value at inverse temperature $\beta=1/(k_{B}T)$ with symmetrically (Weyl-)
ordered arguments. We compare with symmetrically ordered predictions since
expectation values of experimental observables are inferred from realizations
of observables $O^{(i)}$ given by polynomials of complex numbers (cf.
symes_static_2014  for a similar normal-ordered computation). Fourier-
transforming Eq. (29) with respect to the relative coordinate
$\mathbf{x}-\mathbf{y}$, we obtain the structure factor

$\displaystyle\langle\hat{C}^{\dagger}(\mathbf{k})\hat{C}(\mathbf{k})\rangle_{\beta,s}=$
$\displaystyle\;\int
d(\mathbf{x}-\mathbf{y})\,\langle\hat{C}^{\dagger}(\mathbf{x})\hat{C}(\mathbf{y})\rangle_{\beta,s}\,e^{-i\mathbf{k}(\mathbf{x}-\mathbf{y})}$
(31) $\displaystyle=$
$\displaystyle\;\frac{1}{V}\sum_{m,m^{\prime},n,n^{\prime}}c^{\dagger}_{mm^{\prime}}c_{nn^{\prime}}\sum_{\mathbf{p},\mathbf{q}}\langle\hat{a}^{\dagger}_{\mathbf{p}+\mathbf{k},m}\hat{a}_{\mathbf{p},m^{\prime}}\hat{a}^{\dagger}_{\mathbf{q},n}\hat{a}_{\mathbf{q}+\mathbf{k},n^{\prime}}\rangle_{\beta,s}.$
(32)

To total density structure factors a photon shot noise level of 0.6 after
normalization is added.

In the Bogoliubov approximation Eq. (32) simplifies as follows. We replace
zero modes of creation and annihilation operators by numbers,
$\hat{a}_{\mathbf{0},m}=\sqrt{N_{\mathrm{cond}}}\zeta_{m}$,
$N_{\mathrm{cond}}$ the total number of condensate atoms. Contributions from
fluctuating modes $\hat{a}_{\mathbf{k}\neq\mathbf{0},m}$ are computed via
Wick’s theorem. Symmetrically ordered propagators are defined as

$\displaystyle G^{11}_{mm^{\prime}}(\mathbf{k})\equiv$
$\displaystyle\;\langle\hat{a}_{\mathbf{k},m}\hat{a}^{\dagger}_{\mathbf{k},m^{\prime}}\rangle_{\beta,s},\qquad
G^{22}_{mm^{\prime}}(\mathbf{k})\equiv\langle\hat{a}^{\dagger}_{\mathbf{k},m}\hat{a}_{\mathbf{k},m^{\prime}}\rangle_{\beta,s},$
(33) $\displaystyle G^{12}_{mm^{\prime}}(\mathbf{k})\equiv$
$\displaystyle\;\langle\hat{a}_{\mathbf{k},m}\hat{a}_{-\mathbf{k},m^{\prime}}\rangle_{\beta,s},\qquad
G^{21}_{mm^{\prime}}(\mathbf{k})\equiv\langle\hat{a}^{\dagger}_{\mathbf{k},m}\hat{a}^{\dagger}_{-\mathbf{k},m^{\prime}}\rangle_{\beta,s}.$
(34)

The first two of these we refer to as normal propagators, the second two as
anomalous propagators. Any normal propagator evaluated for non-diagonal
momenta such as
$\langle\hat{a}^{\dagger}_{\mathbf{k},m}\hat{a}_{\mathbf{k}^{\prime},m^{\prime}}\rangle_{\beta,s}$
for $\mathbf{k}^{\prime}\neq\mathbf{k}$ and any anomalous propagator evaluated
for non-anti-diagonal momenta such as
$\langle\hat{a}_{\mathbf{k},m}\hat{a}_{-\mathbf{k}^{\prime},m^{\prime}}\rangle_{\beta,s}$
for $\mathbf{k}^{\prime}\neq\mathbf{k}$ equates to zero. We then find for
$\mathbf{k}=\mathbf{0}$,

$\displaystyle\langle\hat{C}^{\dagger}(\mathbf{0})\hat{C}(\mathbf{0})\rangle_{\beta,s}=$
$\displaystyle\;\frac{N_{\mathrm{cond}}(N_{\mathrm{cond}}-1)}{V}\sum_{m,m^{\prime},n,n^{\prime}}c^{\dagger}_{mm^{\prime}}c_{nn^{\prime}}\zeta_{m}^{*}\zeta_{n}^{*}\zeta_{m^{\prime}}\zeta_{n^{\prime}}+\mathcal{O}(N_{\mathrm{cond}}),$
(35)

and for non-zero modes $\mathbf{k}\neq\mathbf{0}$,

$\displaystyle\langle\hat{C}^{\dagger}(\mathbf{k})\hat{C}(\mathbf{k})\rangle_{\beta,s}=$
$\displaystyle\;\frac{N_{\mathrm{cond}}}{V}\sum_{m,m^{\prime},n,n^{\prime}}c^{\dagger}_{mm^{\prime}}c_{nn^{\prime}}\bigg{[}\zeta_{m}^{*}\zeta_{n}^{*}G^{12}_{m^{\prime}n^{\prime}}(\mathbf{k})+\zeta_{m}^{*}\zeta_{n^{\prime}}G^{11}_{m^{\prime}n}(\mathbf{k})$
$\displaystyle\qquad\qquad\qquad+\zeta_{n}^{*}\zeta_{m^{\prime}}G^{22}_{mn^{\prime}}(\mathbf{k})+\zeta_{m^{\prime}}\zeta_{n^{\prime}}G^{21}_{mn}(\mathbf{k})\bigg{]}+\mathcal{O}(1),$
(36)

where we have used that the propagators $G^{ab}_{mm^{\prime}}(\mathbf{k})$
only depend on the absolute value of the momentum, $|\mathbf{k}|$.
Experimental structure factors, normalized analogously, are compared to
$\langle\hat{C}^{\dagger}(\mathbf{k})\hat{C}(\mathbf{k})\rangle_{\beta,s}/n_{\mathrm{cond}}$,
$n_{\mathrm{cond}}$ the total condensate atom density. The photon shot noise
of the absorption imaging is determined to be $0.6$ after normalization and is
added on the thermal prediction of the total density fluctuations.

In the easy-plane ferromagnetic phase thermal propagators
$G^{ab}_{mm^{\prime}}(\mathbf{k}\neq\mathbf{0})$ can be computed from the
Bogoliubov transformations derived in Methods Sec. .4. Bogoliubov
quasiparticle excitations $f_{z}$ and $\pm$ are expressed in terms of
fundamental magnetic sublevel excitations as

$\left(\begin{matrix}\hat{b}_{\mathbf{k},j}\\\
\hat{b}^{\dagger}_{-\mathbf{k},j}\end{matrix}\right)=\sum_{m=-1}^{+1}\left(\begin{matrix}U_{\mathbf{k},jm}&V_{\mathbf{k},jm}\\\
V_{-\mathbf{k},jm}^{*}&U_{-\mathbf{k},jm}^{*}\end{matrix}\right)\left(\begin{matrix}\hat{a}_{\mathbf{k},m}\\\
\hat{a}_{-\mathbf{k},m}^{\dagger}\end{matrix}\right),$ (37)

for $j\in\\{f_{z},\pm\\}$. Inverting this, we obtain kawaguchi_spinor_2012

$\left(\begin{matrix}\hat{a}_{\mathbf{k},m}\\\
\hat{a}_{-\mathbf{k},m}^{\dagger}\end{matrix}\right)=\sum_{j\in\\{f_{z},\pm\\}}\left(\begin{matrix}U_{\mathbf{k},mj}^{\dagger}&-V_{\mathbf{k},mj}^{T}\\\
-V_{-\mathbf{k},mj}^{\dagger}&U_{-\mathbf{k},mj}^{T}\end{matrix}\right)\left(\begin{matrix}\hat{b}_{\mathbf{k},j}\\\
\hat{b}^{\dagger}_{-\mathbf{k},j}\end{matrix}\right).$ (38)

The propagators can now efficiently be computed from the tensor product

$\displaystyle\left(\begin{matrix}G^{12}_{mm^{\prime}}(\mathbf{k})&G^{11}_{mm^{\prime}}(\mathbf{k})\\\
G^{22}_{mm^{\prime}}(-\mathbf{k})&G^{21}_{mm^{\prime}}(-\mathbf{k})\end{matrix}\right)=\langle\left(\begin{matrix}\hat{a}_{\mathbf{k},m}\\\
\hat{a}_{-\mathbf{k},m}^{\dagger}\end{matrix}\right)\otimes\left(\hat{a}_{-\mathbf{k},m^{\prime}},\hat{a}^{\dagger}_{\mathbf{k},m^{\prime}}\right)\rangle_{\beta,s}$
(39)
$\displaystyle=\sum_{j,j^{\prime}\in\\{f_{z},\pm\\}}\left(\begin{matrix}U_{\mathbf{k},mj}^{\dagger}&-V_{\mathbf{k},mj}^{T}\\\
-V_{-\mathbf{k},mj}^{\dagger}&U_{-\mathbf{k},mj}^{T}\end{matrix}\right)\left(\begin{matrix}\langle\hat{b}_{\mathbf{k},j}\hat{b}_{-\mathbf{k},j^{\prime}}\rangle_{\beta,s}&\langle\hat{b}_{\mathbf{k},j}\hat{b}^{\dagger}_{\mathbf{k},j^{\prime}}\rangle_{\beta,s}\\\
\langle\hat{b}^{\dagger}_{-\mathbf{k},j}\hat{b}_{-\mathbf{k},j^{\prime}}\rangle_{\beta,s}&\langle\hat{b}^{\dagger}_{-\mathbf{k},j}\hat{b}^{\dagger}_{\mathbf{k},j^{\prime}}\rangle_{\beta,s}\end{matrix}\right)\left(\begin{matrix}U_{-\mathbf{k},j^{\prime}m^{\prime}}^{*}&-V_{-\mathbf{k},j^{\prime}m^{\prime}}^{*}\\\
-V_{\mathbf{k},j^{\prime}m^{\prime}}&U_{\mathbf{k},j^{\prime}m^{\prime}}\end{matrix}\right).$
(40)

The Bogoliubov quasiparticle modes $\hat{b}_{\mathbf{k},j}$ are occupied
thermally,

$\langle\hat{b}^{\dagger}_{\mathbf{k},j}\hat{b}_{\mathbf{k},j^{\prime}}\rangle_{\beta}=\delta_{jj^{\prime}}n_{\beta}(E_{\mathbf{k},j}),\qquad\langle\hat{b}_{\mathbf{k},j}\hat{b}^{\dagger}_{\mathbf{k},j^{\prime}}\rangle_{\beta}=\delta_{jj^{\prime}}(n_{\beta}(E_{\mathbf{k},j})+1),$
(41)

with the Bose-Einstein distribution $n_{\beta}(E_{\mathbf{k},j})\equiv
1/(\exp(\beta E_{\mathbf{k},j})-1)$. Anomalous propagators of
$\hat{b}_{\mathbf{k},j}$-modes are zero. Insertion of (41) into (40) and using
that $U_{\mathbf{k},mj}$ and $V_{\mathbf{k},mj}$ only depend on $|\mathbf{k}|$
leads to

$\displaystyle G^{11}_{mm^{\prime}}(\mathbf{k})=$
$\displaystyle\;\sum_{j\in\\{f_{z},\pm\\}}\bigg{[}U^{\dagger}_{\mathbf{k},mj}\left(n_{\beta}(E_{\mathbf{k},j})+\frac{1}{2}\right)U_{\mathbf{k},jm^{\prime}}+V^{T}_{\mathbf{k},mj}\left(n_{\beta}(E_{\mathbf{k},j})+\frac{1}{2}\right)V^{*}_{\mathbf{k},jm^{\prime}}\bigg{]},$
(42) $\displaystyle G^{22}_{mm^{\prime}}(\mathbf{k})=$
$\displaystyle\;\sum_{j\in\\{f_{z},\pm\\}}\bigg{[}V^{\dagger}_{\mathbf{k},mj}\left(n_{\beta}(E_{\mathbf{k},j})+\frac{1}{2}\right)V_{\mathbf{k},jm^{\prime}}+U^{T}_{\mathbf{k},mj}\left(n_{\beta}(E_{\mathbf{k},j})+\frac{1}{2}\right)U^{*}_{\mathbf{k},jm^{\prime}}\bigg{]},$
(43) $\displaystyle G^{12}_{mm^{\prime}}(\mathbf{k})=$
$\displaystyle\;-\sum_{j\in\\{f_{z},\pm\\}}\bigg{[}U^{\dagger}_{\mathbf{k},mj}\left(n_{\beta}(E_{\mathbf{k},j})+\frac{1}{2}\right)V_{\mathbf{k},jm^{\prime}}+V^{T}_{\mathbf{k},mj}\left(n_{\beta}(E_{\mathbf{k},j})+\frac{1}{2}\right)U^{*}_{\mathbf{k},jm^{\prime}}\bigg{]},$
(44) $\displaystyle G^{21}_{mm^{\prime}}(\mathbf{k})=$
$\displaystyle\;-\sum_{j\in\\{f_{z},\pm\\}}\bigg{[}V^{\dagger}_{\mathbf{k},mj}\left(n_{\beta}(E_{\mathbf{k},j})+\frac{1}{2}\right)U_{\mathbf{k},jm^{\prime}}+U^{T}_{\mathbf{k},mj}\left(n_{\beta}(E_{\mathbf{k},j})+\frac{1}{2}\right)V^{*}_{\mathbf{k},jm^{\prime}}\bigg{]}.$
(45)

With these expressions thermal structure factors can be readily computed from
Eq. (36).

##  .6 Fitting thermal Bogoliubov theory structure factors

Using a least-squares fitting procedure and Gibbs sampling, systematic as well
as statistical uncertainties on the optimal set of parameters are estimated.
Given experimental structure factors
$S_{\hat{C},\exp}(k)=\langle\hat{C}^{\dagger}(k)\hat{C}(k)\rangle$ with
$\hat{C}\in\mathcal{S}:=\\{\hat{N}_{+1},\hat{N}_{0},\hat{N}_{-1},\hat{F}_{z},\hat{F}_{\perp}\\}$,
we determine an optimal set of parameters $T,q,nc_{1}$ by minimizing

$\chi^{2}(T,q,nc_{1};k_{\max})=\sum_{\hat{C}\in\mathcal{S}}\sum_{k}^{k_{\max}}\frac{(S_{\hat{C},\exp}(k)-S_{\hat{C},\mathrm{Bog}}(k;T,q,nc_{1}))^{2}}{\Delta
S_{\hat{C},\exp}(k)^{2}},$ (46)

with
$S_{\hat{C},\mathrm{Bog}}(k;T,q,nc_{1})=\langle\hat{C}^{\dagger}(k)\hat{C}(k)\rangle_{1/(k_{B}T),s}$
the Bogoliubov theory structure factor computed for parameters $T,q,nc_{1}$,
and $\Delta S_{\hat{C},\exp}(k)$ the standard deviation of
$S_{\hat{C},\exp}(k)$ computed from experimental realizations. Technical
correlations of the absorption imaging are described by real-space signals
convoluted with a Gaussian of r. m. s. width $w=5.0\,\mu\mathrm{m}$, taken
into account by the multiplication of momentum-space structure factors with a
Gaussian of width $2\pi/w$ kunkel_spatially_2018 . Throughout the fitting
procedure we set $c_{0}/c_{1}\simeq-216$ in accordance with stamper-
kurn_spinor_2013 . In the definition of $\chi^{2}$ we did not include the
structure factor of the total density.

We define a distribution of parameters for specific $k_{\max}$,

$W(T,q,nc_{1};k_{\max})\sim\exp\left(-\chi^{2}(T,q,nc_{1};k_{\max})/2\right).$
(47)

We exploit Gibbs sampling to draw $i=1,\dots,100$ approximately i.i.d. samples
$(T^{(i)}(k_{\max})$, $q^{(i)}(k_{\max})$, $nc_{1}^{(i)}(k_{\max}))$ from
$W(T,q,nc_{1};k_{\max})$, which only requires corresponding conditional
distributions normalized individually. For each $k_{\max}$ we compute their
mean
$(\overline{T}(k_{\max}),\overline{q}(k_{\max}),\overline{nc_{1}}(k_{\max}))$.
We repeat this for 5 values of $k_{\max}$ evenly spaced between $0.1\cdot
2\pi/\mu\mathrm{m}$ and $0.2\cdot 2\pi/\mu\mathrm{m}$. Collecting all $5\cdot
100$ samples in a single array, we take the mean values
$\overline{T},\overline{q},\overline{nc_{1}}$ of all samples as final
parameter estimates and their distances to the boundaries of 68% confidence
intervals as corresponding error estimates. Errors include systematic fit
uncertainties. We obtain the final fit parameters

$T=(57.4\pm 2.9)\,\mathrm{Hz},\quad q=(0.30\pm 0.08)\,\mathrm{Hz},\quad
nc_{1}=(-1.17\pm 0.25)\,\mathrm{Hz},$ (48)

such that $nc_{0}=(252\pm 54)\,\mathrm{Hz}$ and $q/(nc_{1})=(-0.26\pm 0.09)$.

##  .7 Drawing Bogoliubov theory samples

Bogoliubov theory being quadratic in fluctuating field creation and
annihilation operators, it is fully described by zero modes and a suitable
covariance matrix of fluctuations. The latter can be constructed from the
propagators $G^{ij}_{mm^{\prime}}(\mathbf{k})$.

Given a one-dimensional real-space lattice
$\\{-\mathcal{N},\dots,\mathcal{N}\\}\cdot a$ with lattice spacing
$a=L/(2\mathcal{N})$, the corresponding momentum-space lattice reads
$\\{-\mathcal{N},\dots,\mathcal{N}\\}\cdot\pi/(\mathcal{N}a)$. With $\Delta
z=z_{1}-z_{2}\in\\{-2\mathcal{N},\dots,2\mathcal{N}\\}$ for
$z_{i}\in\\{-\mathcal{N},\dots,\mathcal{N}\\}$ we compute real-space
propagators via

${\tilde{G}^{ij}_{mm^{\prime}}(\Delta
z)=\frac{1}{L}\sum_{p=-\mathcal{N}}^{\mathcal{N}}G^{ij}_{mm^{\prime}}(p)\exp\left(\frac{2\pi
ip\Delta z}{2\mathcal{N}+1}\right).}$ (49)

We assemble these into magnetic sublevel-specific covariance matrices,

$\displaystyle\;\textrm{Cov}_{mm^{\prime}}=\langle\left(\begin{matrix}\hat{\psi}_{m}(-\mathcal{N})\\\
\vdots\\\ \hat{\psi}_{m}(\mathcal{N})\\\
\hat{\psi}_{m}^{\dagger}(-\mathcal{N})\\\ \vdots\\\
\hat{\psi}_{m}^{\dagger}(\mathcal{N})\end{matrix}\right)\left(\begin{matrix}\hat{\psi}_{m^{\prime}}(-\mathcal{N}),\dots,\hat{\psi}_{m^{\prime}}(\mathcal{N}),\hat{\psi}^{\dagger}_{m^{\prime}}(-\mathcal{N}),\dots,\hat{\psi}^{\dagger}_{m^{\prime}}(\mathcal{N})\end{matrix}\right)\rangle_{\beta,s}$
(50)
$\displaystyle\;=\left(\begin{matrix}\tilde{G}^{12}_{mm^{\prime}}(0)&\tilde{G}^{12}_{mm^{\prime}}(-1)&\dots&\tilde{G}^{12}_{mm^{\prime}}(-2\mathcal{N})&\tilde{G}^{11}_{mm^{\prime}}(0)&\tilde{G}^{11}_{mm^{\prime}}(-1)&\dots&\tilde{G}^{11}_{mm^{\prime}}(-2\mathcal{N})\\\
\tilde{G}^{12}_{mm^{\prime}}(1)&\tilde{G}^{12}_{mm^{\prime}}(0)&\dots&\tilde{G}^{12}_{mm^{\prime}}(-2\mathcal{N}+1)&\tilde{G}^{11}_{mm^{\prime}}(1)&\tilde{G}^{11}_{mm^{\prime}}(0)&\dots&\tilde{G}^{11}_{mm^{\prime}}(-2\mathcal{N}+1)\\\
\vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\ddots&\vdots\\\
\tilde{G}^{12}_{mm^{\prime}}(2\mathcal{N})&\tilde{G}^{12}_{mm^{\prime}}(2\mathcal{N}-1)&\dots&\tilde{G}^{12}_{mm^{\prime}}(0)&\tilde{G}^{11}_{mm^{\prime}}(2\mathcal{N})&\tilde{G}^{11}_{mm^{\prime}}(2\mathcal{N}-1)&\dots&\tilde{G}^{11}_{mm^{\prime}}(0)\\\
\tilde{G}^{22}_{mm^{\prime}}(0)&\tilde{G}^{22}_{mm^{\prime}}(-1)&\dots&\tilde{G}^{22}_{mm^{\prime}}(-2\mathcal{N})&\tilde{G}^{21}_{mm^{\prime}}(0)&\tilde{G}^{21}_{mm^{\prime}}(-1)&\dots&\tilde{G}^{21}_{mm^{\prime}}(-2\mathcal{N})\\\
\tilde{G}^{22}_{mm^{\prime}}(1)&\tilde{G}^{22}_{mm^{\prime}}(0)&\dots&\tilde{G}^{22}_{mm^{\prime}}(-2\mathcal{N}+1)&\tilde{G}^{21}_{mm^{\prime}}(1)&\tilde{G}^{21}_{mm^{\prime}}(0)&\dots&\tilde{G}^{21}_{mm^{\prime}}(-2\mathcal{N}+1)\\\
\vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\ddots&\vdots\\\
\tilde{G}^{22}_{mm^{\prime}}(2\mathcal{N})&\tilde{G}^{22}_{mm^{\prime}}(2\mathcal{N}-1)&\dots&\tilde{G}^{22}_{mm^{\prime}}(0)&\tilde{G}^{21}_{mm^{\prime}}(2\mathcal{N})&\tilde{G}^{21}_{mm^{\prime}}(2\mathcal{N}-1)&\dots&\tilde{G}^{21}_{mm^{\prime}}(0)\end{matrix}\right),$
(51)

having exploited spatial homogeneity. We decompose complex field operators
into real components,
$\hat{\psi}(x)=\frac{1}{\sqrt{2}}(\hat{\psi}_{1}(x)+i\hat{\psi}_{2}(x))$,
translating into the unitary transformation

$\left(\begin{matrix}\hat{\psi}_{m,1}(-\mathcal{N})\\\ \vdots\\\
\hat{\psi}_{m,1}(\mathcal{N})\\\ \hat{\psi}_{m,2}(-\mathcal{N})\\\ \vdots\\\
\hat{\psi}_{m,2}(\mathcal{N})\end{matrix}\right)=A\left(\begin{matrix}\hat{\psi}_{m}(-\mathcal{N})\\\
\vdots\\\ \hat{\psi}_{m}(\mathcal{N})\\\
\hat{\psi}_{m}^{\dagger}(-\mathcal{N})\\\ \vdots\\\
\hat{\psi}_{m}^{\dagger}(\mathcal{N})\end{matrix}\right),\qquad
A=\frac{1}{\sqrt{2}}\left(\begin{matrix}I&I\\\ -iI&iI\end{matrix}\right),$
(52)

with $I$ the $(2\mathcal{N}+1)\times(2\mathcal{N}+1)$-dimensional identity
matrix. We define the final covariance matrix of the theory as

$\textrm{Cov}=\left(\begin{matrix}A\,\textrm{Cov}_{+1,+1}A^{T}&A\,\textrm{Cov}_{+1,0}A^{T}&A\,\textrm{Cov}_{+1,-1}A^{T}\\\
A\,\textrm{Cov}_{0,+1}A^{T}&A\,\textrm{Cov}_{0,0}A^{T}&A\,\textrm{Cov}_{0,-1}A^{T}\\\
A\,\textrm{Cov}_{-1,+1}A^{T}&A\,\textrm{Cov}_{-1,0}A^{T}&A\,\textrm{Cov}_{-1,-1}A^{T}\end{matrix}\right).$
(53)

Finally, we sample $i=1,\dots,N_{\mathrm{sample}}$ field realizations

$\psi^{(i)}=\left(\begin{matrix}\psi^{(i)}_{+1}\\\ \psi^{(i)}_{0}\\\
\psi^{(i)}_{-1}\end{matrix}\right),\qquad\psi^{(i)}_{m}=\left(\begin{matrix}\psi^{(i)}_{m,1}(-\mathcal{N}),&\cdots,&\psi^{(i)}_{m,1}(\mathcal{N}),&\psi^{(i)}_{m,2}(-\mathcal{N}),&\cdots,&\psi^{(i)}_{m,2}(\mathcal{N})\end{matrix}\right)^{T},$
(54)

from the multivariate Gaussian distribution with zero mean vector and
covariance matrix $\mathrm{Cov}$. This corresponds to samples from the Wigner
distribution of the symmetrically ordered Bogoliubov theory of fluctuating
modes at inverse temperature $\beta$. We sample fields in position space
instead of momentum space, since in momentum space, having decomposed the
operators $\hat{a}_{k,m}$ into
$\hat{a}_{k,m}=(\hat{a}_{k,m,1}+i\hat{a}_{k,m,2})/\sqrt{2}$, the components
$\hat{a}_{k,m,j}$ need to satisfy
$\hat{a}_{k,m,j}^{\dagger}=\hat{a}_{-k,m,j}$, such that samples of individual
momentum modes cannot be drawn independently. From the fluctuating
realizations $\psi^{(i)}$ we can compute realizations of the individual spin
sublevel fields in real space,

$\psi^{(i)}_{m}(x)=\frac{1}{\sqrt{2}}[\psi^{(i)}_{m,1}(x)+i\psi^{(i)}_{m,2}(x)]+\sqrt{n_{\mathrm{cond}}}\zeta_{m},$
(55)

with $n_{\mathrm{cond}}$ the condensate density. We explicitely checked that
for increasing sample numbers structure factors computed from samples
$\psi^{(i)}_{m}(x)$ converge towards their expectations
$G^{ij}_{mm^{\prime}}(k)$.

The composite operator histograms displayed in Ext. Data Fig. 2 are computed
from composite profiles of individual realizations given by
$\sum_{m,m^{\prime}=-1}^{+1}(\psi^{(i)}_{m}(x))^{*}c_{mm^{\prime}}\psi^{(i)}_{m^{\prime}}(x)/\sqrt{n_{\mathrm{cond}}}$
with matrices $c$ as denoted in Methods .5.

## References

  * (1) Žutić, I., Fabian, J. & Das Sarma, S.  Spintronics: Fundamentals and applications.  _Reviews of Modern Physics_ 76, 323–410 (2004).
  * (2) Sonin, E. B.  Spin currents and spin superfluidity.  _Advances in Physics_ 59, 181–255 (2010).
  * (3) Stamper-Kurn, D. M. & Ueda, M.  Spinor Bose gases: Symmetries, magnetism, and quantum dynamics.  _Reviews of Modern Physics_ 85, 1191–1244 (2013).
  * (4) Sadler, L. E., Higbie, J. M., Leslie, S. R., Vengalattore, M. & Stamper-Kurn, D. M.  Spontaneous symmetry breaking in a quenched ferromagnetic spinor Bose–Einstein condensate.  _Nature_ 443, 312–315 (2006).
  * (5) Kunkel, P. _et al._ Simultaneous Readout of Noncommuting Collective Spin Observables beyond the Standard Quantum Limit.  _Physical Review Letters_ 123, 063603 (2019).
  * (6) Glauber, R. J.  The Quantum Theory of Optical Coherence.  _Physical Review_ 130, 2529–2539 (1963).
  * (7) Landau, L.  Theory of the Superfluidity of Helium II.  _Physical Review_ 60, 356–358 (1941).
  * (8) Uchino, S., Kobayashi, M. & Ueda, M.  Bogoliubov theory and Lee-Huang-Yang corrections in spin-1 and spin-2 Bose-Einstein condensates in the presence of the quadratic Zeeman effect.  _Physical Review A_ 81, 063632 (2010).
  * (9) Phuc, N. T., Kawaguchi, Y. & Ueda, M.  Effects of thermal and quantum fluctuations on the phase diagram of a spin-1 87 Rb Bose-Einstein condensate.  _Physical Review A_ 84, 043645 (2011).
  * (10) Symes, L. M., Baillie, D. & Blakie, P. B.  Static structure factors for a spin-1 Bose-Einstein condensate.  _Physical Review A_ 89, 053628 (2014).
  * (11) Bloch, I., Dalibard, J. & Nascimbène, S.  Quantum simulations with ultracold quantum gases.  _Nature Physics_ 8, 267–276 (2012).
  * (12) Gross, C. & Bloch, I.  Quantum simulations with ultracold atoms in optical lattices.  _Science_ 357, 995–1001 (2017).
  * (13) Rigol, M., Dunjko, V. & Olshanii, M.  Thermalization and its mechanism for generic isolated quantum systems.  _Nature_ 452, 854–858 (2008).
  * (14) Eisert, J., Friesdorf, M. & Gogolin, C.  Quantum many-body systems out of equilibrium.  _Nature Physics_ 11, 124–130 (2015).
  * (15) Reimann, P.  Typical fast thermalization processes in closed many-body systems.  _Nature Communications_ 7, 10821 (2016).
  * (16) Gogolin, C. & Eisert, J.  Equilibration, thermalisation, and the emergence of statistical mechanics in closed quantum systems.  _Reports on Progress in Physics_ 79, 056001 (2016).
  * (17) Mori, T., Ikeda, T. N., Kaminishi, E. & Ueda, M.  Thermalization and prethermalization in isolated quantum systems: a theoretical overview.  _Journal of Physics B: Atomic, Molecular and Optical Physics_ 51, 112001 (2018).
  * (18) Kaufman, A. M. _et al._ Quantum thermalization through entanglement in an isolated many-body system.  _Science_ 353, 794–800 (2016).
  * (19) Neill, C. _et al._ Ergodic dynamics and thermalization in an isolated quantum system.  _Nature Physics_ 12, 1037–1041 (2016).
  * (20) Clos, G., Porras, D., Warring, U. & Schaetz, T.  Time-Resolved Observation of Thermalization in an Isolated Quantum System.  _Physical Review Letters_ 117, 170401 (2016).
  * (21) Evrard, B., Qu, A., Dalibard, J. & Gerbier, F.  From many-body oscillations to thermalization in an isolated spinor gas.  _Physical Review Letters_ 126, 063401 (2021).
  * (22) Griffin, A., Snoke, D. W. & Stringari, S. (eds.) _Bose-Einstein Condensation_ (Cambridge University Press, 1995).
  * (23) Bloch, I., Hänsch, T. W. & Esslinger, T.  Measurement of the spatial coherence of a trapped Bose gas at the phase transition.  _Nature_ 403, 166–170 (2000).
  * (24) Deng, H., Solomon, G. S., Hey, R., Ploog, K. H. & Yamamoto, Y.  Spatial Coherence of a Polariton Condensate.  _Physical Review Letters_ 99, 126403 (2007).
  * (25) Damm, T., Dung, D., Vewinger, F., Weitz, M. & Schmitt, J.  First-order spatial coherence measurements in a thermalized two-dimensional photonic quantum gas.  _Nature Communications_ 8, 158 (2017).
  * (26) Kollath, C., Giamarchi, T. & Rüegg, C.  Bose-Einstein Condensation in Quantum Magnets.  In Proukakis, N. P., Snoke, D. W. & Littlewood, P. B. (eds.) _Universal Themes of Bose-Einstein Condensation_ , 549–568 (Cambridge University Press, 2017).
  * (27) Jepsen, P. N. _et al._ Spin transport in a tunable Heisenberg model realized with ultracold atoms.  _Nature_ 588, 403–407 (2020).
  * (28) Geier, S. _et al._ Floquet hamiltonian engineering of an isolated many-body spin system.  _Science_ 374, 1149–1152 (2021).
  * (29) Scholl, P. _et al._ Microwave-engineering of programmable XXZ Hamiltonians in arrays of Rydberg atoms.  _arXiv:2107.14459_ (2022).
  * (30) Uchino, S.  Spinor Bose gas in an elongated trap.  _Physical Review A_ 91, 033605 (2015).
  * (31) Navon, N., Smith, R. P. & Hadzibabic, Z.  Quantum gases in optical boxes.  _Nature Physics_ 17, 1334–1341 (2021).
  * (32) Gerbier, F., Widera, A., Fölling, S., Mandel, O. & Bloch, I.  Resonant control of spin dynamics in ultracold quantum gases by microwave dressing.  _Physical Review A_ 73, 041602 (2006).
  * (33) Prüfer, M. _et al._ Experimental extraction of the quantum effective action for a non-equilibrium many-body system.  _Nature Physics_ 16, 1012–1016 (2020).
  * (34) Kawaguchi, Y. & Ueda, M.  Spinor Bose-Einstein condensates.  _Physics Reports_ 520, 253–381 (2012).
  * (35) Köhl, M., Davis, M. J., Gardiner, C. W., Hänsch, T. W. & Esslinger, T.  Growth of Bose-Einstein Condensates from Thermal Vapor.  _Physical Review Letters_ 88, 080402 (2002).
  * (36) Ritter, S. _et al._ Observing the Formation of Long-Range Order during Bose-Einstein Condensation.  _Physical Review Letters_ 98, 090402 (2007).
  * (37) Kunkel, P. _et al._ Detecting Entanglement Structure in Continuous Many-Body Quantum Systems.  _Physical Review Letters_ 128, 020402 (2022).
  * (38) Naraschewski, M. & Glauber, R. J.  Spatial coherence and density correlations of trapped Bose gases.  _Physical Review A_ 59, 4595–4607 (1999).
  * (39) Öttl, A., Ritter, S., Köhl, M. & Esslinger, T.  Correlations and Counting Statistics of an Atom Laser.  _Physical Review Letters_ 95, 090404 (2005).
  * (40) Hodgman, S. S., Dall, R. G., Manning, A. G., Baldwin, K. G. H. & Truscott, A. G.  Direct Measurement of Long-Range Third-Order Coherence in Bose-Einstein Condensates.  _Science_ 331, 1046–1049 (2011).
  * (41) Perrin, A. _et al._ Hanbury Brown and Twiss correlations across the Bose–Einstein condensation threshold.  _Nature Physics_ 8, 195–198 (2012).
  * (42) Cayla, H. _et al._ Hanbury-Brown and Twiss bunching of phonons and of the quantum depletion in a strongly-interacting Bose gas.  _Physical Review Letters_ 125, 165301 (2020).
  * (43) Raman, C. _et al._ Evidence for a Critical Velocity in a Bose-Einstein Condensed Gas.  _Physical Review Letters_ 83, 2502–2505 (1999).
  * (44) Weimer, W. _et al._ Critical Velocity in the BEC-BCS Crossover.  _Physical Review Letters_ 114, 095301 (2015).
  * (45) Kim, J. H., Hong, D. & Shin, Y.  Observation of two sound modes in a binary superfluid gas.  _Physical Review A_ 101, 061601 (2020).
  * (46) Recati, A. & Piazza, F.  Breaking of Goldstone modes in a two-component Bose-Einstein condensate.  _Physical Review B_ 99, 064505 (2019).
  * (47) Cominotti, R. _et al._ Observation of Massless and Massive Collective Excitations with Faraday Patterns in a Two-Component Superfluid.  _arXiv:2112.09880_ (2021).
  * (48) Betz, T. _et al._ Two-Point Phase Correlations of a One-Dimensional Bosonic Josephson Junction.  _Physical Review Letters_ 106, 020407 (2011).
  * (49) Bogoliubov, N.  On the theory of superfluidity.  _J. Phys._ 11, 23 (1947).
  * (50) Kawaguchi, Y., Phuc, N. T. & Blakie, P. B.  Finite-temperature phase diagram of a spin-1 Bose gas.  _Physical Review A_ 85, 053611 (2012).
  * (51) Yu, X. & Blakie, P. B.  Propagating ferrodark solitons in a superfluid: Exact solutions and anomalous dynamics.  _Phys. Rev. Lett._ 128, 125301 (2022).
  * (52) Prüfer, M. _et al._ Observation of universal dynamics in a spinor Bose gas far from equilibrium.  _Nature_ 563, 217–220 (2018).
  * (53) Manz, S. _et al._ Two-point density correlations of quasicondensates in free expansion.  _Physical Review A_ 81, 031610 (2010).
  * (54) Lannig, S. _et al._ Collisions of three-component vector solitons in Bose-Einstein condensates.  _Physical Review Letters_ 125, 170401 (2020).
  * (55) Kunkel, P. _et al._ Spatially distributed multipartite entanglement enables EPR steering of atomic clouds.  _Science_ 360, 413–416 (2018).

Extended Data Figure 1: Measurement of the gap by observation of temporal
oscillations of the $k=0$ mode. a, We measure the gap of the quadratic spin
mode by a global rotation of the spinor phase. We record the resulting
oscillations of the fractional $m=0$ population as a function of evolution
time after the rotation. We fit a sinusoidal function (solid line) to infer
the frequency. b, Extracted oscillation frequency (diamonds) and mean value of
the $m=0$ population (circles). We compare to theoretical expectations for the
easy-plane phase (solid lines; see Methods eq. (1)). The dashed line
extrapolates the expectations to $q<0$ under the assumption of equal
populations of $m=\pm 1$. For the theory curves we use $nc_{1}=1.3\,$Hz.
Extended Data Figure 2: Histograms of local observables in the thermalized
state. Histograms obtained from evaluating the local observations of the
experimental data presented in Fig. 4 (green bars). Here, each local
observable is normalized to the square-root of the local mean of the total
atom number. On top we display theoretical estimates from 1000 samples
generated according to thermal Bogoliubov theory with parameters as in Fig. 4
(grey line; grey band indicates 68% confidence interval including statistical
and systematic uncertainties). The mean value of each histogram is subtracted.
For details on the sampling procedure see Methods. Extended Data Figure 3:
Structure factor close to $q=0$. We show experimental power spectra of
different spin and density degrees of freedom close to $q=0$ (green diamonds).
The grey diamonds represent the fluctuations of a coherent spin state with
comparable atom numbers. We compare to thermal Bogoliubov theory predictions
for the same parameters as displayed in Fig. 4 but with $q=0$ (green line;
grey band indicates 68% confidence interval of statistical and systematic
uncertainties). Experimentally, we find that for momenta in the range of
$0.02\mu\text{m}^{-1}$ to $0.1\mu\text{m}^{-1}$ the fluctuations are higher
than for the thermal predictions for all observables (except the transversal
spin $F_{\perp}$). The length scale of these fluctuations is in accordance
with observable localized long-lived non-linear excitations which are not
present in the thermalized data of Fig. 4.