Hanbury Brown-Twiss correlations and multimode dynamics of quenched spinor Bose-Einstein condensates

We have studied the interaction between multiple, competing spatial modes that are excited by a quantum quench of an antiferromagnetic spinor Bose-Einstein condensate. We observed Hanbury Brown-Twiss correlations and associated super-Poissonian noise in the mode populations. The decay of these correlations was consistent with experimentally observed spin domain patterns. Data were compared with a real-space Bogoliubov theory as well as numerical solution of the coupled Gross-Pitaevskii equations that were seeded by quantum noise via the truncated Wigner approximation. The spatial modes that were both observed experimentally and deduced theoretically are intimately connected to the inhomogeneous density profile of the condensate, which imparts many rich features to the dynamics.


I. INTRODUCTION
Quantum phase transitions have become an important arena of exploration using ultracold atoms. In contrast to the traditional view of these transitions that has focused on equilibrium properties [1], quantum gases afford a dynamical view made experimentally possible by rapidly quenching the system [2]. A major advantage exists for such quench experiments over equilibrium studies very close to the critical point, where the timescales become longer and the equilibration condition harder to fulfill. Instead, in a quench scenario, the transition is crossed very rapidly from well outside its boundaries, and if the system has no time to respond to the change, then a full panoply of dynamical behavior can be observed that is characteristic of the symmetries broken by the transition (see Figure 1). This includes i) instability of the initial state and formation of seeds of the ground state on the opposite side of the critical point, ii) rapid expansion of the seeds to macroscopic proportions, and iii) slow growth to steady-state. Using this approach, we have recently been able to make sub-Hz level (picoKelvin) precision measurements on the transition boundary in an antiferromagnetic spinor Bose-Einstein condensate, although the system temperature was 400 nK [3], well above the energy scales of the phase transition itself.
In the present work we extend our earlier studies on the q = 0 quantum phase transition in antiferromagnetic spinor Bose-Einstein condensate (BEC). We present measurements of Hanbury-Brown and Twiss correlations and study their decay length. While our previous works have confirmed the Bogoliubov predictions for the instability rate [3,5], here we provide many details-both experimental and theoretical-that were missing in the earlier works. For instance, we use statistical analysis to elucidate the variety and wealth of modes that were observed * Electronic address: chandra.raman@physics.gatech.edu FIG. 1: (Color Online). Quantum quenches reveal two phases of an antiferromagnetic spinor BEC. a) Quantum phase transition between polar and antiferromagnetic (AF) phases at zero quadratic Zeeman shift q. (b) Dynamical instability induced by a rapid quench from q1 > 0 to q2 < 0 results in the appearance of a small number of AF domains at short times T1 after the quench, and a quasi-equilibrium between the phases at longer times T2. (c) Measurement of total fraction in the AF phase (the mF = ±1 population fraction) versus time showing the various dynamical phases (figure adapted from [4]). Error bars are the standard deviation of 3 separate measurements.
to form in the experiment, and compare with Bogoliubov predictions for these modes. Similar to the quench experiments of reference [6], we have observed that the maximally unstable modes are localized near the center of the cloud, where the density is highest. We show that this result appears naturally from the real space calculation Typeset by REVT E X arXiv:1806.04231v1 [cond-mat.quant-gas] 11 Jun 2018 of the Bogoliubov eigenvectors, which are different from the momentum modes typical of a homogeneous BEC.
The paper is organized as follows. Section II provides background on the phase transition being studied, and describes the experimental method in detail. Section III contains the experimental data on Hanbury-Brown and Twiss correlations. Section IV explains the theoretical technique and section V compares theory with a number of experimental observations. Section VI is a conclusion and outlook.

A. Hamiltonian
Multi-component spinor BECs offer unique opportunities to perform dynamical studies due to a number of quantum phase transitions that can be accessed within the spin sector of the Hamiltonian [7][8][9][10][11][12][13][14][15]. Consider the Hamiltonian for spin F = 1 Bose-Einstein condensates (see [16] and references therein): The first term is the energy of spin-dependent interactions. The interaction coefficient c 2 = 4π 2 3M (a 2 − a 0 ) arises from spin-changing collisions that can convert two m F = 0 atoms into an m F = ±1 pair and vice-versa, a process constrained by the conservation of angular momentum. Here M is the atomic mass, and a 2,0 are the scattering lengths for atom pairs whose total angular momentum F tot = 2 and 0, respectively. These intrinsic interactions are antiferromagnetic for c 2 > 0 or ferromagnetic for c 2 < 0. In this work we consider a sodium BEC, for which c 2 > 0. The overall density profile n(r) is determined by the chemical potential µ, which is larger in magnitude by a factor of 100 compared with Eqn. (1) for typical values of q. Thus, we have removed the spinindependent terms from the Hamiltonian, which are assumed to be a constant.F,F z are the vector spin-1 operator and its z-projection, respectively. Hereafter in this work, we write m ≡ m F . The second term in Eqn. (1) is the quadratic Zeeman shift due to the external magnetic field, q =qB 2 0 + q M . B 0 is the magnetic field at the trap center,q = 276 Hz/Gauss 2 is the coefficient of the quadratic Zeeman shift for sodium atoms, and g F = 1/2 and µ B are the Lande gfactor and Bohr magneton, respectively [17]. In addition to the static magnetic field, we introduce an additional term q M , which is the shift caused by a microwave magnetic field through the AC Zeeman effect on the F = 1, m sublevels [4,5]. The microwave field generates an additional shift q M < 0, allowing us to access the q < 0 region of the phase diagram without having to change the static magnetic field.
For an antiferromagnetic spinor BEC prepared in an initial state with zero net magnetization F z = 0, the ground state of the above Hamiltonian for q > 0 is a polar condensate consisting of a single component-the m = 0 spin projection that minimizes F 2 z . For q < 0 the ground state maximizes the same quantity through a superposition of two components m = ±1, a so-called antiferromagnetic phase [18]. We study the effect of a quantum quench across the zero temperature quantum phase transition at q = 0.

