Relaxation dynamics of an isolated large-spin Fermi gas far from equilibrium

A fundamental question in many-body physics is how closed quantum systems reach equilibrium. We address this question experimentally and theoretically in an ultracold large-spin Fermi gas where we find a complex interplay between internal and motional degrees of freedom. The fermions are initially prepared far from equilibrium with only a few spin states occupied. The subsequent dynamics leading to redistribution among all spin states is observed experimentally and simulated theoretically using a kinetic Boltzmann equation with full spin coherence. The latter is derived microscopically and provides good agreement with experimental data without any free parameters. We identify several collisional processes, which occur on different time scales. By varying density and magnetic field, we control the relaxation dynamics and are able to continuously tune the character of a subset of spin states from an open to a closed system.


I. INTRODUCTION
The relaxation of closed quantum systems towards equilibrium is a fundamental problem in many-body physics. It is particularly challenging to fully understand this macroscopic process on the basis of microscopic properties [1][2][3][4]. Here, ultracold atomic quantum gases provide an exceptional experimental platform due to the nearly perfect isolation from their environment and the excellent control on a microscopic level. In particular, the possibility to prepare well-defined states far from equilibrium, as well as widely tunable Hamiltonians, has recently attracted a lot of attention, e.g. prethermalization [5][6][7], relaxation in strongly interacting lattice systems [8][9][10][11][12] and the interplay between thermal and condensate fraction of multicomponent bosons have been studied [13].
Spinor quantum gases are of particular interest since the spin offers an additional degree of freedom, giving rise to complex dynamics involving different relaxation processes on different time scales. Ultracold bosonic quantum gases have been intensively studied and exhibit a rich variety of effects such as texture formation and spin dynamics in spinor Bose-Einstein condensates [14][15][16][17][18], which can be well described theoretically using a multicomponent Gross-Pitaevskii-Equation [19][20][21]. Recently, collective spin dynamics was also observed in a thermal Bose gas [22]. Fermions, in contrast, are governed by * Electronic address:ulrich.ebling@icfo.es † These two authors contributed equally to this paper Pauli blocking and reveal a different behavior. So far, most experiments studied spin 1/2 fermions e.g. the BEC-BCS crossover [23,24], thermodynamic and transport properties [25,26], collective excitations [27,28] and magnetic ordering [29][30][31][32]. Spin-related phenomena in multi-component Fermi gases (F > 1/2) have been investigated only recently, in the context of spin-mixing dynamics on individual sites of an optical lattice [33], or collective coherent excitations in a trapped system [34][35][36]. The latter has been proven to be well described within a Boltzmann equation. The question, how a largespin fermionic many-body system reaches an equilibrium state via relaxation involving spin and spatial degree of freedom, has not been addressed.
In this paper, we study the relaxation dynamics of a trapped fermionic quantum gas of 40 K atoms with large spin of F = 9/2. Starting from an initial mixture with only a few spin states occupied, we observe a rich relaxation dynamics leading to a redistribution of the atoms among all available spin states. We study the intermediate regime between the collisionless and the hydrodynamic limit. In the collisionless limit, interactions are weak and can be taken into account on a mean-field level, while the hydrodynamic limit is characterized by stronger interactions which ensure local equilibrium. The dynamics in this intermediate regime is governed by different processes on very different time scales. We identify these processes by deriving a Boltzmann equation from the microscopic Hamiltonian of the large-spin system. This approach describes the time evolution of the system on the level of single particles in contact with the bath of the many-body system [37,38]. This corresponds to the in- tuitive expectation that the system acts as a bath for its own subsystems. We present a detailed comparison between numerical simulations and experimental data and find good agreement. Our analysis includes the dependence of the relaxation on density as well as on magnetic field. Whereas a higher density enhances the spin relaxation, we find a suppression of spin-changing processes at large magnetic fields due to the quadratic Zeeman shift. The latter effect can be used to control the loss of particles from the subsystem defined by the initially occupied spin components into the initially empty spin states. Generally, we observe that the relaxation within a subset of spin states, driven by incoherent spin-conserving collisions, happens on a much faster time scale than the redistribution among the spin components due to spin-changing incoherent collisions. The reason is that the spin-changing collisions are driven by the relatively small part of the interactions that breaks the SU (10) symmetry between the spin states. Thus, we encounter a situation similar to prethermalization [7], where first a prethermal state is reached, approximately conserving the initial occupations of the single spin states, before the redistribution among all spin states due to slight symmetry breaking sets in. This separation of time scales also allows us to monitor the increase of (effective) temperature within the subsystem of the initially populated spin states, as it is caused by dissipation into empty spin states.

