Absence of heating in a uniform Fermi gas created by periodic driving

Ultracold atomic gas provides a useful tool to explore many-body physics. One of the recent additions to this experimental toolbox is the Floquet engineering, where periodic modulation of the Hamiltonian allows the creation of effective potentials that do not exist otherwise. When subject to external modulations, however, generic interacting many-body systems absorb energy, thus posing a heating problem that may impair the usefulness of this method. For discrete systems with bounded local energy, an exponentially suppressed heating rate with the driving frequency has been observed previously, leaving the system in a prethermal state for exceedingly long durations. But for systems in continuous space, the situation remains unclear. Here we show that Floquet engineering can be employed to a strongly interacting degenerate Fermi gas held in a flat box-like potential without inducing excessive heating on experimentally relevant timescales. The driving eliminates the effect of a spin-dependent potential originating from a simultaneous magnetic levitation of two different spin states. We calculate the heating rate and obtain a power-law suppression with the drive frequency. To further test the many-body behavior of the driven gas, we measure both the pair-condensation fraction at unitarity and the contact parameter across the BEC-BCS crossover. At low driving frequencies, the condensate fraction is reduced by the time-dependent force, but at higher frequencies, it revives and attains an even higher value than without driving. Our results are promising for future exploration of exotic many-body phases of a bulk strongly-interacting Fermi gas with dynamically engineered Hamiltonians.


I. INTRODUCTION
The last two decades have witnessed tremendous advance in studying many-body problems with ultracold atomic gases [1].The vast majority of works have been done with static Hamiltonians.Adding periodic driving can generate effective Hamiltonians with completely different properties than the original one -an approach called Floquet engineering [2][3][4].For example, modulation of the barrier between potential wells renormalizes the tunneling rate [5][6][7][8] and can drive quantum phase transitions [9].Floquet engineering can also be used to create artificial gauge fields [10][11][12][13], and give rise to new phases which do not exist at equilibrium [14,15].
An inherent problem with externally driven systems is their tendency to heat up.Apart from integrable and many-body localized systems, generic interacting ensembles will absorb energy from the external force and eventually reach an "infinite temperature" where all states are equally populated [16][17][18].Nonetheless, recent theoretical works suggest that in discrete lattice systems, the energy absorption rate is generally exponentially small in the drive frequency over the energy of a local excitation [19][20][21][22][23][24].These predictions are supported by heating rate measurements done with bosons in a driven optical lattice [25].The exponential suppression of heating in discrete systems relies on the fact that the energy is locally bounded.However, in continuous systems (e.g., bulk quantum gases), there is no such bound.The perti- * Electronic address: yoavsagi@technion.ac.il nent question, in this case, is under what conditions one can obtain "cold" prethermal states exhibiting collective phenomena that are governed by an effective "Floquet engineered" Hamiltonian.
In this paper we address this question using a driven ultracold Fermi gas near unitarity, exhibiting high-T c fermionic superfluidity [26].The driving we apply is intended to create a uniform effective potential across the Fermi gas.Although in situ measurements [27][28][29][30] and spatial selection [31][32][33][34][35][36] can give access to quasihomogeneous observables, it is better to create a uniform gas from the outset.This is essential, for example, to study critical properties and avoid spurious phase separated states [37].Indeed, in recent years, uniform Bose [38,39] and Fermi [40,41] gases have been created in flat optical traps.These traps are formed by several shaped laser beams that create sharp repelling walls enclosing a dark volume.
A significant challenge is posed by the need to offset the gravitational potential, which leads to a substantial energy change in the trap.One obvious solution is to use a shaped optical potential to counter gravity [42].But generating such a potential, smooth on a nano-Kelvin scale, is a formidable task.A simpler approach taken in previous experiments is to use a magnetic field with an appropriate gradient.This works if all particles have approximately the same magnetic dipole moment.However our 40 K Fermi gas is a mixture of two hyperfine states with different (but not opposite) magnetic moments µ ↑ = µ ↓ .An appropriate magnetic field gradient counters the average gravitational potential, while leaving an opposite potential gradient on each of the two species (Fig. 1a).As a result, it is only possible to partially counteract the gravitational potential with a magnetic field gradient, Bz, set according to Eq.( 3).The total external potential Vext,s depends on the spin, s ∈ {↑, ↓}, and consequently, the density distribution of each spin is different and not uniform (color gradient in the top left figure).(b) By adding a resonant rf field that drives rapid spin rotations, we create an effective spin-independent potential, in which the gas becomes homogeneous.Importantly, the intrinsic manybody behavior of the gas is unchanged by this driving.
To counter the residual field gradient we apply an rf field that induces a rapid percession at a Rabi frequency Ω.In the rotating frame of the percessing spins, the static spin dependent potential gradient is translated to a periodic perturbation of frequency Ω.The rest of the interacting fermion Hamiltonian, including the flat spinindependent potential, is invariant to spin rotations and therefore unchanged in the rotating frame (Fig. 1b).Thus we achieve an ultracold uniform state of a spinbalanced gas of 40 K atoms that is useful if the periodic perturbation does not cause significant heating over experimental time scales.
We establish this property by measuring the paircondensate fraction (CF) at unitarity while applying a continuous driving.At low frequencies, driving impairs the gas conditions and reduces the CF.As the frequency increases further, the CF recovers and even surpasses its value without the driving.At high driving frequencies, we do not detect heating or excessive loss of atoms which can be attributed to the drive.Finally, we perform rf spectroscopy with a uniform gas in the BEC-BCS crossover regime, and extract the homogeneous contact parameter as a function of the interaction strength.
The structure of this paper is as follows.In section II we review the theoretical model for radio frequency driving and calculate the expected heating rate in presence of this drive.In section III, we describe the experimental setup and measurement sequence.The results are presented in section IV.We study with time-dependent in situ imaging the relaxation dynamics following the application of the driving field.The temperature of the uniform gas is probed by Raman spectroscopy.The manybody behavior of the uniform driven gas is studied with a pair-projection technique and rf spectroscopy.We study both the frequency dependence and the long time behavior of the driven gas.Section V concludes with a discussion and outlook.