B. Experimental Method
Sodium Bose-Einstein condensates in a single focus optical trap were prepared in the m = 0 state in a static magnetic field. The protocol is described in earlier work [5], but we include relevant details here. The peak density n 0 = 5 × 10 14 cm −3 and axial Thomas-Fermi radius R x = 340 µm were measured to an accuracy of 5%, from which we determined the peak spin-dependent interaction energy c 2 n 0 = h × 120 Hz. The axial and radial trapping frequencies were 7 and 470 Hz, respectively, accurate to 10%. The radial Thomas-Fermi radius was R ⊥ = 5 µm, and thus the aspect ratio of the cigar-shaped cloud was ≈ 70:1. The measured temperature was 400 nK, close to the chemical potential of 360 nK.
Our experiment required precise control over the bias magnetic field. We applied a magnetic field B x aligned with the long axis x of the cigar and tuned the field gradient dB x /dx to cancel ambient field inhomogeneities in our vacuum chamber to within ±10 µG using a procedure that is described in [3]. The microwave field used to control q M was generated by an HP8648B synthesizer with ∼ 100 Hz accuracy. Its frequency was tuned below the "clock" transition, |F = 1, m = 0 → |F = 2, m = 0 at 1.772 GHz [19] by an amount between 260 to 470 kHz, resulting in a negative shift q M < 0 that counteracts the positive shift from the static fieldqB 2 0 . The axial extent of the cloud, R x = 340µm, was substantially larger than the typical domain size of ∼ 30µm. Thus multiple domains can form in this system. Although the system is multi-mode, the spin modes are one-dimensional due to the quench energetics, as can be seen from the following argument. The tight confinement along the transverse dimensions of the cigar implies a large transverse zero point energy for spin excitations. We can estimate these from a two-dimensional circular box model similar to Reference [20]. The single particle eigenfunctions in transverse polar coordinates ρ, φ are ∝ J l (β nl ρ R ⊥ )e ilφ , where J l is the l-th Bessel function with zeros β nl . The box eigenenergies are then nl = 2 β 2 nl /(2M R 2 ⊥ ). The lowest 3 energies are 10 , 11 , 21 whose numerical values are, for our parameters, h× 50,130 and 230 Hz, respectively. For the data collected here very close to the phase transition point, the quadratic Zeeman energy released was h × 8 Hz 10 and thus transverse excitation was impossible. The quench experiment consisted of rapidly switching q from q i > 0 to a final value q f < 0 at t = 0. This was accomplished by sudden turn-on of the microwave field amplitude using a radiofrequency switch with submicrosecond rise time. Following a variable hold time, we switched off the trap and used time-of-flight Stern-Gerlach (TOF-SG) observations [16] to record an absorption image for each experimental run. Due to the highly anisotropic expansion of the 70:1 aspect ratio, cigar-shaped BEC, the x-axis of the image was identical to the axial coordinate within the trap, with one minor difference-each spin component was offset by a fixed amount, ∝ m, that was determined by its Stern-Gerlach separation in time-of-flight. Using this information we could undo the SG separation to extract a onedimensional density distribution in each of the 3 spin components, n i (x); i = 0, ±1, with a spatial resolution of 10 µm [4].
Time-of-flight Stern-Gerlach imaging is a highly sensitive technique, but can be limited by technical, rather than fundamental noise. In our case, technical noise was mostly caused by interference fringes whose position varied from shot to shot. We used powerful post-processing techniques to minimize this effect [21]. This allowed us to observe data at hold times as short as 16 ms after the quench, where the time-of-flight m = ±1 atom signal on the camera corresponded to a mere 9 atoms/pixel 2 . The standard deviation of the resulting image noise could be reduced to a value that was no more than 1.7 times higher than the optical shot noise. The sensitivity could not be further reduced without increasing the fluence of the imaging light to reduce shot noise contributions. However, we found this incurred nonlinear effects and motion induced blurring of the cloud. Future experiments aimed at increasing sensitivity could boost the atom signal at even earlier times by mixing the m = ±1 spins with that of the strong m = 0 cloud through spin rotations, effectively deploying a heterodyne technique.

III. EXPERIMENTAL RESULTS ON HANBURY BROWN-TWISS CORRELATIONS
A. Spin fluctuations Figure 2 encapsulates the multi-mode character of the instability, and reveals several interesting features. The upper panel shows 4 representative Stern-Gerlach images taken on different experimental runs with a hold time of 20 ms after the quench. Spin-exchange collisions are responsible for converting on average 15% of the atoms from m = 0 into m = ±1 pairs, which form localized domains, seen as vertical stripes in the images. Since the pairs were spin-correlated, the location and number of domains was highly correlated between m = 1 and m = −1 atoms, as also observed in earlier work [4]. In the current work we focus on the behavior within each spin cloud m = ±1, where the number and location of domains varied stochastically from run to run. This fluctuation is not due to technical reasons, but is a fundamental feature of the quench, as unoccupied spin excitation modes become rapidly occupied, with a highly variable number and mode distribution.
We analyzed these fluctuations using statistical methods on an ensemble of measurements on identically prepared BECs. Thirty images were collected at each quench hold time, whose mean atom numberN and standard deviation σ N were computed. We suppressed shot to shot atom number fluctuations by filtering only those whose atom number lay inside the rangeN < N <N + σ N , resulting in samples with less than 10% atom number fluctuations. This allowed us to more clearly observe the intrinsic noise.
For this filtered data set, each camera image was reduced to a one-dimensional slice f (x), where x is the axial coordinate within the image. To accomplish this, we converted the absorption images to atomic column density (atoms/pixel 2 ) through the known absorption crosssection for light resonant with the F = 2 → F = 3 transition [22,23]. We then summed each image over the central 75% of the time-of-flight Thomas-Fermi distribution along the radial direction. The domain of the sum was chosen to maximize the number of atoms counted without introducing excess noise from the edges of the Thomas-Fermi distribution where the atom numbers were smaller than image noise. The resulting set of slices f (x), calibrated in atoms/pixel, were processed to derive mean and standard deviations. The lower panel of Figure 2 shows these as dashed blue and solid red lines, respectively. It is clearly visible in the figure that the m = 0 number density fluctuations were no more than 40% of the mean, while those of the m = ±1 cloud were as much as 100% of the mean. As we will show, the two noise sources in fact have the same origin.
These results are intriguing, since a Bose-Einstein condensate is not expected to have 100% variability in atom number. Similar to laser light, a BEC can be described by a coherent state, with Poissonian number density fluctuations (shot noise). In the absence of any quench, an analysis of individual pixels showed that the m = 0 cloud possessed a variance quite close to this limit. For shot noise, the ratio of standard deviation to mean atom number density ∆n /n = 1/ √n is very small for the large atom numbers present in the experiment. By contrast, the m = ±1 clouds show "super-Poissonian" (SP) noise with ∆n /n = 1. We note that SP fluctuations have been observed in related works where spin-exchange collisions are responsible for the statistical fluctuations [24,25].

