## Abstract

A proposed paradigm for out-of-equilibrium quantum systems is that an analog of quantum phase transitions exists between parameter regimes of qualitatively distinct time-dependent behavior. Here, we present evidence of such a transition between dynamical phases in a cold-atom quantum simulator of the collective Heisenberg model. Our simulator encodes spin in the hyperfine states of ultracold fermionic potassium. Atoms are pinned in a network of single-particle modes, whose spatial extent emulates the long-range interactions of traditional quantum magnets. We find that below a critical interaction strength, magnetization of an initially polarized fermionic gas decays quickly, while above the transition point, the magnetization becomes long-lived because of an energy gap that protects against dephasing by the inhomogeneous axial field. Our quantum simulation reveals a nonequilibrium transition predicted to exist but not yet directly observed in quenched s-wave superconductors.

## INTRODUCTION

The challenge faced in understanding out-of-equilibrium systems is that the powerful formalism of statistical physics, which has allowed a classification of quantum phases of matter based on simple principles such as minimization of free energy, does not apply. A diverse range of nonequilibrium phenomena has been observed, including synchronization (*1*–*5*), self-organization (*6*–*9*), quantum chaos (*10*, *11*), Loschmidt echo singularities (*12*), and time crystals (*13*, *14*). A proposed organizing principle is that transitions, reminiscent of those found between thermodynamic ground states, can also be found between dynamical phases (*15*–*19*).

In general terms, a nonequilibrium phase transition is characterized by the existence of a critical point separating phases with distinct dynamical properties in many-body systems. The analog of thermodynamic order parameters is found in long-time average observables, which have a nonanalytic dependence on system parameters. In driven open systems, energy and particle number are not necessarily conserved, and nonequilibrium transitions are typically signaled by different steady states that occur upon varying system parameters such as pump or loss rates (*20*–*22*), independently of initial conditions. In closed systems, dynamics are often initiated by quenching control parameters, with qualitatively distinct behaviors observed below, above, or at a critical point (*15*–*18*) that can depend on the initial state of the system. The label “dynamical phase transition” has been applied not only to the boundary between two dynamical phases but also to the nonanalytic behavior in real-time dynamics of the return probability amplitude (*12*, *23*), which does not require an order parameter to be defined (*24*). The phenomenon under investigation in our work is the former case, which we will refer to as a transition between dynamical phases (TDP) to avoid confusion. The theoretical study of these transitions has encompassed a broad range of platforms including collective spin models (*25*–*27*), nonequilibrium phases of superconductors (*28*–*30*), interacting fermions and bosons on the lattice (*15*–*17*, *31*, *32*), and quantum field theories (*18*, *33*); however, experimental investigations have so far been restricted to self-trapping transitions in bosonic systems (*34*–*38*) and the transverse-field Ising model realized with trapped-ion chains (*19*).

Here, we report the observation of a transition between two dynamical phases of a quantum degenerate Fermi gas. The sample under investigation consists of neutral potassium atoms (^{40}K) confined in a harmonic optical trap and cooled to nanokelvin temperatures. The controllable interactions of this closed quantum system enable a broad search for nonequilibrium phenomena that arise from the interplay of atomic contact interactions, quantum statistics, and motion. Using collective magnetization as an order parameter, system dynamics are observed directly and compared to theoretical models at various levels of approximation.

We understand and analyze our system through a mapping of the single-particle eigenstates of the harmonic trap onto a lattice in mode space, as depicted in Fig. 1. By tuning the interaction strength to suppress collisions that would change the occupancy of the modes, the atoms become pinned on the conceptual lattice, enabling the description of our system with a spin model (*39*–*41*); here*i*th atom and *X*, *Y*, and *Z* denote orientations in Bloch space. This is the collective Heisenberg model, a canonical model for magnetism (*42*) in which the nonlocal spin-spin couplings *J _{ij}* compete with an inhomogeneous axial field,