II. THEORETICAL MODEL
Effective Hamiltonian.-Weconsider fermions in two possible spin states, denoted by |↓ and |↑ , where the energy of the latter is larger by ω 0 .The two particles are placed in an external potential and coupled by an rf field with a frequency ω rf .The Hamiltonian is a sum of three terms Ĥ = Ĥ0 + Ĥint + Ĥrf that account for the single-particle kinetic and potential energy ( Ĥ0 ), the interaction Hamiltonian ( Ĥint ) and the coupling to the external rf field ( Ĥrf ).In the frame rotating with the ↑ spin, they are given by [43,44] Ĥrf = 2 d 3 r Ωe iω0t e iω rf t + e −iω rf t Ψ † ↑ (r) Ψ↓ (r)+h.c., (1c) where Ω is the Rabi frequency of the rf field, V σ (r) is the external potential for spin σ, and Ψσ (r) are fermionic field operators obeying the anti-commutation relation { Ψσ (r), Ψ † σ (r )} = δ σσ δ(r − r ).Note that here we consider a general spin symmetric (and translationally invariant) two body interaction.The interaction potential u(r − r ) represents the microscopic interaction, rather than a low energy limit of the T -matrix (pseudopotential).In particular we will be interested in a unitary gas for which the T -matrix is explicitly energy dependent in the low energy limit.
In our experiment, the rf field is resonant with the bare energy difference ω rf = ω 0 , and ω 0 Ω.Thus, within the rotating wave approximation, the rf field is seen as a large magnetic field along S x , Ĥrf = 2 d 3 rΩ Ψ † ↑ (r) Ψ↓ (r) + h.c..The initial state is presumed to be an ultracold Fermi gas with equally populated spin components N ↑ = N ↓ , which in presence of the field precesses around the x axis at the Rabi frequency Ω.To clearly see the effect of the external rf field, we eliminate the field by a unitary transformation, Û = e i Ĥrf t , into a reference frame that rotates with the spins.
The Hamiltonian transforms in a simple way under Û .The kinetic energy and interaction Hamiltonians are both invariant under spin rotations and are therefore unchanged by the time-dependent transformation.As first noted by Zwierlein et al. [45], the invariance of contact interactions under rf rotations is the reason for the absence of a spectroscopic shift in the transition frequency between the spins.The external potential, on the other hand, can be decomposed into spin symmetric and antisymmetric parts where n(r) = n↑ (r) + n↓ (r) is the number density, Ŝz (r) = (n ↑ (r) − n↓ (r))/2 the spin density, Only the spin symmetric part of the potential, which couples to the number density, is invariant under the transformation Û and does not change with time.
In our experiment, the external potential is given by where the first term is the flat optical potential, the second term is the gravitational potential, and the last term describes the interaction of a spin with a magnetic moment µ σ with the external magnetic field, which is linear in height z.Thus, the spin symmetric part of the potential is . By tuning the magnetic field gradient to the value we obtain a flat total potential.The gravitational potential is thus eliminated by the static part of the Hamiltonian.
The spin anti-symmetric part of the external potential, which couples to the spin density, gives rise to a time-dependent coupling in the rotating frame Û .Specifically, the spin density rotates as Û Ŝz Û † = cos(Ωt) Ŝz − sin(Ωt) Ŝy , leading to a Zeeman field rotating at frequency Ω in the Ŝz , Ŝy plane.With a simple redefinition of the spin axes we write the residual time-dependent part of the Hamiltonian as a rotating field in the Ŝx , Ŝy plane Here Ŝ+ (r) = Ψ † ↑ (r) Ψ↓ (r), θ(z) is the Heaviside step function, and L is the height of the box potential.
, where we use Eq.( 3) to tune a flat static potential.
Heating rate.-Theheating rate due to irreversible transitions induced by the residual time-dependent term (4) in the effective potential can be calculated from the linear response of the system to this perturbation.For this purpose it is convenient to rewrite the timedependent perturbation in momentum space Linear response to the residual timedependent potential.In the diagram, the incoming and outgoing dashed lines represent the operators Ŝ+ q and Ŝ− q attached at the vertices with the couplings Vq and V−q respectively.The vertices connecting four solid lines represent the contact s-wave interactions.
where Ŝ+ q = d 3 r Ŝ+ (r)e −iq•r and v q is the spatial Fourier transform of the time-dependent effective potential The Gaussian fall off at large q is due to convolution with the optical resolution that limits the sharpness of the potential features.In our experiment the resolution λ 0 ≈ 3µm happens to be close to the Fermi wavelength λ F ≈ 2.6µm, so for simplicity we shall identify the two scales and use λ F also as the resolution limit.
The transition rate Γ q induced by the perturbation at wave-vector q is directly related to the spectral function of the operator Ŝ+ q at the Rabi frequency through Without interactions the transition rate is exponentially small in the high frequency Ω because the only way the perturbation can excite a resonant transition is to give the particle extremely high momentum q ≈ √ m Ω, which is suppressed by the Gaussian resolution limit in the matrix element (6).In presence of interactions, on the other hand, a transition can satisfy energy conservation even for small momentum-transfer by utilizing the interaction term to create another particle-hole pair.The diagram of the lowest order process that can lead to a transition is shown in Fig. 2. Pictorially, the perturbation creates a particle hole pair with small momentum q (and concomitantly low energy).The excited particle with momentum k + q then scatters on another particle with momentum k inside the Fermi sea in a collision with large momentum transfer p k F (k F is the Fermi wave-vector) so that the pair of particles emerging from the collision are near the required final energy Ω.
The transition rate in this process can be framed as a Fermi golden rule calculation, equivalent to the imaginary part of the diagram in Fig. 2: Here the intermediate states are the particle hole states labeled by the hole momenta k and the accessible final states consist of two particles and two holes with S z = ±1, such as The effective interaction vertex merits a brief discussion.By a slight abuse of notation we write the contact using the same symbol used for the interaction Hamiltonian (1b), where V is the volume.However, it should be noted that the effective contact interaction is really the low energy limit of the T -matrix.In most cases it is given by the pseudopotential u = 4π 2 a/m with a the swave scattering length.But at unitarity the dependence on the energy of the scattered states is important as it cuts off the divergence at any non vanishing energy.The T -matrix appropriate for the vertex in Fig. 2 from the Lippman-Schwinger equation at unitarity is where n in the particle density and ρ( ) the single particle density of states per unit volume (DOS).We also define an interaction energy scale U = un.Substituting into the Fermi golden rule expression we obtain: The factors n k (1 − n k+q ) constrain the sum over k to a fraction ∼ |q|/k F of the Fermi surface volume.Physically, this is the phase space available for creating the low momentum particle-hole pair intermediate state.The other constraints are obeyed automatically due to the required large momentum transfer p in the collision.In addition, because we assume Ω F we shall neglect the small particle-hole energy ξ k − ξ k+q in the denominator and the hole energies −ξ k and −ξ k in the δ-function.Thus we obtain where N is the particle number.The total transition rate is the integral over q = qẑ up to the resolution cutoff The factor (k F L) −1 stems from the suppressed phase space for excitation of the fermi sea with low momentum transfer, while the logarithmic correction is due to the integration over the 1/q behavior of the function |v q | 2 q for q π/L.We remark that the scaling with Ω can be understood without a calculation by noting that the transition rate is a second order process with the intermediate states constrained to low energy compared to Ω due to the momentum cutoff on the perturbation.Thus the energy denominator due to the virtual transition is ∼ Ω and the transition rate must scale as Γ ∼ VgU Ω 2 ρ(Ω).Plugging in the DOS in 3d we get the correct scaling with Ω.The small momentum transfer q limits the phase space for intermediate states giving the |q|/k F factor in Γ q .
Finally we can convert the transition rate to a heating rate by multiplying it with the average energy increase per particle per unit time.In our experiment, the energy Ω/2 given to the excited particle pair is sufficient to escape the trap.Therefore, the heating is due to the holes left in the Fermi sea, leading to the heating rate ˙ = F Γ/N .To get a dimensionless measure of the heating rate we define η = ˙ / 2 F .η gives the relative energy change per particle δ / F , which occurs in the characteristic timescale of the system τ F = / F .If η 1, we can observe phenomena slow on the scale of F without being affected by the heating.At unitarity we get In our experiment, we can reach driving frequencies Ω/2π > 10kHz where this parameter is η < 10 −4 , so heating is not expected to be a problem.Before proceeding we note that a semi-classical calculation of the heating induced by periodic forces acting on an ultracold gas was reported in Ref. [46].