B. Role of thermal atoms in the quench
Our data in the inset to the lower panel of Figure 2 also reveal that thermal atoms play a negligible role in the quench. The m = 0 cloud can be seen to possess a bimodal distribution in space, with the central, sharp peak corresponding to the Bose-Einstein condensate and a lower density, more diffuse pedestal due to the thermal atoms whose occupation fraction was 40%, corresponding to a temperature of T ∼ 400 nK.
In spite of such a large thermal population, at short hold times the m = 0 super-Poissonian spin fluctuations clearly occur only in the condensate, and not in the normal gas, as the fluctuations (the red curve) drop sharply to zero at the Thomas-Fermi radius. At longer hold times t = 48 ms (not shown in the figure), the m = ±1 populations had grown substantially to comprise a fraction 0.6 of the total condensate. By this time, a thermal gas of m = ±1 atoms had become populated via interactions between the thermal cloud in m = 0 and various condensed spin components [26]. Nonetheless, the m = −1 spin fluctuations (and therefore, spin relaxation dynamics) still continued to occur only in the condensate, and not in the normal gas.
Our results can be summarized as coherent spin evolution for short times, followed by thermally assisted spin redistribution occurring at long hold times. The separation of timescales poses the intriguing possibility to cool the sample via spin-changing collisions and selective spin state removal of the thermal cloud.

C. Hanbury Brown-Twiss correlations
Turning now to the condensed components, the super-Poissonian fluctuations mentioned earlier lead to Hanbury Brown-Twiss correlations in our spatially extended one-dimensional system. Such correlations have been widely observed from the domain of radio-frequencies where they were initially observed and used to determine the angular diameter of stars [27] to the optical domain [28]. For massive particles, these correlations have been observed with ensembles of ultracold atoms [29][30][31], as well as in the sub-atomic realm, where one can use them to extract information about nuclear structure from collisions [32]. In our experiment with Bose-Einstein condensates they reveal the coherence length associated with the non-equilibrium state created by the quench.
We examined the second order spin density correlation function g 2 (x) with relative coordinate x between points of observation: Here n(x 0 ) is the density of a particular spin component, for example, the m = −1 atoms, and x 0 represents a spatial position within the cloud. The inner and outer brackets refer to ensemble and spatial averages, respectively. Since the Thomas-Fermi density profile breaks translation invariance in the system, we define ensemble and spatial averages to be distinct quantities. First ensemble averages were taken over a restricted set of images with reduced atom number fluctuations as discussed earlier, after which a spatially averaged g 2 was computed over all values of x 0 . For each data set we performed bimodal fits to the central, m = 0 cloud to determine the condensate Thomas-Fermi radius R T F . From the mean of the m = ±1 slices we computed the central position x c of the m = ±1 clouds in order to undo the Stern-Gerlach expansion. From x c and R T F we could generate distributions n −1 (x) within the cloud. Since one divides by the mean value to compute g 2 , its value blows up as one approaches the distribution edges-to avoid this, we restricted our analysis to −R T F /2 < x < +R T F /2 and used periodic boundary conditions to compute the averages over x 0 [50].
In optics, the equivalent of the spatial variable x is the delay time τ between light paths in an optical system. In this realm, a pure coherent state has Poisson fluctuations which result in g 2 (x) = 1 for all values of x [28]. Excess fluctuations are evidence of "photon bunching"-a tendency of photons to arrive at the same time due to their bosonic nature, resulting in g 2 (0) > 1. Thermal, or chaotic light fields such as are emitted by a lamp or other low coherence source, possess g 2 (0) = 2, with g 2 (x) → 1 as the delay time exceeds the coherence time.
Our experiment clearly demonstrates the correspondence with thermal light, although, as shown earlier, the thermal cloud itself plays no role in the instability. In Figure 3a we have evaluated the normalized g 2 (x) for our experimental data using the m = −1 and m = 0 atom distributions at a time 16 ms after the quench. Corresponding data for m = +1 was similar to −1. The data clearly show a very strong "spin bunching" effect, as individual m = −1 atomic spins tend to be co-located, with g 2 (0) = 1.9 > 1. This was observed very early in the quench, when the m = +/ − 1 atom numbers were very small, and the corresponding fluctuations proportionally large. The measured value of g 2 (0) close to 2 is consistent with a thermal state and the super-Poissonian noise in the population.
In contrast to the above bunching phenomenon, the m = 0 spins exhibited very little bunching-the measured g 2 (x) was very close to 1 for all values of x, consistent with a Bose-Einstein condensate in a coherent state. Upon closer examination, we determined that g 2 (0)−1 ≈ 0.01. This small excess in the normalized variance can be largely explained by technical noise caused by atom number fluctuations as well as spin-exchange noise in the m F = 0 population due to the production of m = ±1. The former (latter) had a standard deviation of 7% (5%), and were thus of the same magnitude at this early hold time. For only slightly later times of t ≥ 20 ms, the stochastic fluctuations of the quench became 40% as noted earlier in Figure 2, and dominated over technical noise.
For the m = −1 atoms, Figure 3b shows that g 2 (x) rapidly decays in space, exhibiting damped oscillations that approach a value of 1. Defining the correlation length, x 1/2 , through the formula This correlation length is substantially smaller than the condensate Thomas-Fermi radius R T F = 340µm, and reveals the average spatial extent of the spin modes excited by the quench. The oscillations in the g 2 function indicated the formation of multiple domains simultaneously, and is further evidence of the multi-mode character that we explore in the next section. Individual images revealed domains between 15 − 45µm in size (see Figure 6). The middle panel of Figure 3 shows that for slightly longer hold times where the ratio of m = ±1 to m = 0 populations became appreciable, g 2 (0) → 1 indicating the formation of a more stable, non-equilibrium state. In spite of this, the population in ±1 had not yet reached its equilibrium value. The lower panel of Figure 3 shows that the correlation length shrinks by a factor of about 2 with hold time (although an oscillation in the population seen in Figure 1 causes it to momentarily increase at 34 ms). Thus the system transferred energy from long to short wavelength modes as the quench progressed.