II. RELAXATION PROCESSES IN A LARGE-SPIN SYSTEM
We perform measurements in a quantum degenerate gas of fermionic 40 K, which has total spin F = 9/2 in its hyperfine ground state, yielding ten spin states m = −9/2 . . . + 9/2, as depicted in FIG. 1(a). We prepare an atomic sample with two spin states occupied (see Appendix A 1 for details), confined in a spinindependent dipole trap. Due to the broken SU (N ) symmetry in 40 K resulting from spin-dependent scattering lengths, spin-changing collisions can occur. A microscopic collision process is depicted in FIG. 1(b): Two particles collide and exchange both spin m and momentum k: (|m 1 , k 1 + |m 2 , k 2 → |m 1 , k 1 + |m 2 , k 2 ). The total spin S, the total magnetization M = m 1 + m 2 as well as the total momentum k 1 + k 2 have to be conserved in this process. As a particular fermionic feature, the Pauli exclusion principle has to be obeyed, i.e. m 1 = m 2 and m 1 = m 2 . The interplay between the differential quadratic Zeeman energy ∝ m 2 1 + m 2 2 − m 2 1 − m 2 2 and interaction energy determines whether spin-changing collisions are likely or suppressed. In the presence of spinchanging collisions, the atoms will in general relax into a steady state with population in all ten spin states. Hence, preparation of an initial non-equilibrium state with only a few spin states populated will lead to complex dynamics, in which more and more spin states are gradually occupied (see FIG. 1(c)). It is a compelling question how the system relaxes towards a steady state.
In FIG. 2 we show exemplarily an experimentally obtained time evolution of the spin occupations in our system. Here, the initial spin configuration is a superposition of all ten spin states created by rotating a mixture of the states m = ±1/2 using rf-pulses [36]. We can clearly identify three different processes occurring on three different time scales: (i) We observe coherent spin-changing oscillations with a periodicity on the order of hundred ms. (ii) These oscillations are damped with a rate on the order of several hundred ms. (iii) Beyond this, we observe a slow redistribution among the ten spin states on a much longer time scale on the order of tens of seconds.
In the following, we derive a Boltzmann equation, which reproduces the experimentally observed effects and The abovementioned Boltzmann equation includes all these collision processes and captures the nonequilibrium dynamics in a general fashion, applicable to trapped weakly-interacting gases with arbitrary spin. In this approach, the single-particle dynamics is treated as an open system in contact with the environment represented by all the other particles. An approach for deriving a Boltzmann equation was applied successfully to the description of spin dynamics in liquid Hydrogen and Helium [39][40][41][42][43][44][45] and later in spin 1/2 Fermi gases [46,47]. In this paper we generalize this approach to one-and three-dimensional systems with large spin, accounting for the quadratic Zeeman effect (QZE). In general, a kinetic equation or Boltzmann equation is used to describe the time evolution of the single-particle density matrixρ. It has the form Here,Ĥ 0 denotes the single-particle Hamiltonian which contains the kinetic energy, the harmonic trapping potential and the quadratic Zeeman splitting Q induced by a homogeneous magnetic field. The term I coll [ρ] on the right hand side of (1) is called the collision term and is derived from two-particle contact interaction. Due to total spin conservation, collisions are best described in the basis of total spin |S, M = m1,m2 |m 1 , m 2 m 1 m 2 |SM with the short notation for Clebsch-Gordan coefficients m 1 m 2 |SM ≡ F, m 1 ; F, m 2 |S, M , which we use throughout this paper. In general, scattering in each channel of total spin S depends on a different s-wave scattering length a S . Due to antisymmetrization of the total wave function, s-wave scattering with odd S is forbidden. Thus for 40 K there are five different scattering lengths present for S = 0, 2, 4, 6, 8 [33]. In each collision channel defined by |S, M , particles interact with a contact interaction of strength g 3D S = 4π 2 M a S , which is used in all 3D calculations. We also consider a 1D system, in which the motion in two transversal directions is frozen out completely by a tight trapping potential (characterized by radial frequencies ω x,y ) such that the effective 1D contact interaction parameter is given by g 1D S = 2 √ ω y ω z a S (see appendix A 1 for details). The notation g S ≡ g 1D S for this quantity is used throughout this paper. In 40 K, the a S range from 120 to 170 Bohr radii.
We obtain an explicit expression for the collision term in (1) using the method originally developed by Lhuillier and Laloë [43,44,46] for transport properties in Helium. In this approach, collisions are treated as single "atomic beam" experiments, where the colliding particles are assumed to be uncorrelated before and after a collision, reminiscent of Boltzmann's original molecular chaos hypothesis, but the scattering process itself is treated on a full quantum level. This approximation is valid for dilute gases where the mean time between collisions is long and the particle number is large. In this regime, binary collisions can be described by the T -matrix, which connects the two-body density matrices before and after a collision. Subsequently, the description is reduced to a single-particle level by tracing out the second particle, similar to tracing out a thermal bath in studies of collisional decoherence [37,38].
We calculate the kinetic equation (1) in its phase-space representation, where the single-particle density matrix ρ mn (x, x ) is expressed by the Wigner function Note that we performed the transformation only with respect to the spatial degrees of freedom. With respect to spin, denoted by the indices, it retains the form of a single-particle density matrix. The derivation is carried out in detail in section IV and involves a semiclassical gradient expansion of the Wigner function in position and momentum space leading to an equation in matrix form given by where the collision integral reads Here, we define the coupling constants and denote energy shifts induced by the quadratic Zeeman splitting as The infrared cutoffs are given by 1 Equation (4) contains several terms, each describing a different dynamical process. The free particle motion in the trap is described by The leading interaction term appears in the commutator [·, ·]. The commutator drives coherent spin dynamics through the interplay of the quadratic Zeeman effect and a spindependent mean-field potential resulting from forward scattering: a function of the density N (x) = dpW (x, p). In our large-spin system, described by several scattering lengths a S , it is helpful to decompose the mean-field potential (8) into two contributions. The first contribution is symmetric with respect to all N spin states and proportional to the mean scattering length. It conserves the occupations of the different spin components. The second term contains that part of the interactions that breaks the SU (N ) symmetry between the spin states and describes spinchanging processes. It depends on differences of scattering lengths only and is, thus, typically much smaller than the symmetric term. The commutator in Eq. (4) vanishes unless the Wigner function possesses off-diagonal elements indicating spin coherence. Moreover, the symmetric spin-conserving mean-field interactions can only contribute if the Wigner function describes an inhomogeneous spin state.
The mean-field potential also appears in the anticommutator {·, ·}. This term results from the subleading order of the semiclassical gradient expansion (where also the spin-independent trapping potential appears) and it is generally smaller than the commutator. It describes spin-dependent forces, that modify the kinetics in the trap. The collision integral (5) describes effects beyond mean-field that result from non-forward scattering and generates a dynamics that appears incoherent on the  (4) and (dots) from a 1D version of the single-mode approximation (9). Axial trapping frequency is ωx = 2π × 84 Hz and radial frequencies are ωy,z = 2π × 47 kHz, particle number N = 100 per tube at temperature T /TF = 0.2 and magnetic field B = 0.12 G. Inset: The system approaches a steady state for longer times. (b) Experimental data (circles) in a 3D configuration compared to calculations (lines) in single-mode approximation (9), ω = 2π × (33, 33, 137) Hz, N = 1.3 × 10 5 and T /TF = 0.15 at B = 0.34 G.
level of a single-particle description. It is quadratic in the scattering lengths. Again we have to distinguish between SU (N ) symmetric spin-conserving collision processes on the one hand and spin-changing collisions on the other. The latter processes are described by those terms for which the quadratic Zeeman shifts ∆ are nonzero and they are much smaller than the former, since they depend on the relatively small differences between the scattering lengths only.
The 1D equation (4) (4), which is linear in differences of scattering lengths and describes forward collisions. (ii) The damping of coherent phenomena arises from the spin-conserving part of the collision integral (5), which is quadratic in the scattering lengths. We have checked that spatial dephasing is not responsible as it is suppressed by the dynamically induced long-range nature of mean-field interactions induced by the rapid particle motion in the trap [36,48]. (iii) The long-term relaxation originates from the spin-changing non-forward collisions in the collision term, quadratic in differences of scattering lengths. Spin-conserving forward scattering does not play a role in the dynamics, it only has a noticeable effect, if spatial symmetry is broken by a magnetic field gradient, as in studies of spin-waves [28,35,46], which is not the case in our setup. The anticommutator in Eq. (4) leads to a mean-field driven correction to the trapping potential which is however negligible in the experiments considered here.
The collision integral (5) enables us to determine whether our system is in the collisionless, hydrodynamic or an intermediate regime. The average collision time in the 3D setup τ 3D ∼ (4πa 2 n p v T ) −1 [47], with the relevant scattering length a, peak density n p and thermal velocity v T = k B T /M , ranges from ∼ 30 ms to ∼ 150 ms for spin-conserving collisions and ∼ 3 s to ∼ 15 s in the spin-changing case. Compared to the average trapping frequency ofω = (ω x ω y ω z ) 1/3 ≈ 2π × 58 Hz, we obtain values forωτ 3D between 11 and 55 for the spin-conserving collisions and between 1100 and 5500 in the spin-changing case. The lowest and highest values of ωτ 3D are reached for the lowest and highest densities shown in FIG. 5 respectively. This means we may approach the hydrodynamic regime, where the collision rate is larger thanω and local equilibrium can be established. On the other hand our system becomes almost collisionless regarding the spin-changing collisions. Generally we are in an intermediate regime. In the 1D case, collision times τ 1D ∼ (n p ω y ω z a 2 /v T ) −1 are on the order of 0.35 ms and 35 ms respectively, meaning that ωτ 1D ∼ 0.2 and ωτ 1D ∼ 20, concerning spin-conserving and spinchanging collisions respectively. Hence with respect to the former, the system would be hydrodynamic. However, it is still in an intermediate regime regarding the redistribution of particles among the spin states driven by spin-changing collisions.