III. EXPERIMENT
Our experiments are performed with a quantum degenerate gas of 40 K atoms, prepared in an incoherent, spin-balanced mixture of the two lowest energy states, |↓ = |9/2, −9/2 and |↑ = |9/2, −7/2 , with the notation |F, m F .The flat trap, V trap , is created by three laser beams with a wavelength of 532nm [38,47] (see Fig. 1); a 'tube' beam is created by a wide Gaussian beam (125µm waist radius) that has a circular hole at its center, created by a digital mirror device [48].The other two 'end-cap' beams are created by two highly elliptical Gaussian beams with waist radii of 5.5µm and 180µm.Together, they generate a dark cylindrical volume with an approximate height of 39µm and a diameter of 55µm, defined by the full width at half maximum of the atomic density.The cylinder symmetry axis is parallel to the gravitational force.
The experimental sequence starts by cooling the gas to quantum degeneracy in a crossed optical dipole trap [49].To improve the loading efficiency into the flat trap, we have added a second crossing beam to the optical trap described in Ref. [49].This yields a harmonic trap with trapping frequencies of ω r = 2π × 236(1)Hz and ω z = 2π × 27(2)Hz, in the radial and axial directions, respectively.After forced evaporation, there are N ≈ 5 × 10 5 atoms at T /T F ≈ 0.24 in this trap, where N is total atom number in both spin states and T F the Fermi temperature.
To load the flat trap, the tube beam is ramped to 30mW already at the beginning of the evaporation in the harmonic trap.The magnetic field gradient that counteracts gravity is ramped to its final value, as given in Eq.( 3), in 0.5s, 1s before the harmonic trap is turned off.In our system, we have an additional undesirable small magnetic gradient of dBz dy ≈ 0.68G/cm in the transverse direction, which we compensate with another pair of coils.The sequence continues with a ramp up of the caps and tube beams power to 50mW and 150mW, respectively.The two traps are held overlapping for 50ms, and then the harmonic trap is ramped down in 200ms.We typically load around 32% of the atoms into the flat trap.Finally, the atoms are cooled in the flat trap by evaporation, forced by ramping down the power of the caps and tube beams in 2s to 20mW and 50mW, respectively.To ensure the cloud has reached equilibrium, we wait for an additional 0.8s before performing a measurement.The final typical conditions in the flat trap with the rf field are N ≈ 60 × 10 3 atoms with T /T F ≈ 0.15.The typical Fermi energy, F /h ≈ 940Hz, is determined independently from an in situ density measurement of the gas (see Appendix A).The magnetic field is tuned around the Feshbach resonance, at 202.14G, determining the strength of interactions, 1/k F a. It is ramped adiabatically to its final value in 10ms, where it is typically kept for 400ms.In the last 200ms of the experimental sequence, the rf pulse is turned on with a typical Rabi frequency of around Ω/2π ≈ 10.5kHz.