D. Boltzmann statistics for the spin
How does a thermal state appear within a Bose-Einstein condensate, whose dynamics are governed by quantum mechanics? A closely related, and more general question is how isolated quantum systems thermalize when placed out of equilibrium [2]. One answer that has emerged in recent years addresses the similarity between quantum and thermal fluctuations, particularly when one . Each graph uses 20 bins covering the region of non-zero data, approximately 0 < n < 3n(x0). Good agreement is found with the Bose-Einstein distribution (Eqn. 3, solid curves) which uses the local sample averagen(x0) but otherwise has no fitting parameters. looks at the system locally. That is, even though the entire system is quantum, and has a pure case density matrix ρ, for which Tr(ρ) = 1, a subsystem A will have a reduced density matrix ρ A =TrĀ(ρ) that is mixed. Here, the trace is taken overĀ, the part of the system not in A. Measurements made within A will be indistinguishable from those made on a global thermal ensemble, since the entanglement that is generated between A and the rest of the system by the quench is not detected. This notion has recently been tested in site-resolved optical lattices, where the sub-system consists of a finite number of sites [33].
In our case "local" refers to a sector within the spin space of all particles. In particular, since the m = 0 atom number is 10 6 − 10 7 , these atoms behave classically. The quantum behavior is restricted to the |m| = 1 sector, within which there is entanglement between +1 and −1 spins generated by the quench. A measurement of both spins together shows strong quantum correlations, and has been observed previously [4,24,25,34]. Measurement of a subspace consisting of just one of the spins should result in a mixed case density matrix. Thus our experiment realizes quantum thermalization within spin space, analogous to the real space thermalization of Kaufman et al. [33].
A theoretical prediction of Mias et al. [35] elaborates upon this idea. Using a Bogoliubov treatment they showed that in a quench experiment, if one observes either of the two spin states, +1 or −1, the result is a Bose-Einstein probability distribution for the number of atoms n k in the k-th spatial mode, wheren k is the mean number of atoms in that mode, a number that grows exponentially with time subsequent to the quench. In the above formula the latter approximation holds forn k 1, which holds for all of our experimental data. From Eqn. 3, we can see that the distribution of just one of the spins should obey Boltzmann statistics with an effective temperature T ∝n k .
In Figure 4 we make a direct experimental comparison with the predicted probability distribution, Eqn. 3, at a hold time of 20 ms after the quench. Shown are probability histograms for the number of detected m = −1 atoms at 3 different spatial locations within the cloud, x 0 = −130µm, 0 and +130µm, where the mean atom number,n, varied due to the inhomogeneous Thomas-Fermi density distribution. We used different spatial locations x 0 spaced by much more than the correlation length of Figure 3c, in order to demonstrate that the probability distribution is a local quantity, and varies throughout the cloud. The theoretical prediction from Eqn. 3 is plotted as a solid line. It uses this sample mean as its only adjustable parameter. The agreement between the data and theory in each case is quite good. To generate sufficient statistics to generate an entire probability distribution from our limited data set, we used a 15 pixel (100µm) wide sample centered at x = x 0 , and 10 experimental runs whose atom number fluctuations had been filtered to < 5%, as discussed earlier. Thus each graph had 150 data points, and the mean atom numbern(x) was evaluated at x = x 0 . By averaging over spatial pixels we necessarily included data with different values of n, by ±50% for x 0 = ±130µm and ±10% for x 0 = 0. In spite of this, the local exponential character of the distribution clearly persists, and reveals the thermal statistics of the spin states produced by the quench.