III. DISSIPATIVE REDISTRIBUTION OF SPIN OCCUPATIONS
In the following, we focus on the long-term spin relaxation shown in FIG. 2(b) and FIG. 3, while recent experiments have studied spin oscillations and their damping [36]. In order to restrict the dynamics to this process, we initially prepare a spin mixture consisting of only the spin states m = ±1/2 without coherences. In this case, the coherent oscillations driven by the commutator in Eq. (4) are absent and the Wigner function W mn remains diagonal at all times. In the following we investigate theoretically and experimentally this spin relaxation dynamics in 3D as well as 1D systems. For a direct comparison between theory and experiment, we realize a 1D system employing a deep 2D optical lattice, which confines the atoms into tight elongated tubes [8,49] as described in Appendix A 1. As shown in FIG. 4(a), the system gradually occupies all spin states and evolves towards a state of almost equal spin populations on a time scale of milliseconds. As a key result, we can well reproduce the experimentally observed dynamics using the full 1D Boltzmann equation without free parameters.
For harmonically trapped 3D systems, where all trap frequencies are about equal, we derive the full 3D version of equation (4) as well [see Eq. (28) in section IV]. However, numerical simulations of this equation are too demanding numerically. Nevertheless, the trap-induced motion of the particles is considerably faster than meanfield or relaxation dynamics, which averages the spatial dependence of the interaction via dynamically created long-range interactions. The Wigner function then approximately separates into a product W mn (x, p, t) ≈ M mn (t)·f 0 (x, p) [22,36,48]. The spatial part is assumed to be time independent and given by the initial equilib- We substitute this expression into the 3D kinetic equation (28) with the appropriate collision term. Hence, for negligible magnetic fields we find an equation for the matrix M mn (t) given by where The 3D coupling constants are given by Note that in the 1D case the single-mode approximation has a similar form given by T abcd For both single-mode approximations, the quadratic Zeeman shift has been neglected in the above equations, which are thus valid for small magnetic fields only (See Appendix G for full equations). In FIG. 4 we compare results from a single-mode approximation with experiments in a harmonically trapped Fermi gas yielding a surprisingly good agreement without free parameters. Note the qualitatively comparable behavior on different time scales of milliseconds for 1D and seconds for 3D. On the contrary the damping of the coherent spin oscillations visible in FIG. 2(a) and FIG. 3 is not described by this approach. This results from the assumption for the single-mode approximation that it completely neglects the multi-mode character of the fermionic many-body system and thus cannot account for spatial redistribution via lateral scattering events.
The high degree of control over all crucial parameters allows for a detailed investigation of the spin redistribution. To obtain further insight into the relaxation mechanisms, we measure the relaxation dynamics (as exemplarily shown in FIG. 4) for different densities while keeping T /T F constant. With higher density the collision rate increases and the relaxation process accelerates, as shown in FIG. 5. The measured rates correspond to the redistribution of the initially populated components m = ±1/2 into m = ±3/2, ±5/2 and are well reproduced using the single-mode approximation (9). The rate of spin-changing collisions increases with increasing density in accordance with the density dependence of the integral λ (10). As a second important parameter of the system, we investigate the influence of the magnetic field on the relaxation process. As the Zeeman energy of an atom pair changes during a spin-changing collision, a strong magnetic field suppresses this process by increasing the energy difference between the initial and final spin configuration. In FIG. 6(a), we depict the experimentally obtained populations of the spin components after 2 s as a function of the magnetic field strength (see Appendix A 2 for details) and compare them to single-mode [ FIG. 6(a)] and 1D calculations [FIG. 6(b)] after 2 ms. In both cases, the general behavior is very similar and shows a suppression of spin-changing collisions for large magnetic fields. Spin configurations with high values of |m| are energetically significantly separated from the initially populated m = ±1/2 and are only occupied at very low field strengths. By changing the magnetic field we can thus tune the magnitude of spin-changing collisions relative to the unaffected spin-conserving collisions up to a complete suppression. This gives us the possibility to view the m = ±1/2 subsystem as a dissipative twocomponent Fermi gas with a tunable loss mechanism.
We have further investigated the time evolution of the temperature of this subsystem exposed to losses induced by these spin-and momentum changing collisions. We compare two experiments: On the one hand, we perform an experiment at a high magnetic field (B = 7.6 G), where spin-changing collisions are suppressed. On the other hand, we perform a second experiment at a low magnetic field (B = 0.12 G) with strong spin relaxation.
In both cases the system is in the hydrodynamic limit with respect to external degrees of freedom due to the comparatively large spin-conserving interactions. Hence we make the assumption, that at each time the subsystems are close to an intermediate equilibrium state with a well-defined temperature. On the other hand, spinchanging collisions are at least two orders of magnitude weaker and very slowly change the particle number in each subsystem. This situation is reminiscent of prethermalization, first a "prethermal" state is reached under the assumption of conserved quantities, which on a much longer time scale are actually not fully conserved due to a "slightly broken symmetry", leading eventually to full thermalization [7]. Here, the role of the nearly conserved quantities is played by the occupation numbers of the ten spin states, which change only on a very long time scale. We measure the time evolution of the temperature of the initially populated m = ±1/2 components and compare the temperatures for both cases described above. For large magnetic fields we observe a small heating rate, which we mainly attribute to inelastic photon scattering. However, at low magnetic fields, the heating rate is significantly increased. In FIG. 7(a), we plot the temperature difference to extract the heating contributions solely generated by spin-changing collisions. This additional increase in temperature is due to hole creation in the Fermi sea [50] by scattering into the unoccupied spin states. We initially prepare a very cold two-component Fermi sea, with only few unoccupied trap levels below the Fermi energy. Losses through spin-changing collisions "perforate" this Fermi sea with holes, such that the experimentally obtained temperature increases. Numerical simulations using the Boltzmann equation (4) confirm the experimentally observed heating induced by redistribution [see FIG. 7(b)].