IV. RESULTS
Relaxation dynamics.-Priorto turning the rf field on, the atomic densities of the two spins are not uniform (Fig. 3a), because the magnetic field gradient given by Eq.( 3) over (under) compensates gravity for state |↓ (|↑ ) by 5.9%.Once the resonant rf field is turned on, the densities start to equilibrate.We study this relaxation process by preparing the gas with only spin |↓ atoms, and imaging them from the side of the cylinder after different waiting times (Fig. 3a-c).To quantify the non-uniformity of the gas, we plot the slope of the optical depth (OD) at the center of the trap, normalized by its maximal value at t = 0 (main part of Fig. 3).The gas re-laxes to a uniform density with damped oscillations.The oscillation frequency is roughly given by the time it takes an atom with a Fermi velocity to traverse the trap height back and forth.The density reaches a steady-state for rf pulse duration longer than 100ms.The residual density inhomogeneity due to the finite steepness and imperfections of the trapping potential is discussed in Appendix A.
Since spin-polarized fermions do not interact via swave scattering, a question that may arise is how the relaxation process actually occurs.While the rf rotation creates a |↑ component on a short time-scale of the inverse Rabi frequency, the gas remains spin polarized, albeit with the spin oriented in a different direction.However, small spatial inhomogeneities in the rf and magnetic fields lead to spin dephasing, and eventually, through atomic diffusion, also to decoherence.Therefore, the spin-polarized gas becomes a balanced spin mixture on a relatively short timescale of 10ms [50].Thanks to the invariance of the interaction Hamiltonian, the rf resonance frequency does not change as the gas transforms from being non-interacting to strongly-interacting [45].
Momentum distribution of the uniform gas.-Animportant issue to consider is heating, which may occur during the initial relaxation phase and during the continuous operation of the rf pulse.We obtain the temperature of the gas by measuring its momentum distribution.Ordinarily, this is done by letting the gas expand ballistically either in free space or in a harmonic trap [51].Due to the relatively large initial size of the cloud, the free expansion requires particularly long expansion times, which are not always feasible.Expansion in a harmonic trap, on the other hand, is done for a quarter of the trap period, but is sensitive to anharmonicity of the trap [40,41,51,52].
Here, we take a different approach and use Raman spectroscopy, which has the advantage that it can be applied to a trapped gas [53,54].The technique relies on a linear relation between the two-photon Raman detuning and the velocity of the atoms which are transferred from state |↑ to the initially unoccupied state |9/2, −5/2 .In the experiment, the two Raman beams are pulsed after the application of a 200ms-long rf pulse.By scanning the relative frequency between the beams, we obtain a spectrum that is directly proportional to the one-dimensional momentum distribution [53].A typical result with dynamically-driven uniform Fermi gas is shown in Fig. 4. The number of atoms in state |9/2, −5/2 is measured by selectively capturing them in a magneto-optical trap (MOT) and recording their fluorescence [49,53].To improve the detection, we separated the wavelength of the MOT, which is close to the D 2 transition, from that of the collected scattered photons [55].To this end, we added a dedicated probe beam, tuned to the D 1 transition, and filtered the recorded image with an ultra-narrow, 1nm, optical band-pass filter [56].The intensities of the two Raman beams are actively stabilized and programmed to follow a 1ms-long Blackman FIG. 4. The one-dimensional momentum distribution of a uniform periodically-driven Fermi gas.The distribution, measured by Raman spectroscopy [53], is fitted with three models: a numerical model that provides a realistic description of the flat trap in the experiment (solid blue line, see Appendix A), an ideal homogeneous Fermi-Dirac distribution given by Eq.( 18) (dotted black line), and a harmonically trapped gas (dashed red line).The numerical and homogeneous models yield a better fit to the data than the harmonic one, based on the χ 2 values.The differences between the uniform and numerical models are not significant.The temperature extracted from the numerical model is T /TF = 0.15 (1) with T = 7(2)nK.The results of the other fits are discussed in the main text.This measurement was performed at a magnetic field of B = 209.18G,where the atoms are very weakly interacting.The rf pulse duration is 200ms with a Rabi frequency of Ω/2π ≈ 10.5kHz.Error bars represent one standard deviation of the measured values.
pulse [57].The one-photon Raman detuning is around 46.1GHz below the D 1 transition.To reduce unwanted single-photon scattering, which constitutes most of the background signal, we incorporated a temperature stabilized etalon after the Raman laser to filter the broadband amplified spontaneous emission.
We analyze the momentum distribution by fitting it with three different models (see Fig. 4).The first one is a numerical local-density approximation model of a gas in a realistic flat trap that is used in our experiment (solid blue line).In this model, we account for the finite steepness of the trap walls, which is calibrated using in situ density images (see Appendix A).The free parameters are the reduced temperature (T /T F ) and the background offset.For comparison, we fit the data with two ideal models of harmonically trapped gas (dashed red line) [53] and ideal uniform gas (dotted black line).The doubly-integrated momentum distribution of the latter is given by: where k z is in the direction of the two-photon momentum transfer [53], and ζ is the fugacity with the implicit , with Li n (z) being the Polylogarithm function.Notice that due to the doubleintegration, there is no sharp Fermi surface in this functional even at T = 0.
The numerical and uniform models fit the data markedly better than the harmonic one, as evident by comparing their χ 2 fit values (see legend of Fig. 4).The temperature extracted from the numerical model is T /T F = 0.15(1) with T = 7(2)nK.The uniform model yields a very similar result of T /T F = 0.16 (2).In contrast, the harmonic model gives a much lower temperature of T /T F = 0.08 (2).In the numerical model, k F is calculated directly from in situ density images (see Appendix A), while in the uniform and harmonic models, it is left as a free fitting parameter.We find that the k F extracted from the uniform model is only 6(2)% higher than the one calculated directly, while that of the harmonically trapped model is higher by 30(2)%.
To test whether the rf-induced spin rotation causes heating, we repeat the Raman measurement at different rf pulse durations, shown in Fig. 5. Within the experimental accuracy, we do not observe an increase of the temperature.This result is corroborated by the measurements of the condensate fraction versus time, described below.We therefore conclude that with our experimental parameters, satisfying Ω F / , the heating rate is too small to be detected.
Pair condensation.-Wenow turn to probe the many-body properties of a dynamically driven Fermi gas.When a spin-balanced Fermi gas is cooled below the critical temperature, T c , atoms with opposite spins pair and condense, forming a fermionic superfluid [26,[58][59][60].The value of T c depends on the interaction strength.The survival of superfluidity is a stringent test of our Floquet engineering scheme, since this phase is extremely sensitive to heating and differential forces acting on the spins.
In these experiments, we cool the gas below the su-perfluid transition at unitarity (1/k F a ≈ 0).At the end of the cooling stage, the magnetic field is ramped in 10ms from 203.5G (weak interactions) to unitarity.There, it is held for 5ms, during which the atoms pair-up and condense.The magnetic field gradient and rf pulse are present during the last 200ms, long enough to ensure equilibrium.Since during this time the magnetic field is changing, we program the rf frequency to track the resonance transition.
We characterize the survival of the superfluid phase by measuring the condensate fraction, using the pairprojection technique [59,61].To this end, the trap is abruptly turned off and at the same time the magnetic field is ramped rapidly (40µs) to the BEC side of the resonance (199.8G)[59].This procedure projects the looselybound pairs onto tightly-bound molecules.We then let the gas expand for 24ms and measure the distribution using absorption imaging.For the imaging to work, we dissociate the molecules by ramping back the magnetic field to unitarity just before taking the image.When the gas is superfluid, the recorded density distribution is bimodal, with condensed pairs appearing as a pronounced central peak (see Appendix B).
In Fig. 6, we plot the condensate fraction (lower panel) and total number of atoms (upper panel) at unitarity, extracted from the images of the expanded gas, as a function of the driving frequency.The overlap between the two spin distributions is large enough even without the rf field to yield a condensate fraction of around 0.11.We distinguish between three frequency regimes with qualitatively different behavior.At very low frequencies, we observe a quasi-static behavior with local equilibrium and no apparent change in the conditions of the gas.As the frequency is increased, we cross to the second regime, where the driving generates spin currents and micromotion that clearly harms the superfluid.Initially, at around Ω ∼ 100Hz, the heating does not lead to a loss of atoms because the energy of the excitation is smaller than the trap depth.In fact, as the condensate fraction decreases, the number of detected atoms increases.This surprising behavior is a result of a partial correlation between the number of detected atoms following the pairprojection technique and the number of pairs.It exists since during the magnetic field ramping a small fraction of the pairs are lost, most likely to deeply bound molecular states.At higher frequencies, atoms have enough energy to leave the trap, but heating is gradually suppressed due to the scaling with Ω of the second-order process (see section II).In the third regime, defined by η < 1, the loss decreases and the condensate revives.In the high frequency limit, where η 1, the atom number returns to its initial value while the condensate fraction reaches an even higher value, an effect we attribute to a better spatial overlap between the spins in a uniform gas.
We now return to the question of heating at high driving frequency.As shown in Fig. 5, we do not observe a rise of the temperature, as extracted from the momentum distribution.However, in the superfluid phase, the con- Relative heating rate, η FIG. 6. Condensate fraction and total atom number at unitarity as a function of the driving frequency.At very low frequencies, the spin oscillation is slow enough such that both observables are the same as without the driving, which are marked by horizontal lines with shading representing the uncertainty.At higher frequencies, the driving has an adverse effect on the superfluid.The vertical dashed line marks the frequency at which η = 1 (see Eq.( 17)).Above it, heating starts to be suppressed (the color represents the value of η).At even higher frequencies, the atom number returns to its initial value, and the condensate fraction reaches an even higher value than for a stationary gas (inset).The conditions are measured after 5ms at the Feshbach resonance magnetic field, and the Rabi frequency is varied by changing the rf pulse power.The observables extraction procedure is discussed in Appendix B. The error bars combine the 1σ mean confidence interval of the fit with statistical errors over 10 repetitions.
densate fraction is a much more sensitive thermometer.In Fig. 7, we plot the total number of atoms and condensate fraction versus the waiting time at unitarity (black circles).We employ an rf field with a relatively high Rabi frequency (Ω/2π ≈ 10.5kHz), where the density is already uniform and the condensate fraction reaches its high value (see Fig. 6).To distinguish between loss and heating due to the rf driving and other sources, we repeat this measurement without the rf field (red squares).The data, taken with a waiting time of up to more than 1s, do not point to any heating, as there is no reduction of the condensate fraction, even when the total number of atoms has decreased by approximately a factor of two.The decay in the number of atoms is almost identical with or without the rf pulse.We analyze the data using the following loss model [62][63][64] where n is the total atomic density.K 1 = 1/13.5s −1 is the single-body loss rate, determined by the rate of collisions with the residual gas in the vacuum chamber, and .Time dependence of the atom number and condensate fraction at unitarity.Data is taken after different waiting durations both with (black circles) and without (red squares) the rf driving.In both cases, the loss has a similar trend, which is well-fitted by the model in Eq.( 19) (solid lines).Inset: The condensate fraction (same marks) is plotted together with the weighted average (solid lines), and its standard deviation (shades), shows no decrease.Error bars are determined as in Fig. 6.
measured independently.K 2 and K 3 are the two-and three-body loss rate coefficients.Previously, these parameters were measured only with harmonically trapped gases, which complicated the analysis due to the nonlinear density dependence in this model.Here, we benefit directly from the uniformity of the gas and from the fact its shape and volume are almost unchanged as the atom number diminishes.Fitting the data taken with the rf pulse with both coefficients as free parameters (black solid line) yields K 3 = 9(1) × 10 −25 cm 6 s −1 and K 2 = 0.This shows that the loss is mainly due to three-body recombination.Since in the data without the rf pulse the density is not homogeneous, and in fact, differs between the two spin components, we do not use it to extract loss coefficients.Qualitatively, however, it is still fitted well by the model of Eq.( 19) (red solid line).Our value for K 3 is 10 times higher than the one measured in a harmonic trap and at a significantly higher temperature in Ref. [63].We note, however, that their maximal value of K 3 , which was measured on the BEC side of the resonance, agrees with our measurement at unitarity.The contact parameter of a uniform gas.-We now turn to a measurement of the contact parameter in the BEC-BCS crossover regime.This parameter is central to a set of universal thermodynamic and energetic relations [65][66][67][68][69][70], many of which have been tested experimentally [71][72][73][74].Previous works determined the value of the contact with harmonically trapped gas at different temperatures and interaction strengths [71,75,76].Local measurements resolved the contact of a quasihomogeneous sample [33,35,36,[77][78][79][80].Until now, the contact of a truly uniform gas was measured only at unitarity [81].
We determine the contact from the power-law tail of rf line-shapes taken with the uniform gas [33,49,71,81].In contrast to the condensate fraction experiments, where the condensate was formed after ∼ 2ms, here we observed that it takes at least 100ms for the tail to fully develop.For this reason, we wait for 400ms in the final magnetic field before measuring the rf line-shape.The number of atoms and temperature are similar to the condensate fraction experiments.The spin-rotation rf field, with Ω/2π ≈ 10.5kHz, is turned on for the last 200ms.It is turned off 0.5ms before we probe with a 1ms square pulse of a second rf field, whose frequency we scan near the |↑ → |9/2, −5/2 transition.The atom number in state |9/2, −5/2 is again detected with fluorescence imaging [49].A typical rf lineshape is shown on a logarithmic scale in the inset of Fig. 8.A universal power-law tail over two decades is clearly visible.To extract the contact, we work in natural Fermi units and normalize the spectrum to 1/2.The tail is then fitted with C/(2 2/3 π 2 ν 3/2 ) (black line in the inset), where C is the contact parameter in units of N k F .Owing to the high sensitivity of our fluorescence detection scheme, we keep the rf power constant for all detunings, while the maximal transferred fraction is no more than 8%.The systematic error in the determination of the contact due to the remaining density inhomogeneity in our realization of the flat trap is estimated by calculating the densityweighted average of a theoretical contact [82] using our calibrated model of the trap (see Appendix A).We find that the density-averaged contact, representing our measured contact, is lower than the homogeneous contact by 5% in the BCS side (1/k F a = −1), and higher by 7% in the BEC side (1/k F a = 1).Near unitarity, the large scattering length reduces the deviation to less than 1%.As a comparison, for harmonically trapped gas at the same average density, the systematic errors are 19% and 25% in the BCS and BEC sides, respectively.
In Fig. 8 we plot the contact of a uniform Fermi gas at various interaction strengths in the BEC-BCS crossover.Starting from the BCS side (a < 0), the contact increases monotonically towards the BEC side of the crossover, where it converges to the asymptotic behavior of a molecule, C BEC = 4π/k F a [83].We find that already above 1/k F a ≈ 0.8, our data are very close to C BEC .In contrast, the weak-coupling BCS limit of the contact, C BCS = 4 (k F a) 2 /3 [83], is not attained even at 1/k F a = −1.We compare our data to several theories and previous measurements.On the BEC side, there is a pronounced difference between the T = 0 and T = T c predictions [82].Our data, which were taken slightly below T c , agrees with the T = 0 T-matrix calculation.We also find a good agreement with the Gaussian pair fluctuations (GPF) calculation [84,85] and fixed-node diffusion Monte Carlo simulation (FNDMC) [79,86], especially in the BEC region.A Luttinger-Ward calculation [83] is slightly below our data on the BEC side.FIG. 8.The contact of a uniform Fermi gas in the BEC-BCS crossover.The contact is extracted from the tail of rf line-shapes taken at different interaction parameters.As an example, the inset shows the line-shape at (kF a) −1 = 0.75(2) together with its fit (black solid line).Upper panel: the theoretical prediction for the contact in the BCS (BEC) limit is shown as a loosely (densely) dash-doubledot line.We compare our data (blue circles) with a non-selfconsistent T-matrix model at T = Tc improved by Popov (dash-dotted line) and at T = 0 (solid line) [82], a Luttinger-Ward calculation (dashed line) [83] and a GPF calculation (dotted line) [84,85].The FNDMC line [79,86] is indistinguishable from the GPF line on this scale.Lower panel: for a better quantitative comparison, we defined a shifted contact as C = C − 3 exp (1.4/kF a).In this plot we also add previous measurements done with a quasi-homogeneous gas (green diamonds) [35].Error bars represent 1σ confidence interval of the fit.

V. DISCUSSION
In this work, we have demonstrated that Floquet engineering can be used with a continuous interacting Fermi gas without affecting its intrinsic many-body behavior.Specifically, we have employed the technique to eliminate the effect of a spin-dependent potential and achieve a flat trap.Our experiments are done with a driving frequency that is much higher than all other relevant experimental scales.In this regime, we have found no detectable heating during the experiment due to periodic driving.Measurements of the condensate fraction and contact parameter show the same behavior as expected in a stationary uniform Fermi gas.Furthermore, our dynamical levitation scheme can be used to generate a uniform density of other spin mixtures.
The full Hamiltonian of Eq.( 1) depends explicitly on time and therefore does not conserve energy.In contrast to many-body localized systems [87][88][89][90], our gas is ergodic and thus it is not protected from heating [16][17][18].Nevertheless, interacting many-body systems can attain long lived prethermal states, following a quench of the Hamiltonian parameters [91][92][93][94][95] or upon initiation of periodic driving [20,21,96,97].In particular discrete lattice systems show an exponentially slow energy absorption rate at high frequencies leading to Floquet prethermal states that persist for a time exponentially long with the drive frequency Ω [19][20][21][22][23]25].The Fermi gas in the continuum does not benefit from such exponential suppression of the heating rate.Nonetheless, we have shown in section II that the smooth spatial structure of the periodic perturbation leads to parametric suppression of the heating rate in a unitary Fermi gas.The suppression is controlled by the small values of V g /( Ω), with V g being the strength of the time-dependent perturbation, F /( Ω), and 1/k F L. This allows to obtain prethermal states with lifetimes much longer than the characteristic many-body timescale.In this work, we have focused on demonstrating this fact with measurements of the momentum distribution, the contact parameter, and the condensate fraction.It will be interesting to investigate the heating process and thermodynamic properties of the driven gas in the future.
Our work should also be placed in the framework of quantum information processing, where ultracold atoms are proposed as a resource for quantum memory [98][99][100].A common cause of decoherence is spatially inhomogeneous spin-dependent potentials.As an example, the energy difference between two internal states of optically trapped atoms usually varies in position due to differential light shifts.In a classical ensemble, each atom can be treated independently with the rest of the ensemble acting as a fluctuating bath [101,102].These fluctuations lead to decoherence of the qubit stored in this atom.Dynamical decoupling [103,104], a generalization of the celebrated Hahn echo technique [105], can substantially slow this relaxation process by applying multiple spin rotations [106].
Dynamical decoupling has been applied successfully in NMR [107][108][109], photonic systems [110], trapped ions [111,112], electron spin in solids [113][114][115][116], ultracold atoms [106,117], and Bose-Einstein condensates (BEC) [118,119].In all cases, the decoupled system was weakly interacting and could be treated in a mean-field approach.Dynamical decoupling was not applied before to a strongly interacting ensemble with the aim of preserving its many-body behavior.In our work, the spindependent potential originates from the magnetic cancellation of gravity.The spin-rotation rf pulse we apply is a continuous version of a simple dynamical decoupling scheme.The absence of heating we observe is promising for future explorations of more sophisticated sequences tailored to generate specific local and global symmetries [24].For example, the realization of tilted Fermi-Hubbard chains where the rf dressing was used to tune a relative tilt difference of two different spin states [120].
|9/2, −5/2 , which is detuned by 92MHz from the optical transition.The last pulse ensures that no artifacts are introduced to the imaging due to large atom number in off-resonant states.The two pulses are completed within less than 70µs, ensuring the density is unchanged during this procedure.
The spin-polarized gas is essentially non-interacting and can be described by a Fermi-Dirac distribution.To fit the two-dimensional integrated density image, we calculate the density in a local density approximation [121]: A1) where λ dB is the de Broglie wavelength, β = 1/(k B T ) with k B being the Boltzmann constant, U (r, z) is the trapping potential and µ is the chemical potential, which is set by the total number of atoms.The model potential of the tube beam is parametrized by a power law function, while the potential of the cap beams is taken as a Gaussian function: Here r is the radial coordinate relative to the symmetry axis of the tube (denoted by z), and U r is the potential barrier of the cylindrical wall.The tube radius σ r = 32µm and the power-law exponent p = 13.6 are extracted from a direct measurement of the laser beam that generates the potential [55].U z is the potential barrier of the cap beams, σ z = 7.5µm is their waist radii in the z direction, also measured directly.z 0 = 24µm is half the separation between the two cap beams, measured by imprinting the caps profile on a dilute expanded cloud of atoms [55].The temperature, T , is found selfconsistently together with momentum distribution measurements.The only free parameters in the fit are U r , U z and the center position of the fit.A two-dimensional density function for fitting the in situ image is generated by integrating the three-dimensional density n (r, z) along one of the Cartesian axes.An example of a typical calibration is shown in Fig. 9. Once the model potential is calibrated, we calculate F from the peak density, E F = 6π 2 n F (0, 0) 2/3 2 /2m.The density-averaged F is smaller by 12%.This difference is considerably smaller than in the harmonically trapped gas before loading the flat trap, in which the density-averaged F is smaller by 44% than the peak F .We quantify the density homogeneity in the form of a probability distribution function P (n).We calculate the axially-integrated probability distribution P (n 2D /n 2D ) as an histogram of absorption images taken along the symmetry axis (inset (b) of Fig. 9), where n 2D is the two-dimensional atomic density and n2D is its weighted average.The probability distribution is plotted in inset (e) of Fig. 9, together with the probability distribution of a harmonically trapped gas.A direct measure of uniformity is the standard deviation of P (n 2D /n 2D ), which  is 0.23 and 0.42 in units of n2D for the flat and harmonic traps, respectively.
The knowledge of the trapping potential enables us to calculate the one-dimensional momentum distribution n (k z ) and use it as a fitting function for the data acquired in the Raman spectroscopy experiment (Fig. 4).n (k z ) is obtained by integrating the semi-classical distribution [121]  approach to filter out the background noise in the images [122].We have verified that this noise removal procedure does not change significantly the reported CF values and only reduces the uncertainty.The recorded density can be roughly considered as dissociated pairs that have either non-zero or zero center-of-mass momentum.The latter are the condensed pairs which constitute the central peak of the image (see Fig. 10).Each of the two populations expands differently, according to their respective momentum distribution.The total atom number is extracted by a direct integration of the OD image, and the error bars indicate statistical standard error only.
Following the time of flight, the non-condensed part of the gas is characterized by a wider expansion with respect to the trap dimensions, while the condensed part spreads only slightly beyond its original size, set by the tube beam diameter [40].To separate between the condensed and non-condensed parts, we mask out the central region of the image (dashed line in Fig. 10) and fit only the tail of the azimuthally-averaged signal with the momentum distribution of a thermal gas of non-interacting bosons [121] (pink shading in Fig. 10), where the thermal de Broglie wavelength λ T = h/ √ 2πmk B T and the fugacity z, are fitted under the normalization constraint N th = 2π n th (k)kdk, and N bg accounts for the background signal.We extract the condensed population from the signal that lies above the fit (yellow shading in Fig. 10).The mask radius, R mask , should be large enough to leave only the thermal wings for fitting.When we analyze data taken with no condensate, we observe that the width of the fitted distribution is almost independent of the mask radius up to around R mask ≈ 80µm.For larger radii, the signal in the remaining thermal wings is too weak, and the fit exhibits a systematic deviation.Therefore, we set the mask radius to R mask = 75µm.Finally, we note that the CF values we obtain are close to those reported in Ref. [81], taking into account our measured reduced temperature.Lower panels present the corresponding azimuthally-averaged signals (blue lines).We fit Eq.(B1) (red line) to the thermal wings at radii > R mask (dashed black line).The condensate fraction is defined by the integrated signal above the fit extrapolation (yellow shading) over the total integrated signal (pink plus yellow shadings).To make a fair comparison, each of the two distributions is normalized by its total number of atoms, 61.0×10 3 and 68.3×10 3 in the right and left examples, respectively.