IV. THEORY
Having directly generated the modes experimentally, we turn now to their theoretical description. We focus our theory on both q < 0 and q > 0 with zero magnetization. Our effort closely parallels that of other experimental observations of spinor instabilities, with important differences. For example, Bogoliubov theory was applied to a finite q > 0 instability of ferromagnetic F = 1, 87 Rb spinor BEC [9], as well as to the q = 0 instability [20] and other instabilities [13] of antiferromagnetic F = 2 spinor BEC. Broadly speaking, these works have identified instabilities arising either through bulk modes with a finite wavevector, as in [9,13], or a specific mode or set of modes that are resonantly excited at certain values of the quadratic Zeeman tuning parameter [20]. Our studies, by contrast, explore an intermediate regime. A bulk analysis assuming spatial homogeneity fails to capture essential features of our observations, particularly the localized instability near the trap center. However, neither is our experiment dominated by the discrete mode structure of the trap, as the relevant modes along the long axis of the cigar are too closely spaced for us to resolve. Instead, in our specific experimental geometry, the inhomogeneous density profile plays an important role in shaping the unstable modes. We uncover these modes by solving the Bogoliubov equations directly in coordinate space.
Bogoliubov theory was first applied to multicomponent (spinor) BEC separately by Ho [36] and Ohmi and Machida [37]. Following their approach and others [18], one linearizes the spinor Gross-Pitaevskii (GP) equations (or the corresponding Heisenberg equations of motion for the field operators) about an initial state that is classical. In our case and several of the examples above, this is a state ψ 0 consisting of all atoms in the m = 0 sublevel, with only small corrections δψ ±1 describing the populations in m = +1 and −1. Due to the small ratio c 2 /c 0 , we assume that the spin instabilities do not couple strongly to density fluctuations, and thus we can neglect fluctuations in the m = 0 state. The resulting spinor wavefunction may be written as and the resulting linearized spinor GP equation for m = ±1 is [38]: In the above we have simplified the notation to ψ m ≡ δψ m , and assumed a one-dimensional description with axial coordinate x, so that U (x) = c 2 n 0 (x) = c 2 |ψ 0 (x)| 2 .
To model the data presented in this paper we set the linear Zeeman term p(x) = 0 and allow the quadratic Zeeman shift q to vary as a free parameter. We exclusively study the regime very close to the phase transition, i.e., −U 0 q < 0, where U 0 = c 2 n 0 (x = 0) and n 0 (x = 0) is the peak m = 0 atom number density at the trap center.
We expand the wavefunctions in a basis of spin excitations, with spatial mode index k and frequency ω k : Note that since the above is actually two equations, one each for m = ±1, there will be two spin modes associated with each spatial mode k. Putting Eqn. (5) into Eqn. (4) and equating terms with equal time-dependence, we get the Bogoliubov-de Gennes equations:  [39]. For the spin wave case, however, the dynamics are much slower than for sound waves by the factor c 2 /c 0 (≈ 8 for sodium atoms), as first elucidated by Ho [36] and Ohmi and Machida [37]. Thus the lowest excitation frequencies are typically in the range of a few Hz to 10s of Hz. We find the eigenvalues and eigenvectors of Eqns. (6) using a straightforward matrix diagonalization in MAT-LAB [40]. Figure 5 shows the numerical solutions for the energy eigenvalues for some representative parameters. For q > 0 the eigenvalues are real. If q is chosen to be negative, one or more collective modes in Eqn. (5) has an energy eigenvalue E k which crosses into the complex plane. This triggers an exponential growth in the population of those modes, which are linear combinations of the spin states ±1. Thus the populations ψ † m ψ m , for m = ±1, also grow exponentially with time, similar to a parametric amplifier [41]. Although we have not written down the Bogoliubov expansion in terms of the field operators ψ m , it is straightforward to do so, and all quantum effects and correlations can be calculated in a straightforward manner [18].
Before turning to solutions to the equations, we point out some differences between the coordinate and momentum representations. For uniform systems, Bogoliubov theory is best described in momentum space, using plane wave modes. One can then write the annihilation operator for a boson with momentum k in terms of corresponding operators for quasiparticles with momenta k and −k. The Bogoliubov transformation contains within it, therefore, a direct correlation between quasi-particles of opposite momenta. This correlation is similar to that obtained in the Bogoliubov diagonalization of the Hamiltonian of a single component weakly interacting Bose gas [42]. In addition, in the multicomponent Hamiltonian, Eqn. (1), the correlation is between quasiparticles of opposite spin.
In our experiment, where the m = 0 condensate has a Thomas-Fermi spatial density profile, an expansion in momentum eigenstates is not useful. Instead, we have followed the approach of Ruprecht et al. in the analysis of collective excitations of a scalar BEC in a trap [43]. In that case, the Thomas-Fermi density profile led to collective mode functions that were spatially varying, and which represented modes located inside of or near the Thomas-Fermi surface.
We will find the same to be true of the collective spin modes for a spinor BEC under harmonic confinement. An alternate way to view these modes is in terms of standing wave solutions u k,m (x), v k,m (x) for the small excitations ±1 that are created within the Thomas-Fermi boundaries of the m = 0 cloud (see Figure 6A for an example). Thus, rather than momentum correlations, as expected for a uniform system, the Bogoliubov analysis reveals the spatial correlations for particles of opposite spin m = ±1, as noted in Figure 2 and our earlier work [4]. The correlations only exist, however, within the domains defined by those modes.

A. Uniform density profile
The case of a uniform m = 0 density, U = constant, is a useful point of reference since the solution can be analytically obtained. The energy spectrum in this case is where k = 2 k 2 /(2M ), k = π/L × n, n = 1, 2, 3, ..., for excitations in a box of length L = 2R T F , where R T F is the axial Thomas-Fermi radius. The Bogoliubov eigenfunctions are box modes for |x| < L/2. For our parameters the ground state energy of the box 1 h× 0.05 Hz is smaller than our experimental resolution, so we can assume a quasi-continuous spectrum. Thus for q > 0 all eigenvalues are real, while for q < 0 imaginary eigenvalues define unstable modes. For short times after the quench, the amplitude of these unstable modes grows exponentially in time with a rate Γ = |Im(E k )|/ . The maximally unstable mode is defined to be the one whose imaginary component is the largest, i.e., Γ = Γ max . For −U < q < 0, maximizing Γ yields the mode It has a wavevector k max = π/L, i.e. a wavelength twice the Thomas-Fermi length of the condensate. Neglecting the zero point energy, the corresponding instability rate is

Thomas-Fermi density profile
To solve the Bogoliubov Eqns. (6) in the inhomogeneous case, we expand u, v in the box basis (Eqn . (8)), which for p = 0, yields the matrix equation where the basis size was held to N elements, u ν , v ν are the box basis coefficients for ν = 1, 2, 3, ...N , and a summation over ν is implied in the matrix product. In this basis, the matrix elements of the operators in Eqn. (6) are H 0 µν = ( µ + q)δ µν and U µν = U 0 The box eigenenergies are µ = µ 2 1 . The numerical problem consisted of diagonalizing a square matrix of order 2N , with N pairs of eigenvalues E, −E corresponding to the pair of spin modes discussed earlier. For each value of q and given the values of 1 , U 0 , the eigenvalues and eigenvectors were found numerically using MATLAB. The routines were tested against the exact solutions, Eqns. (7) and (8), by fixing the density to be a constant. Typically, the ground state energy was found to converge to 10 −5 using 150 basis elements. Unless otherwise stated, the parameters used were U 0 /h = 96 Hz and 1 /h = 0.0047 Hz. These values corresponded closely to typical experimental parameters.