*h*(

_{i}*39*–

*41*). Similar treatments of fermionic systems using a spin model have been used successfully in optical lattice clocks (

*39*,

*43*–

*46*) in the microkelvin regime. However, in those experiments, undesirable inelastic collisions limited the number,

*N*≲ 50, and prevented the study of transitions between well-defined dynamical phases. The stability of low-lying hyperfine states and control over interactions in our experiment allow us to explore many-body dynamics in macroscopic ultracold samples of

*N*≃ 3 × 10

^{4}alkali atoms at nanokelvin temperatures and test the spin model in this new regime.

We find that below a critical interaction strength, the magnetization of an initially polarized gas quickly decays, while above this critical point, the magnetization becomes long-lived and is protected by an energy gap against inhomogeneous field-induced dephasing. These observed dynamical phases are a manifestation of an emergent property in a many-body dynamical quantum system. The transition would be absent for small particle number since the emergence of a sufficiently strong gap from weak two-body collisions relies upon a collective nonlinearity. We then implement a many-body echo sequence (*4*, *47*, *48*), which is a direct test of reversibility. This allows us to identify the boundaries of the parameter regime in which the complex far-from-equilibrium dynamics of interacting fermions is quantitatively described by the collective Heisenberg model. The successful mapping allows us to implement quantum simulations of the nonequilibrium phases predicted to exist in quenched s-wave superconductors by Richardson-Gaudin models (*28*–*30*) but not yet directly observed, given the need for ultrafast probes (*49*).

## RESULTS

The simulation cycle begins with a noninteracting sample fully polarized in the lower spin state ∣↓〉, which ensures that no site in the mode lattice is doubly occupied. Nonequilibrium dynamics are initiated by a fast radio-frequency (rf) pulse that rotates the collective magnetization into the *XY* plane. The time evolution of transverse magnetization is probed using a Ramsey sequence: Following the initial π/2 pulse, atoms evolve for a variable time *t*, after which a second π/2 pulse is applied, and the total populations in the ∣↓〉 and ∣↑〉 states are measured with a Stern-Gerlach technique. Shot-to-shot field drifts on the microtesla scale prevent a reproducible accumulated phase in the Ramsey sequence. We estimate the magnitude of the transverse coherence by repeating the sequence at least 10 times and using a maximum-likelihood estimator that assumes a randomized interferometric phase (see Materials and Methods for further details). This procedure measures the total transverse magnetization *N*, with

The optical confinement creates a potential that is approximately harmonic, with frequencies ω = {ω* _{x}*, ω

*, ω*

_{y}*} = 2π × {395,1140,950} Hz along three spatial directions. Spin-dependent curvature in the confinement potential produces a further shift ±Δω in the oscillator frequency between the ∣↓〉 and ∣↑〉 states. Since the resultant energy shift depends linearly on the index of the single-particle motional eigenstates, labeled by*

_{z}Interactions are proportional to the s-wave scattering length *a* of the colliding atoms, which are tuned by a magnetic Feshbach resonance near 20 mT (*50*). We factor *J _{ij}* as

*U*

*i*th and

*j*th particles (see the Supplementary Materials). Because of the extended nature of the motional wave functions, the

*p*=

*x*,

*y*,

*z*. The TDP is observed in a weak-scattering regime, where atoms remain frozen in their initial modes and dynamics involve only spin degrees of freedom (

*3*,

*4*,

*51*). The precise control of

*a*required an improved determination of the zero-crossing field,

*B*

_{zc}, at which

*a*= 0. As described in Materials and Methods, we find

*B*

_{zc}= 20.907(1) mT.

### Dynamical phase diagram

Figure 2 shows an exploration of the nonequilibrium phase diagram using total magnetization at 100 ms, S(*t* = 100 ms), as the order parameter. Simulations were run with scans of mean interaction strength *J* = 〈∑_{i,j} *J _{ij}*/

*N*

^{2}〉

*(Fig. 2, A to C) or scans of axial field inhomogeneity*