V
ext,# (a) < l a t e x i t s h a 1 _ b a s e 6 4 = " b 6 6 O A P A 0 b Y H / Z b 6 K U k 2 9 T z 7 7 F 0 k = " > A A A C L H i c b V D L S g M x F M 3 4 t j 6 r 7 t w E B 8 F V m Y i g S 8 G N S w V r h c 4 g d 9 L b G p p k h i S j 1 G m / w K 1 + h l / j S s S t 3 2 E 6 d u H r Q O B w z r n k 3 p P m U l g X R a / B 1 P T M 7 N z 8 w m J t a X l l d W 2 9 v n F p s 8 J w b P J M Z u Y q B Y t S a G w 6 4 S R e 5 Q Z B p R J b a f 9 k 7 L d u 0V i R 6 Q s 3 y D F R 0 N O i K z g 4 L 5 3 f X 6 + H U S O q Q P 8 S N i E h m e D s u h 5 s x Z 2 M F w q 1 4 x K s b b M o d 0 k J x g k u c V S L C 4 s 5 8 D 7 0 s O 2 p B o U 2 K a t N R 3 T X K x 3 a z Y x / 2 t F K / T 5 R g r J 2 o F K f V O B u 7 G 9 v L P 7 n t Q v X P U p K o f P C o e Z f H 3 U L S V 1 G x 2 f T j j D I n R x 4 A t w I v y v l N 2 C A O 1 9 O L d Z 4 x z O l Q H f K u I 9 u 1 G Z J O Q x Z b E D 3 q q O + B 1 I D V S C W l U t D N h z 5 H t n v 1 v6 S 1 n 6 D H T Q Y O z 8 I j 4 8 m l S 6 Q b b J D 9 g g j h + S Y n J I z 0 i S c I H k g j + Q p e A 5 e g r f g / S s 6 F U x m N s k P B B + f m B u n G w = = < / l a t e x i t > z < l a t e x i t s h a 1 _ b a s e 6 4 = " N j d I 7 s m x u o Q t c s d w v c d x K 1 T Dv Y g = " > A A A C L H i c b V D L S g M x F M 3 4 t r 4 f O z f B Q X B V J l L Q Z c G N S w V r h c 4 g d 9 L b G p p k h i S j l G m / w K 1 + h l / j S s S t 3 2 E 6 d u H r Q O B w z r n k 3 p P m U l g X R a / B z O z c / M L i 0 n J t Z X V t f W N z a / v K Z o X h 2 O K Z z M x1 C h a l 0 N h y w k m 8 z g 2 C S i W 2 0 8 H p x G / f o b E i 0 5 d u m G O i o K 9 F T 3 B w X r r o 3 2 y G U T 2 q Q P 8 S N i U h m e L 8 Z i v Y j b s Z L x R q x y V Y 2 2 F R 7 p I S j B N c 4 r g W F x Z z 4 A P o Y 8 d T D Q p t U l a b j u m B V 7 q 0 l x n / t K O V + n 2 i B G X t U K U + q c D d 2 t / e R P z P 6 x S u d 5 K U s 4 P w u O j S a U 1 s k W 2 y S 5 h 5 J A c k 1 N y R l q E E y Q P 5 J E 8 B c / B S / A W v H 9 F p 4 L J z A b 5 g e D j E w t j p s w = < / l a t e x i t > + < l a t e x i t s h a 1 _ b a s e 6 4 = " m C Z F 6 b b s 4 P w u O j S a U 1 s k W 2 y S 5 h 5 J A c k 1 N y R l q E E y Q P 5 J E 8 B c / B S / A W v H 9 F p 4 L J z A b 5 g e D j E w t j p s w = < / l a t e x i t > + box < l a t e x i t s h a 1 _ b a s e 6 4 = " 4 f 3 1