IV. MICROSCOPIC DERIVATION OF A LARGE-SPIN BOLTZMANN EQUATION
In this section, we describe the derivation of the 1D Boltzmann equation (4) in more detail. The reader not immediately interested in the details may skip this section and proceed directly to the conclusions. We follow previous work on the theoretical description of spinpolarized systems of H or He [43], called the Lhuillier-Laloë transport equation. We extend it to describe a 1D system with large spin, several scattering channels and a quadratic Zeeman effect. We consider this approach to be suitable for our purpose for a couple of reasons. The entire equation is derived from a microscopic collisional approach, so the collision term is not based on phenomenological assumptions. We avoid the use of a relaxation time approximation, widely used to describe damping in bosonic and fermionic systems [51][52][53], where the collision term is approximated by I coll [W ] = Weq−W τ , with a relaxation time τ . The reason is that for a multicomponent system determining the equilibrium state W eq is very challenging. Also due to the interplay of many different scattering lengths we expect not one but many different relaxation times for each spin component. Our approach allows to better understand the relaxation process itself, rather than merely its effect on other pro-cesses. Furthermore, from a technical point of view, our approach remains quadratic in the Wigner function so it can be numerically simulated using the same standard techniques as the collisionless case [28,47,48].
The idea behind the approach of Lhuillier-Laloë is to interpret the collision integral as the change rate of the state of a single particleρ →ρ due to binary collisions Here ∆t is the elapsed time interval, which is short compared to any relevant macroscopic dynamics of the system, but nevertheless is longer than the duration of a single collision, which is thus considered to be effectively instantaneous. This quantity will drop out and not appear in the final kinetic equation. With this in mind, we treat collisions in the asymptotic limit, where they are described by the Heisenberg S-matrix. It relates the twobody density matrix of both scattering particles before a collisionρ(1, 2) with the one after a collisionρ (1,2) .
Here (1,2) label the quantum numbers of particles 1 and 2 in first quantization. We obtain In order to arrive at a single-particle description we trace out particle 2 later. We next assume that particles involved in a collision are uncorrelated,ρ(1, 2) =ρ(1) ⊗ ρ(2), both before and after the collision, an assumption justified for a system with a large number of particles. This assumption in fact corresponds to Boltzmann's original molecular chaos hypothesis (Stosszahlansatz). For the desired single-particle density matrices before and after a collision we obtain where we account for the indistinguishability of particles with the operator P ex exchanging the quantum numbers of particles 1 and 2. Due to fermionic statistics it comes with a minus sign. This ansatz yields the following expression for the collision integral The S-matrix is related to the T -matrix viaŜ = 1 1−2πiT such that Eq. (17) becomeŝ This expression then has to be evaluated in the phasespace representation. Before performing the trace operation, we compute the two-body Wigner transform of the expression in braces in Eq. (18). This is a very lengthy exercise we demonstrate in detail in Appendix B, as well as the subsequent trace, shown in Appendix C. In the course of these calculations, we require the elements of the Smatrix. In the center-of-mass system they are given by where k, k denote the incoming and outgoing wavevectors of the particles, and m, n; m , n incoming and outgoing spins respectively. The delta-function assures energy conservation, k = 2 k 2 2µ denotes kinetic energy with reduced mass µ = M 2 and Q abcd ≡ Q(a 2 +b 2 −c 2 −d 2 ) the shift in the quadratic Zeeman energy induced by a spin-changing collision. The on-shell T -matrixT (k, k ) depends formally on the relative wave-vectors k, k , but for our case of s-wave scattering, the dependence is only on the modulus. As they are related by energy conservation |k | = k 2 + Q abcd , effectively the dependence is only on k or k . The QZE-shift vanishes for spinconserving collisions, hence it is absent in the spin 1/2 case and has the effect that scattering processes with a large Q abcd are suppressed, because the T -matrix decays ∼ 1/|k | for large |k |. For 40 K, the splitting is given by the quadratic part of the Breit-Rabi-Formula [54], , with the Bohr magneton µ B , nuclear and electronic g-factors g I , g J and hyperfine structure coefficient a hfs [55].
To account for the spin-dependent interactions we separate the T -matrix into channels of total spin S and magnetization M and obtain elements 1 : a, 2 : b|T (k, k )|1 : c, 2 : d ≡ T abcd (k, k ) For a 1D system, the expression for a T -matrix in the channel with a coupling constant g S is In each scattering channel we perform a low energy expansion in the coupling constant g S → 0 up to second order, to maintain the unitarity of the S-matrix. Since for the 1D case, the expansion in powers of g S is accompanied by factors of (k ) −1 , we artificially create a singularity in the imaginary part of the T -matrix in this procedure. We remedy this problem by choosing a cutoff |k | < M g S 2 2 , which is the distance between zero and the maximum of the imaginary part of the T -matrix (see Appendix D for more details). This step is unnecessary in the 3D case discussed by Lhuillier and Laloë [43,44,46]. The result is then given by or respectively The leading terms linear in g S correspond to forward scattering, the quadratic terms to backward (lateral in higher dimensions) scattering processes. We do not explicitly denote the spin-dependent cutoff in equations (5), (26) and (27), where it is used in the integrals over q.
The expansion of the T -matrix is performed in addition to a semiclassical gradient expansion of the Wigner function (see Appendix E for details) to first order in the linear terms and to zero order in the quadratic expressions. During this procedure we encounter squares of delta functions whose interpretation is described in Appendix F. Finally we obtain a collisional integral, consisting of three parts I coll mn (x, p) = I mf mn (x, p) + I T mn (x, p) + I T 2 mn (x, p), (24) where I mf is linear in a S and contains the forward scattering part of collisions leading to phase-shifts, It coincides with the interaction term obtained from the treatment of the same problem on a simpler mean-field level [28,48]. Due to its effect as a non-linear modification of the trap and magnetic field, in Eq. (4) we have separated this term from the collisional integral and added it to the free motion of the particles in the external fields. The quadratic terms, which form Eq. (5) contain backward scattering, including momentum exchange between particles. They appear as dissipation on the single-particle level and are given by The integration domain cutoffs around q = 0, coming from Eq. (23) are given by 1 The corresponding result for a full 3D calculation (see [46] for the case of spin 1/2) reads The mean-field potential is given by V mf mn (r) = 2 dp ab U mnab W ab (r, p) and the collision term reads where p = e Ω q 2 + ∆ mnlabcd and e Ω denotes the unit vector corresponding to solid angle dΩ.
A physical interpretation of this expression can be obtained by looking at the origin of the individual terms. The upper two lines of Eq. (29) originate from the second order of the expansion of the T -matrix (23), describing the intensity shift in the forward scattered wave [43,46], while the first order only describes a phase shift and appears in the mean-field terms in (28). The bottom line of (29) contains all lateral scattering processes, hence the explicit angular dependence. In our formalism, the coupling constants U ,Ũ include all particle indistinguishability and exchange contributions, which are discussed in greater detail in [43].