_{T}*i*and

*j*run over

*N*populated modes. The thermal average, 〈⋯〉

*, is performed by averaging over different realizations of populated modes drawn from a Fermi-Dirac distribution. Several distinct regions appear: a dynamical ferromagnet with high persistent S at large positive or negative*

_{T}*NJ*and a paramagnetic phase with low S at smaller ∣

*NJ*∣.

Measurements are compared to a mean-field treatment of Eq. 1, where the *j*th atom experiences an effective magnetic field, **B*** _{j}*, which depends only on the local field,

*h*, and the average magnetization of the other atoms in the ensemble

_{j}*i*and

*j*run over a set of

*N*populated modes drawn from a finite temperature Fermi-Dirac distribution, and

**Z**is a unit vector. The TDP is a consequence of the opening of an interaction energy gap between the fully polarized manifold and the remainder of the Hilbert space, as further discussed below. Numerical solutions of the corresponding 3

*N*nonlinear Bloch equations show that above some critical interaction strength,

*h*. In contrast to thermodynamic ferromagnetism, which occurs only for

_{i}*J*> 0, the dynamical ferromagnet is a spin-locked state that can be stabilized by either repulsive or attractive exchange interactions.

The red lines in Fig. 2 show S(*t* = 100 ms) calculated by this model using ab initio determination of *J*, a set of field curvatures that match all observations in this manuscript (see the Supplementary Materials), and a decay envelope [*e*^{−Γ(a)t}] discussed in detail further below. The white solid line in Fig. 2D shows the TDP steady-state (*t* → ∞ ) phase boundary that sets *t* = 100 ms ≈ 40 (ω* _{x}*/2π)

^{−1}is sufficiently long to capture the nonequilibrium phase diagram as a function of

*J*and

### Transition between dynamical phases

To gain understanding of the scaling expected near the TDP, one can use a simplified “all-to-all” model, in which coupling constants are replaced by their mean value, *J _{ij}* →

*J*. In this limit, Eq. 2 becomes integrable and maps to a Richardson-Gaudin model, a model used to describe fermionic superconductors when the Bardeen Cooper Schrieffer Hamiltonian is expressed in terms of Anderson pseudospins (

*52*). Borrowing the methodology developed for dealing with dynamical phases in superconductors (

*30*,

*53*,

*54*), one can obtain the frequency spectrum ruling the nonequilibrium dynamics from the roots of

*L*

^{2}(

*u*), where

**L**(

*u*) is the Lax vector of the auxiliary variable

*u*(see the Supplementary Materials). The roots can be found using the property that

*L*

^{2}(

*u*) is an integral of motion of the dynamics and can be evaluated for convenience at time

*t*= 0. When the roots compress in the neighborhood of the real axis, the long-time limit of S(

*t*) relaxes to a zero. This corresponds to the normal phase or “phase I” in the language of superconductors (

*30*). On the other hand, the appearance of a pair of complex conjugate roots determines the TDP critical point and the emergence of a phase characterized by a nonzero steady-state order parameter, S(∞) > 0, i.e., “phase II.” While the precise scaling of the order parameter near the TDP can be complex, since it is determined by the spectrum of the

*h*, we find that above the critical point in our system, it can be approximated by the analytic expression

_{i}*J*

_{eff}=

*J*) for the case of a one-dimensional (1D) system Δω

_{y, z}= 0 at zero temperature. To account for noncollective interactions, higher dimensions, and finite temperature, we introduce renormalization parameters α and

*J*

_{eff}= β

*J*. The critical interaction strength is

*t*) exhibits transient oscillations at the gap frequency Ω, which slowly damp as it reaches S(∞). The gap frequency goes to zero at

The collective nature of these phenomena is emphasized by an alternative interpretation of the gap. The initial π/2 pulse can be said to create a superposition of ∣*S* = *N*/2, *m _{s}*〉 Dicke states, where

*S*and

*m*are eigenvalues of the collective operators