FIG. 1 .
FIG.1.Creating a uniform Fermi gas by periodic driving.(a) The gas, composed of two spin states (marked by red and blue colors and by opposite arrows), is trapped in a box-like optical potential.The two spins have different magnetic dipole moments.As a result, it is only possible to partially counteract the gravitational potential with a magnetic field gradient, Bz, set according to Eq.(3).The total external potential Vext,s depends on the spin, s ∈ {↑, ↓}, and consequently, the density distribution of each spin is different and not uniform (color gradient in the top left figure).(b) By adding a resonant rf field that drives rapid spin rotations, we create an effective spin-independent potential, in which the gas becomes homogeneous.Importantly, the intrinsic manybody behavior of the gas is unchanged by this driving.

FIG. 3 .
FIG.3.Relaxation to a uniform density.The gas is prepared in a spin-polarized |↓ state without the rf field.The magnetic field is set according to Eq.(3), overcompensating gravity for this state by 5.9%.At t = 0, the rf field with Ω/2π ≈ 15.7kHz is switched on and flattens the time-average potential.The in situ density of gas is recorded from the side of the cylindrical trap after variable rf pulse duration (insets showing 3ms, 17.5ms and 110ms).The colorbar represents the optical depth (OD) in the absorption image.The main figure shows the density relaxation dynamics, quantified by the average change in the OD from side to side of the image.The errorbars represent 1σ deviation from linear fit to the OD slope.We fit the data with g (t) = ae −t/τ cos (2πf t + φ) + b (solid line), and obtain τ = 16(1)ms and f = 33.6(7)Hz.The reduction of the OD for longer pulse duration is due to decoherence to a spin-balanced mixture while measuring atoms only in |↓ state.(see main text).This measurement is done at a magnetic field of ∼ 203G, where the scattering length between states |↓ and |↑ is a ≈ −1152a0 (a0 is Bohr radius).