C. Eigenvalue Spectrum
We first discuss the numerical solutions for the eigenvalues before proceeding to the eigenvectors. The left panel of Figure 5 shows the eigenvalue spectrum computed for p = 0, for q = +5, 0 and −5 Hz, corresponding to stable, critically stable and unstable regimes. The eigenvalue spectrum in the case of a homogeneous density gas at the same peak value of n 0 is plotted on the right panel for comparison. In both cases we have plotted the square of the energy eigenvalue, E 2 , rather the eigenvalue E itself, since according to the Bogoliubov equations this quantity is always real for p = 0. The data in Figure 5 highlight important differences between the two cases, as well as some of their similarities.
For q > 0 (stable regime, E 2 > 0), the graph shows that the lowest eigenvalue (the mode with index n = 1) is finite for both inhomogeneous and homogeneous densities, and thus the system always has an energy gap. For a homogeneous system this gap is given by where is the lowest box mode. For the parameters used, E 1 = h × 31Hz. For the inhomogeneous case the numerical data show that the energy gap is reduced by a factor of nearly 3, to h×12Hz. As we will see later in Section V A, the lowest energy Bogoliubov eigenfunctions have a very different spatial profile for the two cases-for the inhomogeneous case, these modes are sharply localized near the Thomas-Fermi boundary (see Figure 10), whereas for the homogeneous case they are mostly near the center of the cloud (see Eqn. (9)), which raises the homogeneous energies due to repulsive interactions with the m = 0 atoms. With modes localized near the boundary rather than at the center, the energy difference between modes of opposite parity becomes negligible since they both have nearly the same overlap with the density profile U (x). Thus, the low-lying modes all come in nearly degenerate pairs that we term a "parity doublet" (the odd parity modes remain slightly higher in energy than the even parity ones, but the difference is very small for n < 20). These doublets can be seen in the data as a series of pairs of dots. For the homogeneous case, by contrast, the near degeneracy between even and odd states disappears, so there is no parity doublet.
For low n the dispersion relation-the variation of E with n-is also considerably different in the two cases. For example, an inflection point can be seen at n 25 in the dispersion relation for the spin modes for the inhomogeneous distribution. For n < 25 the dispersion relation has negative curvature, similar to surface modes of a single component BEC [39]. This stands to reason, as the modes are localized near the Thomas-Fermi surface. For n > 25 the curvature is positive, suggesting a transition to nearly free-particle behavior, although the eigenvalues are still smaller than the height of the potential U 0 . What determines the crossover point is still unclear, although the energy eigenvalue at the inflection point was observed empirically to increase with q.
As q becomes smaller, the entire spectrum of E 2 shifts to smaller values, until the lowest eigenvalue reaches zero at a point very close (within ) of q = 0, which defines the boundary of the unstable region. However, as we have plotted E 2 rather than E, the transition from stable to unstable regimes becomes a more smooth and continuous one, with the number of imaginary eigenvalues (E 2 < 0) increasing as q decreases. For example, at q = −5 Hz 32 eigenvalues have become imaginary, while for n ≥ 33 they still remain real. Due to the tiny value of , it plays no role, and the phase transition occurs at the same point for both inhomogeneous and inhomogeneous cases.

V. COMPARISON OF THEORY WITH EXPERIMENT
We can separate the temporal dynamics into two phases. In the early, growth phase, the population fraction in ±1, f ±1 , increases with time, but always remains small compared to that of the m = 0 state, f 0 . Bogoliubov theory can be used to study this phase. The second, dynamical phase, occurs when all 3 components, f +1 , f −1 , f 0 , have the same order of magnitude, and interact strongly with one another. The dynamical phase is not captured by the Bogoliubov theory, but can be observed experimentally and compared with numerical simulations. We discuss these two phases of the dynamics separately below.

Growth phase
In the growth phase, we computed the complex eigenvalues E n . These were sorted by the magnitude of their imaginary component, with n = 1 mode having the largest imaginary component, n = 2 the second largest, and so on. Figure 6A shows the numerically obtained solutions for the maximally unstable eigenvector, n = 1, as well as that of a less unstable eigenvector, n = 6, at a quadratic shift of q = −4.2 Hz, corresponding to the value used in the experiment. The boundaries of the plot are the Thomas-Fermi surface, x = ±R T F = ±340µm. The maximally unstable mode is localized near the center of the cloud, with a wavefunction resembling a Gaussian profile. This is to be contrasted with Eqn. 9. For larger n the number of nodes increased, as did the spatial domain over which the mode function was non-zero.
We define an rms width for each unstable mode as where u n is the n-th mode, and the upper and lower limits of integration are the Thomas-Fermi surface, ±R T F , respectively. For the n = 1 mode, x rms = 0.05R T F , nearly four times smaller than the homogeneous case, Eqn. 9, whose rms width is 0.18R T F . Thus the inhomogeneous density profile had a profound effect on the instability, causing the nucleation of localized spin domains. This tendency can also be envisioned by applying a local density approximation to the rate Eqn. (7), which reflects the inhomogeneous gain profile for the spinexchange process. Since the m = 0 state is dynamically unstable for q < 0, excitations develop at a rate that depends upon the local value of its density, c 2 n 0 (x). Thus the spin domains become localized near the cloud center, where this rate is highest. A similar effect was observed in reference [6], as noted earlier.
Another, equivalent way to view this is through consideration of how the instability amplifies spin noise, which can be represented in any basis. Using plane waves with momentum p, for example, and if the system contained noise that was uniformly distributed at all spatial frequencies up to p max / = 2π/ξ sp , where ξ sp = 1.5µm is the spin healing length, it would begin in a delocalized state where the average amplitude was roughly the same everywhere. This picture is consistent with the Truncated Wigner Approximation (TWA), as we discuss later. As time develops, only the particular superposition of momentum eigenstates that reproduces the state u 1 (x) ≡ u M AX (x) shown in Figure 6A would be amplified significantly, and the spin distribution would develop into something that is spatially localized. For longer times where the m = 0 component becomes depleted, this spatially localized state is no longer stable, but expands into a multi-domain structure.
Indeed, we can see from our experimental data in Figure 2 that the instability creates spin structures that are spatially localized near the center of the Thomas-Fermi region. Through the second order correlation function we determined the average mode size to be 33µm, which is in good agreement with the rms width of 35µm predicted for the maximally unstable mode u M AX (x). The latter is shown as a horizontal line in Figure 6G, right panel.
From Bogoliubov theory at q = −4.2 Hz, we estimate that 30 modes have an imaginary component. Therefore, the maximally unstable mode is not the only active mode in the problem. While the average domain size, as measured by the normalized second-order correlation 2 Hz shows an average domain size that is smaller than that of the lowest mode, suggesting the involvement of higher lying modes n > 1. function g 2 (x), is well captured by theory, we still must understand the varying number and size of the domains measured in the experiment. To uncover these multimode effects, we look at the spread in experimentally observed domain sizes using an rms domain width analysis. Figure 6, panels (B-G), show mode profiles measured at the onset of the instability, t = 20 ms, when the mean population in the ±1 states was 15%. We show 5 separate instances of the experimental quench sequence that are representative of the variations observed. Shot to shot fluctuations reflected the stochastic dynamics discussed earler. By observing peaks in the data, we could determine the rms size of domains associated with those peaks. We determined these sizes using a simple criterion-the size was the minimum distance from the peak where the data crossed a threshold value of = e −1/2 of its peak value, as would be expected for the rms width of a Gaussian function. The threshold is shown as a solid line in the plots. The preponderance of multi-domain structures makes it clear that multiple modes are present. Figure 6G explores this multi-mode character of the instability by comparing the domain widths quantitatively. On the right panel we have computed the rms width of the maximally unstable mode obtained from the Bogoliubov calculation, u M AX (x), versus the final quadratic Zeeman shift |q| = −q. On the left panel we show a histogram of the observed rms sizes of the domains for 30 runs of the experiment at q = −4.2 Hz. On the higher end, the distribution cuts off at a domain size that corresponds well with x rms = 35µm for the lowest mode, u = u M AX . Nonetheless, domains as small as 15µm were observed, not much larger than our spatial resolution of 10µm, indicating that the quench excited many higher order modes. The tail in the distribution above 35µm is most likely caused by mistaken identification of multiple overlapping domains as a single, larger domain, an example of which is shown in panel (C).