_{S}*m*states have the same energy in the rotating frame of the rf pulse. A finite energy gap ℏΩ inhibits the production of spin waves (generated by the inhomogeneous

_{S}*h*) and keeps the dynamics within the collective Dicke manifold: Flipping a single spin would reduce

_{i}*S*to (

*N*/2 − 1) and change the exchange energy, proportional to

*NJ*.

Figure 3 (A to E) shows the qualitative change in dynamical behavior as *J* crosses *J _{c}* for fixed

*t*) decays monotonically in time (Fig. 3A), but above the transition, S(

*t*) oscillates around a nonzero magnetization (Fig. 3, C to E). All observations can be reproduced by the same theoretical model shown in Fig. 2 (red lines) if

*55*). The TDP is seen in three observables (Fig. 3, F to H): first, by a departure from ungapped dynamics; second, by looking for a jump in S at 100 ms; and last, by a finite value of Ω.

The χ^{2} measure in Fig. 3F compares both data and calculations toS_{J=0}(*t*), the calculated time evolution for *J* = 0. The sharp increase near *NJ*/2π ≈ 10 Hz in both experiment and theory indicates the qualitative deviation of the dynamics from the paramagnetic phase. The magnetization at *t* = 100 ms (in Fig. 3G) also shows an increase near *NJ*/2π ≈ 10 Hz. Numerical solutions of the mean-field dynamics with thermal averaging at finite time (red band) and without thermal averaging in steady state (solid black line) agree well with the experimental data. The simplified all-to-all model (Eq. 3, dotted black line) agrees only qualitatively. The gap frequency Ω in Fig. 3H is found from a fitting function that uses a damped sinusoid for later times and S_{J=0} for early times (see the Supplementary Materials). By fitting the gap parameter to the analytic formula in Eq. 4, we extract a nonzero critical interaction strength *NJ _{c}*/2π = 7.8(1.1) Hz. Using the location of the step in χ

^{2}(Fig. 3G), we exclude time sequences with ∣

*a*∣ < 3

*a*

_{0}, where the oscillation frequency diverges. The excellent agreement of all three measures with theory based on Eq. 1 confirms that our quantum simulator probes the TDP in the collective Heisenberg model.

The significance of the observation is further clarified in Fig. 3 (G and H) by comparison to various approximation levels. Finite-time effects and thermal averaging play a minor role, validating our interpretation of S at sufficiently large *t* as the steady-state order parameter. Inhomogeneous coupling (*J _{ij}* ≠

*J*) plays a significant role for S but less so for Ω. Comparisons to the exact Lax vector analysis [insets in Fig. 3 (G and H)] show the close similarity of the observed TDP to the phase I–to–phase II transition in dynamical superconductors (

*28*–

*30*).

### Effective time reversal

Figure 4 describes a further set of simulations that probe the limits of validity in which our system is described by a spin-lattice model. Figure 4 (A and B) shows that S(*t* = 100 ms) decreases at sufficiently large *a* despite a larger gap. This is accompanied by a breakdown in microscopic reversibility, as seen by comparing to a sequence with a many-body reversal of the spin model (Fig. 4C), in which *4*, *47*, *48*). Signatures of time reversibility within the window ∣*a*∣≲ 20*a*_{0} are seen from the nearly *J*-independent dynamics of S in both the gapped and ungapped phases. The reversibility of demagnetization in our system in this regime is a significant validation of Eq. 1 since the many-body echo sequence does not reverse all terms (e.g., the spin-independent harmonic oscillator term) in the full Hamiltonian.