V. CONCLUSION
We have presented a novel approach to study relaxation dynamics in a closed quantum system, exploiting the unique properties of a large-spin Fermi sea. For this system, we have derived a multicomponent kinetic equation without phenomenological assumptions nor prior knowledge of the equilibrium state. As a key result, we find that this approach is well suited for the quantitative description of weakly interacting fermionic manybody systems with large spin. Both, the comparison of numerical simulations with full spatial resolution to a 1D experiment as well as the comparison of a simplified single-mode approximation to a 3D experiment yield a very good agreement without free parameters. We identify different collisional processes on different time scales and identify spin relaxation as the slowest dynamical process of the system. A variation of the density and the geometry of the system changes the respective spin relaxation rates by several orders of magnitude, ranging from a few milliseconds to several seconds. By tuning the magnetic field, we can precisely control the coupling strengths of individual collision channels, allowing to tune the character of a subsystem of two spin components within the large-spin Fermi sea continuously from an open to a closed system. The spin relaxation manifests itself in a perforation of the Fermi sea accompanied with a temperature increase.
Our results broaden the understanding of many-body relaxation dynamics. In particular, the fermionic char-acter of the system underlines its model character for various systems in nature. The possibility to monitor different spin components individually allows to employ the large-spin Fermi sea for novel studies of decoherence and relaxation processes in quantum many-body systems. Furthermore, spin relaxation dynamics might play an important role for proposed fermionic large-spin phenomena, e.g. quantum-chromodynamic-like color superfluidity or large-spin texture formation. We sympathetically cool spin-polarized 40 K atoms in the state |F = 9/2, m = 9/2 down to a temperature of typically 0.1 T F in a magnetic trap, using bosonic 87 Rb as a buffer gas. Subsequently we transfer the atoms into a crossed circular-elliptical optical dipole trap operated at a wavelength of λ = 812 nm. Using radio-frequency (rf) pulses and rf-sweeps, we create a spin mixture, which we evaporate to quantum degeneracy by lowering the power of the dipole trap exponentially in 2 s. This results in a sample with particle numbers of the order of N ∼ 10 5 at temperatures of T = 0.1 − 0.2 T F . After the evaporation we compress the trap again to avoid particle loss during the experiments, realizing typical trapping frequencies of ω = 2π × (33, 33, 137) Hz. By varying the evaporation sequence and including additional waiting times, we can control the initial temperature and particle number independently in the same trap geometry. This allows us to modify the density while keeping T /T F approximately constant. Typically, a balanced mixture of atoms in spin states m = ±1/2 is used as initial state throughout this paper. To study the spin-changing dynamics, we switch the magnetic field to low values. In FIG. 2, beyond this, a coherent superposition of several spin states is prepared by applying subsequently a rf-pulse at low magnetic field corresponding to a spin rotation of θ = 0.44 (see [36] for more details). The 1D configuration used in FIG. 4  (a) is realized by adiabatically ramping up a 2D optical square lattice over 150 ms. The lattice is created by two orthogonal retro-reflected laser beams at wavelength λ = 1030 nm with a 1/e 2 radius of 200 µm detuned with respect to each other by several tens of megahertz. The lattice depth is 25 λ . This creates an array of 1D tubes, where a single tube can be described as a harmonically trapped system with frequencies ω x = 2π × 84 Hz and ω y,z = 2π × 47 kHz. With a particle number of N ≈ 100 and E F = 2π × 37 kHz, the radial trapping frequencies fulfill ω y,z > E F and at a temperature k B T = 0.2 E F , we can neglect a possible population of excited radial modes, hence we create a true 1D system. The extension of the radial ground state is around 1378 Bohr radii and thus one order of magnitude larger than the scattering lengths [33]. We thus neglect the possibility of a confinement-induced resonance [56], and use the effective coupling constants g 1D S = 2 √ ω y ω z a S .