FIG. 5 .
FIG.5.Reduced temperature as a function of the rf pulse duration.The temperature is extracted from onedimensional momentum distribution measurements fitted by a homogeneous Fermi-Dirac distribution.We observe no increase in T /TF up to 500ms, and this holds true also for T alone.Note that the data in this figure were taken at slightly different conditions than the one in Fig.4, with N ≈ 131×10 3 atoms at T /TF ≈ 0.23.The measurement was performed at a magnetic field of B = 209.18G,where the atoms are very weakly interacting.The Rabi frequency is Ω/2π ≈ 10.5kHz.Error bars represent one standard error of the fit.
FIG.7.Time dependence of the atom number and condensate fraction at unitarity.Data is taken after different waiting durations both with (black circles) and without (red squares) the rf driving.In both cases, the loss has a similar trend, which is well-fitted by the model in Eq.(19) (solid lines).Inset: The condensate fraction (same marks) is plotted together with the weighted average (solid lines), and its standard deviation (shades), shows no decrease.Error bars are determined as in Fig.6.

FIG. 9 .
FIG. 9.In situ density analysis.An example of an in situ atomic density measured by absorption imaging from the side (a) and from the top (b) (the cylinder symmetry axis is vertical).A vertical cuts along the side and top images (circles) together with a fit to the numerical model (solid lines) are shown in insets (c) and (d), respectively.The flat part of the cut is where the density is uniform.Inset (e) shows a 2D density probability distribution of the flat trap (blue line) and a harmonically trapped gas (red line).The density probability function of the flat trap is peaked around the average density value, indicating that most of the density is close to the average.Images (a) and (b) are an average of 5 and 10 experimental repetitions, respectively.

FIG. 10 .
FIG.10.Extraction of the condensate fraction.Upper panels are absorption images of high (left) and low (right) condensate fraction averaged over 10 experimental repetitions.Lower panels present the corresponding azimuthally-averaged signals (blue lines).We fit Eq.(B1) (red line) to the thermal wings at radii > R mask (dashed black line).The condensate fraction is defined by the integrated signal above the fit extrapolation (yellow shading) over the total integrated signal (pink plus yellow shadings).To make a fair comparison, each of the two distributions is normalized by its total number of atoms, 61.0×10 3 and 68.3×10 3 in the right and left examples, respectively.