Two processes prevent full reversibility in our simulation: stray magnetic field gradients and collisional processes. These are quantified by introducing an empirical dephasing rate Γ(*a*) = Γ_{0} + (*a*/*a*_{0})^{2}γ to the transverse magnetization such that S^{X, Y}(*t*) → S^{X, Y}(*t*)*e*^{−Γ(a)t}. Figure 4 (B and C) compares data with calculations using a best-fit ^{−1} = 600 s. Here, Γ_{0} accounts for the single-particle mode-changing processes generated by magnetic field gradients, which are enhanced by a spin reversal (*1*, *40*). γ parametrizes mode-changing collisions that take place at a rate that increases quadratically with *a* and takes a value anticipated by kinetic theory for our experimental density, temperature, and polarization (see the Supplementary Materials). At larger ∣*a*∣, coherence is lost and a quantum picture becomes unnecessary, as shown by the success of a semiclassical picture to describe diffusive transverse demagnetization (*3*, *4*, *56*–*60*). Figure 4A shows an extended phase diagram, where it can be seen that trapped fermions simulate the collective Heisenberg model only in a restricted parameter window of many-body quantum coherence.

## DISCUSSION

In summary, we demonstrate the existence of a TDP in a neutral Fermi gas in a regime of reversible dynamics near the zero crossing of a Feshbach resonance. We outline a direct connection to nonequilibrium phases in the Richardson-Gaudin models for superconductivity, thereby extending experimental observations of TDP’s beyond prior manifestations in Josephson- and Ising-type systems (*19*, *34*–*38*). Moreover, the excellent agreement between spin-model calculations and a two-axis exploration of the dynamical phase diagram with >10^{4} spins provide experimental evidence of the scaling behavior and universal character of TDPs.

The collective nature of the dynamics observed here could protect many-body states in other systems of interest for applied quantum technologies. For example, gap protection would increase the coherence time in optical lattice clocks operated in the quantum degenerate regime (*61*). Furthermore, using the effective time-reversal capability demonstrated here, together with technical improvements to magnetic field stability and homogeneity, our system could provide a fruitful platform to measure out-of-time order correlations and scrambling of quantum information (*48*) or test spectroscopic protocols that use time reversal to relax the detection resolution required for spectroscopy beyond the standard quantum limit (*62*).

## MATERIALS AND METHODS

The neutral atomic sample was prepared using laser cooling, magnetic and optical trapping, and evaporative and sympathetic cooling in an atom-chip apparatus described previously (*63*, *64*). The hyperfine states used to encode spin information are the *F* = 9/2, *m _{F}* = −9/2 (for ∣↓〉), and

*m*= −7/2 (for ∣↑〉) states. At the beginning of a simulation sequence, the fermionic ensemble had a typical temperature of several hundred nanokelvin, with data taken in the range of

_{F}*T*~ 0.3 to 0.5

*E*

_{F}/

*k*

_{B}, where

The magnetic field and its gradients were controlled using a combination of microfabricated wires on the atom chip ≈200 μm from the atoms and macroscopic coils external to the vacuum system. The field was calibrated through rf spectroscopy of the ∣↓〉 to ∣↑〉 transition. During a typical experimental run, drifts are ≲1 μT.

Taking advantage of the symmetry of the many-body dynamics with *a*, we determined the magnetic field corresponding to *a* = 0 by observing the magnetization after *t* = 100 ms at various magnetic fields and then fitting the resulting profile to find the center. The mean result of 10 such trials taken across 6 months results in *B*_{zc} = 20.907(1) mT, which is an improvement in uncertainty by an order of magnitude over previous measurements (see the Supplementary Materials for data and literature comparison).

Magnetic field gradients, which lead to periodic oscillation of the spin clouds (see the Supplementary Materials), were measured by displacing the trap center and repeating rf spectroscopy. In the optimal configuration, gradients are {13(1),12(1),2(5)} μT/m. The differential displacements Δ* _{i}* resulting from these gradients are small compared to the harmonic oscillator length

*/*

_{i}*a*= {2.1,0.4,0.1} × 10

_{i}^{−2}. In this Δ

*≪*

_{i}*a*regime, the more general XXZ spin model (

_{i}*41*) reduces to the Heisenberg model used here.

The effective axial field *h _{i}* in the Heisenberg Hamiltonian is the potential energy differential between the ∣↓〉 and ∣↑〉 states, as sampled by the occupied motional eigenstates. Since one can always subtract a constant potential term, dynamics depend only on the inhomogeneity in

*h*, quantified through its root mean square spread

_{i}*h*. The leading-order magnetic field contribution is curvature in real space. Direct spectroscopic data can bound this to ≲500 μT/m

_{i}^{2}, which is compatible with the 200 μT/m

^{2}best-fit curvature in the model. The optical contribution to

*g*(δ

_{F}_{FS}/3δ), where P is the polarization of the light,

*g*is the

_{F}*g*-factor of the hyperfine sublevels, δ is the frequency detuning, and δ

_{FS}is the fine-structure splitting of the electronic excited state. For a full range of polarization P ∈ [ −1,1] of the optical dipole trapping beam propagating along the

*z*axis, the resulting VLS should modify the trap frequency along

*x*by roughly ±0.06 Hz or ±0.015%. This control over

The magnetization 2S/*N* was determined using repeated measurements and a maximum-likelihood estimator, as follows. Each experimental image measures *N*_{↑}, *N*_{↓} at the end of a Ramsey sequence. The fraction of spin-↑, *f* = *N*_{↑}/*N*, relates to the collective spin as *f* = 1/2 + (*N*)cosθ + (*N*)sinθ, where θ is the phase lag of the second π/2 pulse. However, the typical evolution times (100 ms) exceed the clock coherence time (1.5 ms) so that the accumulated phase ϕ randomizes the orientation of the transverse spin, *f* = 1/2 + (S/*N*) cos (θ′), where θ′ = θ − φ is a random phase. To reconstruct the offset *O* and amplitude *A* of a Ramsey fringe *F* = *O* + *A*cos (θ′) from a set of fractions *f _{i}* acquired in several experimental runs with the same conditions, we assumed a probability distribution

*n*fraction measurements,

*A*and fringe offset

*O*, as well as confidence intervals. Where

*O*falls outside the range of 0.50 ± 0.05, we discarded the data; otherwise, the peak-to-peak amplitude 2

*A*was taken as the best estimate of 2S/

*N*.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/5/8/eaax1568/DC1

Fig. S1. Order parameters predicted by Lax vector analysis in a 1D system.

Fig. S2. Approximate form of the Lax vector in a 3D system.

Fig. S3. Effective mean-field potential.

Fig. S4. Magnetization dynamics for noninteracting particles.

Fig. S5. Determination of the Feshbach zero crossing.

Fig. S6. Spin-echo amplitude near the Feshbach zero crossing.

Fig. S7. Fits to time series.

Fig. S8. Equilibrium scattering rate versus temperature.

Fig. S9. Nonequilibrium scattering rate.

Table S1. Determination of ^{40}K Feshbach resonance parameters.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank A. Koller and C. Luciuk for early work on this project and V. Gurarie, B. Lev, J. Thompson, M. Foster, and D. Stamper-Kurn for discussions.

**Funding:**This work is supported by NSERC; the Air Force Office of Scientific Research grants FA9550-13-1-0063 and FA9550-18-1-0319 and its Multidisciplinary University Research Initiative grant (MURI); the Army Research Office grants W911NF-15-1-0603, W911NF-16-1-0576, and W911NF-19-1-0210; the Defense Advanced Research Projects Agency (DARPA); the NSF grant PHY1820885; JILA-NSF grant PFC-173400; and the National Institute of Standards and Technology.

**Author contributions:**The work was conceived by A.M.R., J.H.T., and S.T. Experiments were performed by S.S., B.A.O., H.S., K.G.J., and S.T. Data were analyzed by S.S., P.H., and B.A.O. Theoretical models and simulation were done by P.H., J.M., J.H.T., and A.M.R. All authors contributed to manuscript preparation.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials.The datasets generated and analyzed during the current study are available from the corresponding authors upon reasonable request.

- Copyright © 2019 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).