Measurement
The relative populations of spin components are measured as follows: We release the atoms from the trap in an inhomogeneous magnetic field to separate the spin components during a time-of-flight expansion of typically 18.5 ms. We count the number of atoms in each spin component with resonant absorption imaging. For comparison, we measure the total number of particles as well as the temperature independently without the Stern-Gerlach field to avoid distortions of the particle cloud during the time-of-flight. The numbers given in this paper correspond to the initial temperature and particle number. In FIG. 7(a), in order to extract the change in temperature over time, we determine the temperature only in one spin component, circumventing deviations associated with the imbalance of the spin mixture. For instance, to measure the temperature in m = 1/2 we apply a sequence of linearly polarized microwave pulses with a duration of 50 µs to transfer all significantly occupied spin components m = 1/2 into the F = 7/2 hyperfine manifold of 40 K. In the other hyperfine manifold the atoms are not resonant with the detection light and are thus obscured during the absorption imaging process. and W T 2 ijmn (r, R, p, P ) = π dκ dk 1 dk 2 dr e iκr e i(k2−k1)r term dominate and we perform a Taylor expansion for the spatial coordinate (E1) therefore the expansion of the product of Wigner functions in (C2) reads Together with the expansion of the T -matrix above we must be careful to expand the collision term in two small parameters in a meaningful way. One small parameter is the coupling constant proportional to the s-wave scattering length. The other one is related to the gradient expansion. Its magnitude is determined by the Fermi or thermal wavelength compared to the variation of the Wigner function determined by the system size.
To maintain the unitarity of the S-matrix, we expand the T -matrix to second order. This means we will obtain terms linear in a S from I T ij (x, p) and quadratic terms from I T ij (x, p) and I T 2 ij (x, p). We expand the terms linear in a S up to first order in gradients and the terms quadratic in a S to zero order, keeping only the local term. This amount to a semi-classical approximation of the theory. In this case we substitute W ac (x 1 − r−r 2 , p 1 )W bd (x 1 − r+r 2 , p 2 ) ≈ W ac (x 1 , p 1 )W bd (x 1 , p 2 ) into (C2), which means that further delta functions dr e i(k2−k1)r = 2πδ(k 2 − k 1 ), dre iκr = 2πδ(κ) appear. We introduce renamed variables k ± → k, k 1 → k , x 1 , p 1 → x, p and p ± ≡ p − (k ± k ) the local parts of the collision integral become and Appendix F: Squares and products of delta functions In scattering theory, the square of a delta function of energy appears frequently, when terms quadratic in the T -matrix are involved. A well-known interpretation of this is [δ(E)] 2 ≈ ∆t 2π δ(E), where ∆t denotes the elapsed time interval, which is quasi-infinite when compared to the duration of a single scattering event but nevertheless short compared to other relevant dynamics, like relaxation or the trapping period. This approximation is obtained by using the Fourier representation of the delta function This can also be applied to products of the form δ( k − k )δ(k − k ) since We modify this approximation to take into account the shift Q in the quadratic Zeeman energy after a spinchanging collision. In our calculations two situations appear. In the first, coming from (F3), there is only one shift and we must be careful that only the delta-function with the shift comes from a T -matrix where we can approximate the integration area with the interval ∆t: In the second case, both delta-functions originate from the energy conservation of the T -matrix )u e i 2 (Q2−Q1)u ≈ ∆t 2π δ k − k + 1 2 (Q 1 + Q 2 ) sinc The time interval ∆t that appears in front cancels with the one introduced at the beginning (18) and for the sincfunction we assume it to be small such that sinc → 1.