Dynamical phase
In this later phase, we observe a coarsening of the domains generated in the growth phase. The upper panel of Figure 7 shows the crossover from growth to dynamical phases in the one-dimensional m = −1 density profiles, which have been plotted against their axial position normalized to the axial Thomas-Fermi radius, R T F . Similar data were obtained for the m = +1 profile. Each curve is an average of 3 experimental runs normalized to the peak value at each time step, with each curve displaced by 1 for clarity. The average reduces the effect of stochastic fluctuations associated with spontaneous domain formation, allowing us to focus on the coarsening trend-there is a growth in the overall size of the m = ±1 clouds with time. As seen in the figure, for short times, when Bogoliubov theory is still applicable, the density profile grows from the center of the cloud, forming a localized hump at a time when f pm ≈ 0.05. As time increases, this hump grows in size to envelop the cloud.
Although we can no longer use the unstable Bogoliubov eigenmodes to analyze the dynamical phase, we can compare our data with numerical simulations. To this end we compute a different rms width, this one pertaining to the entire m = −1 cloud, which is shown in the lower panel of Figure 7. Here |u(x)| 2 is equal to the measured density profiles shown in the upper panel of the figure. The width of the density hump is seen to increase with time as the system evolves into the dynamical phase, exhibiting an overdamped oscillation before reaching a steady-state value of about 0.4R T F . At t = 16 ms the super-Poissonian noise was larger than at later times, and imperfect averaging led to a larger error bar. Although the coarse cloud size in each of the three interpenetrating quantum fluids, m = 0 and m = ±1, has reached a steady-state, the dynamics have not yet halted, as the fraction f ± continues to steadily increase, as seen earlier in Figure 1C. In this later phase of the non-equilibrium behavior micro-domains still exist and move throughout the cloud. As noted in our earlier work, one can observe a small, local magnetization density M (x) = n +1 (x) − n −1 (x). Thus the m = ±1 clouds eventually separate from one another [4].