Appendix G: Single-mode approximation with QZE and in 1D
Taking the quadratic Zeeman effect into account, the expressions for the single-mode approximation (9) become slightly more complicated. The equation of motion is now given by where the two now separate integrals λ (1,2) are spin-dependent and given by λ (1) abcd = 1 N dr dp dq q 2 + ∆ abcd f 0 (r, p)f 0 (r, p − q), and λ (2) mnlabcd = 1 N dr dp dq q 2 + ∆ mnlabcd f 0 (r, p)f 0 (r, p − q), respectively, where the infrared cutoff described in Appendix D must be employed. The single-mode approximation can also be applied to the 1D system, from the Boltzmann equation (4). Equation 1 N dx dp q 2 > 1,2 dq f 0 (r, p)f 0 (r, p − q) and λ (2) mnlabcd = 1 N dx dp q 2 > 3 dq f 0 (r, p)f 0 (r, p − q) A comparison of 1D single-mode results with the full 1D Boltzmann equation shows good agreement for pure spin relaxation as shown in FIG. 4. A further inclusion of the coherent oscillations described by the commutator in Eq. (4) into the singlemode equation shows that while the oscillations themselves are reproduced with high accuracy [36], the damping of coherent oscillations such as in FIG. 3 is not captured well. We attribute this to the fact that damping is driven by the much stronger spin-conserving collisions and affects more strongly the individual phase-space distributions of the spin states, which are taken to be constant in time in the single-mode approximation. Thus we consider it necessary to use the full 1D Boltzmann in these cases.