Numerical results
Numerical simulations allowed us to bridge the growth and the dynamical phases of evolution, and to probe the early part of the growth phase where little experimental data could be obtained. We performed one-dimensional simulations of the 3 coupled spinor Gross-Pitaevskii (GP) equations seeded with noise according to the Truncated Wigner Approximation, or TWA. The population dynamics derived from the simulations were reported in our earlier work [4]. Here we provide more details including the numerically obtained wavefunctions and their temporal evolution. Our numerical procedure is a straightforward forward time propagation of the equations of motion using a time-splitting spectral method (see, for example, reference [44] and references therein). The TWA approximation is expected to be valid for both short and long times, as long as the initial condition is a classical state [45][46][47][48]. To implement these simulations, as discussed in [4], we assumed a BEC initially at zero temperature and obtained the initial wavefunction for the m = 0 component numerically. Vacuum noise in the m = ±1 states was simulated as classical noise, and we computed the average density, ψ † m ψ m , as an ensemble average over 30 separate simulations using different random initial conditions. Vacuum modes with wavelength less than ξ spin are not expected to contribute to the spin instability. Therefore, we imposed a cutoff energy of c 2 n 0 , which resulted in N v 700 virtual particles, while the condensate contained 5 × 10 6 particles, similar to the experimental conditions. To study the early time behavior in the simulations it was also essential to subtract a constant from the average density = N v /(2R T F ) equivalent to the sum of all virtual particles added, which was done according to the Weyl representations of the field operators [46].
As mentioned earlier, the experimental data measure the integrated column densityñ(x) = n(x, y, z)dydz.
To compare this with a one-dimensional simulated density profile n(x), we first posit a solution that is separable in space between axial (x) and transverse (y, z) coordinates: where ψ m (x, t) is the simulated wavefunction for spin state m. Since the quench is one-dimensional, we may assume that the transverse mode function ξ, normalized as |ξ| 2 dydz = 1, is time-independent as the dynamics are frozen. With the approximation Eqn. (12) and the threedimensional density distribution n m (x, y, z) = |Ψ m | 2 , the measured column density becomes n m (x) = n m (x, y, z)dydz = |ψ m (x, t)| 2 and is identical to the one-dimensional density profile obtained from the simulation. However, for a Thomas-Fermi BEC, the solution is not strictly separable, as the transverse Thomas-Fermi radius depends on axial position x. In effect, Eqn. (12) assumes a Thomas-Fermi cylinder rather than a cigar. If the cigar aspect ratio is very large (70 in our case), the difference between cylinder and cigar is insignificant, particularly since most of the important quench dynamics occur near the cloud center. However, care must be paid when one approaches the axial Thomas-Fermi radius, x = ±R T F , as there are likely to be quantitative differences.
On the left panel of Figure 8, we show the result of the numerical simulations, and on the right pane is shown the prediction for the ensemble averaged density profile from Bogoliubov theory. Here we have computed |u k | 2 e Γ k t where the sum runs over all unstable modes. The growth factor Γ k = Im[E k /h]. The curves have been normalized in the same manner as for the experiment. At t = 0 the exponential factor is 1 for all modes, resulting in a nearly uniform initial density profile.
Both Bogoliubov and full numerical theories agree with one another during the growth phase of the dynamics. Moreover, the theory confirms the local density picture discussed earlier, where a uniform distribution comprised of vacuum fluctuations in all modes becomes narrower with time, eventually becoming localized near the cloud center. For the Bogoliubov results, at later times the maximally unstable mode, Γ M AX , begins to dominate, and the curves begin to peak around this mode function, u M AX , which is shown as the narrow distribution in the uppermost plot of the right pane. Even at a time  5) is shown on the right, for times t/τMAX = 1, 2, ...10, with the maximally unstable eigenvector, |uMAX | 2 , shown in the uppermost graph for comparison. Here τMAX = 1/Im(EMAX ) is the time scale associated with uMAX . In both panes, each curve has been normalized to its peak value and displaced for clarity of presentation. t = 10/Γ M AX , however, the Bogoliubov density distribution has not fully converged to u M AX , but remains broader. Our experimental data shows the formation of the localized structure, but not its precursor, the uniform phase, which is hidden in experimental noise. The uniform phase is, however, captured by the TWA simulations (see the lowest traces of the left panel).
An interesting artifact in the simulations can also be observed in Figure 8, one which illustrates some of the limitations of the Bogoliubov analysis. The TWA initial condition was taken to be a sum of Bogoliubov eigenmodes prior to the quench (see the next section for details of these stable modes), with random coefficients. Since these modes are defined on x ∈ [−R T F , +R T F ], their amplitude goes exactly to zero at the Thomas-Fermi radius.
However, the numerical simulations are not restricted to the Thomas-Fermi volume, but capture the full details of the cloud's surface structure, even for |x| > R T F . Thus at short times in the simulation, the repulsive interaction between atoms redistributed the m = ±1 density from inside to outside of R T F such that as x approached ±R T F from within the cloud, the Weyl correction was no longer accurate and yielded a negative density, seen in the lowest traces. Therefore, at short times the m = ±1 density should be even more uniform than the simulation suggests. This artifact had no bearing on the simulation at longer times, since the Weyl correction, which counts only the vacuum fluctuations, was insignificant in comparison with the number of real particles. Nonetheless, it illustrates the difficulty in describing the details of the dynamics near the Thomas-Fermi surface. Finally, an agreement between all 3 density profiles-the experimental data, TWA simulations, and Bogoliubov theory-was found at a time approaching the crossover from growth phase to dynamical phase. This data is shown in Figure 9. The experiment and TWA curves were taken at a time when the ±1 fractions were similar, having reached f ±1 = 0.15 and 0.17, respectively. This allowed us to circumvent a factor of 4 difference observed in the absolute timescale for the quench dynamics between the two, as noted in [4]. These data contain residual oscillations due to imperfect averaging. The Bogoliubov theory was taken at a time t = 10τ M AX . All 3 curves are broader than the maximally unstable eigenmode u M AX . Bogoliubov theory should converge to the maximally unstable eigenmode; however, this only occurs after a sufficiently long time. We observed convergence at a time t ≈ 100Γ M AX , by which time in the experiment the m = 0 cloud would have been significantly depleted, violating the Bogoliubov approximation. This provides further confirmation that our experiment is in a multi-mode regime even during its growth phase, when Bogoliubov theory is valid.

A. Stable Modes
For completeness, we include a discussion of the collective excitation spectrum for q > 0, although no experimental data was taken in this regime. Nonetheless, it provides additional insights into the difference between homogeneous and inhomogeneous cases, particularly as q → 0, the phase transition point. Here the eigenvalues are all real and positive, and n is a mode index by which they are sorted in increasing order. Similar to the unstable modes discussed in the previous section these modes also depart from the homogeneous Bogoliubov solutions, Eqn. (8). Their spatial profile also depends upon q in a manner that was determined numerically.
The lower graph of Figure 10 shows the numerically obtained mode functions, u 1 (x) and v 1 (x), for the lowest energy eigenvalue E 1 , at a quadratic Zeeman shift of q = +5 Hz. These functions are sharply localized near the Thomas-Fermi boundary at x/L T F = ±1/2, in stark contrast with the homogeneous modes that are delocalized throughout the Thomas-Fermi region. For increasing n the modes penetrate further into the cloud-for comparison, the n = 9 mode, with even parity, is shown in the upper graph.
We can understand the mode structure for q > 0 in terms of the total potential appearing in the Bogoliubov equations: where n 0 is the density averaged over the 2 radial directions of the optical trap. Beyond the Thomas-Fermi radius we can set U (x) = ∞, since the harmonic potential increases very rapidly in comparison with the energy scale c 2 n 0 . Figure 10 shows U (x) as a dashed line. It resembles the box potential obtained in the homogeneous case, but with an additional bump at x = 0 that shifts the excitations away from the cloud center and toward the Thomas-Fermi boundary. This repulsion was discussed earlier as being due to antiferromagnetism-for c 2 > 0 the spin m = 0 and spin |m| = 1 quantum fluids repel one another. Since the excitations are ±1 atom pairs, their minimum energy configuration is a localized state at the edge of the cloud.
Crossing the phase transition results in a transformation from stable to unstable behavior of the eigenmode. This effect is explored in Figure 11 for the lowest energy state. Only within a tiny region near the critical point having a width of 0.1 Hz, does the mode function be- come delocalized. By contrast, in the homogeneous case the mode sin π(x + L/2)/L remains the same on both sides of the phase transition, and is always maximum at x = 0.

VI. CONCLUSION
For spatially extended quantum systems, the study of relaxation toward equilibrium naturally involves the dynamics of many modes and the flow of energy between them. We have used the second-order correlation function, g 2 (x), and statistical analysis of domain widths to reveal the richness of this multi-mode behavior in quenched antiferromagnetic spinor Bose-Einstein condensates. These approaches, combined with a portfolio of theoretical tools, have allowed us to span the data from the early, growth phase, to the later, dynamical phase. For the former case, Bogoliubov theory is a semi-analytical approach which provides much physical intuition. However, for the latter case we have relied on numerical simulations in order to explain the data. Future work will explore possible dynamical universality in the long time dynamics, which could provide new analytical insights beyond the numerical work that has been done here [49].
This work was supported by NSF grant No. 1707654.