Optical cooling and trapping of highly magnetic atoms: The benefits of a spontaneous spin polarization

From the study of long-range-interacting systems to the simulation of gauge fields, open-shell Lanthanide atoms with their large magnetic moment and narrow optical transitions open novel directions in the field of ultracold quantum gases. As for other atomic species, the magneto-optical trap (MOT) is the working horse of experiments but its operation is challenging, due to the large electronic spin of the atoms. Here we present an experimental study of narrow-line Dysprosium MOTs. We show that the combination of radiation pressure and gravitational forces leads to a spontaneous polarization of the electronic spin. The spin composition is measured using a Stern-Gerlach separation of spin levels, revealing that the gas becomes almost fully spin-polarized for large laser frequency detunings. In this regime, we reach the optimal operation of the MOT, with samples of typically $3\times 10^8$ atoms at a temperature of 15\,$\mu$K. The spin polarization reduces the complexity of the radiative cooling description, which allows for a simple model accounting for our measurements. We also measure the rate of density-dependent atom losses, finding good agreement with a model based on light-induced Van der Waals forces. A minimal two-body loss rate $\beta\sim 2\times10^{-11}\,$cm$^{3}$/s is reached in the spin-polarized regime. Our results constitute a benchmark for the experimental study of ultracold gases of magnetic Lanthanide atoms.


Introduction
Open-shell Lanthanide atoms bring up new perspectives in the field of ultracold quantum gases, based on their unique physical properties. Their giant magnetic moment allows exploring the behavior of long-range interacting dipolar systems beyond previously accessible regimes [1,2,3,4,5,6,7]. The large electronic spin also leads to complex low-energy scattering between atoms, exhibiting chaotic behavior [8,9]. The atomic spectrum, which includes narrow optical transitions, further permits the efficient production of artificial gauge fields [10,11,12].
For two-electron atoms like Ca, Sr and Yb [22,23,24], the absence of electronic spin simplifies the operation and theoretical understanding of the MOT. On the contrary, the large value of the electron spin for Lanthanide atoms greatly complicates the atom dynamics, the modeling of which a priori requires accounting for optical pumping effects between numerous spin levels. Previous works on narrow-line MOTs with Lanthanide atoms mentioned a spontaneous spin polarization of the atomic sample [19,1,2], but its quantitative analysis was not explicitely described.
In this article, we present a study of Dy magneto-optical traps operated on the 626-nm optical transition (linewidth Γ = 2π × 136 kHz [25]) [15]. We measure the spin populations using a Stern-Gerlach separation of the spin levels. We observe that, for large and negative laser detunings (laser frequency on the red of the optical transition), the atomic sample becomes spin-polarized in the absolute ground state |J = 8, m J = −8 . This spontaneous polarization occurs due to the effect of gravity, which pushes the atoms to a region with a relatively large magnetic field (on the order of 1 G), leading to efficient optical pumping [19,1,2]. In the spin-polarized regime, the system can be simply described with a two-level atom model, in close relation with previous works on Sr magneto-optical traps [26] and narrow-line MOTs of magnetic Lanthanides [19,1,2].
We show that the spin-polarized regime corresponds to optimal operating parameters for the magneto-optical trap, leading to samples with up to N 3 × 10 8 atoms and temperatures down to T 15 µK. This observation is supported by a study of density-dependent atom losses triggered by light-induced Van der Waals interactions between atoms. We observe that minimal loss rates are reached in the spin-polarized regime, as predicted by a simple model of atom dynamics in attractive molecular states. ( 1 P 1 )(8,1) 9 (J = 9) and 4f 10 ( 5 I 8 )6s6p( 3 P 1 )(8,1) 9 (J = 9). (b) Scheme of the magneto-optical trap arrangement, with the six MOT beams (red arrows) and the pair of coils in anti-Helmoltz configuration. Gravity is oriented along −z. (c) Typical in situ absorption image of an atomic sample held in the 'compressed' magneto-optical trap. The cross indicates the location of the quadrupole center. The cloud is shifted by ∼ 1 cm below the zero due to gravity. We attribute the cloud asymmetry with respect to the z axis to the effect of an ambient magnetic field gradient. (d) Atom number N captured in the magneto-optical trap as a function of the loading time t (see text for the MOT parameters).

Preparation of magneto-optically trapped Dysprosium gases
The electronic states of Dysprosium involved in our study are represented in figure 1a, together with a schematics of the magneto-optical trap in figure 1b. An atomic beam is emitted by an effusion cell oven, heated up to 1100°C. The 164 Dy atoms are decelerated in a Zeeman slower, which is built in a spin-flip configuration and operates on the broad optical transition at 421 nm (linewidth 2π × 32 MHz [27]). The flux of atoms along the Zeeman slower axis is enhanced using transverse Doppler cooling at the entrance of the Zeeman slower, also performed on the 421-nm resonance [28].
The magneto-optical trap, loaded in the center of a steel chamber, uses a quadrupole magnetic field of gradient G = 1.71 G/cm along the strong horizontal axis x. The MOT beams, oriented as pictured in figure 1b, operate at a frequency on the red of the optical transition at 626 nm [15], the detuning from resonance being denoted ∆ hereafter. This transition connects the electronic ground state 4f 10 6s 2 ( 5 I 8 ), of angular momentum J = 8 and Landé factor g J 1.24, to the electronic state 4f 10 ( 5 I 8 ) 6s 6p( 3 P 1 ) (8,1) 9 , of angular momentum J = 9 and Landé factor g J 1.29. Its linewidth Γ = 2π × 136 kHz, corresponding to a saturation intensity I sat = 72 µW/cm 2 , allows in principle for Doppler cooling down to the temperature T D = Γ/(2k B ) 3.3 µK. Each MOT beam is prepared with a waist w 20 mm and an intensity I = 3.7 mW/cm 2 on the beam axis, corresponding to a saturation parameter s ≡ I/I sat 50.
The atom loading rate is increased by artificially broadening the MOT beam frequency: The laser frequency is sinusoidally modulated at 135 kHz, over a total frequency range of 6 MHz, and with a mean laser detuning ∆ = −2π × 4.2 MHz. From a typical atom loading curve (see figure 1d) one obtains a loading rate of 6(1)×10 7 atoms/s at short times and a maximum atom number N = 3.1(5) × 10 8 . The atom number is determined up to a 20% systematic error using absorption imaging with resonant light on the broad optical transition, taking into account the variation of scattering crosssections among the spin manifold expected for our imaging setup.
After a loading duration of ∼ 6 s, we switch off the magnetic field of the Zeeman slower, as well as the slowing laser. We then compress the magneto-optical trap by ramping down the frequency modulation, followed by decreasing the saturation parameter s and the laser detuning ∆ over a total duration of 430 ms. At the same time, the magnetic field gradient is ramped to its final value. For most of the MOT configurations used for this study, we do not observe significant atomic loss during the compression. A typical absorption image of the atomic sample after compression is shown in figure 1c. The curved shape of the gas and its mean position (about 1 cm below the magnetic field zero) reveal the role of gravity in the magneto-optical trapping [23,26,19].

Spin composition
An immediate striking difference of narrow-line MOTs with respect to alkali-metal ones is the strong dependence of the MOT center position on detuning [26]. As shown in figure 2a, we indeed observe a drop of the MOT position when increasing the laser detuning, the amplitude of which largely exceeds the cloud size (see figure 1c). This behavior can be explained by considering the MOT equilibrium condition, which requires mean radiative forces to compensate for gravity. When the laser detuning is increased, the MOT position adapts so as to keep the mean amplitude of radiative forces constant. This picture is supported by the calculation of the local detunings ∆ (m J →m J ) loc of optical transitions |J = 8, m J → |J = 9, m J at the MOT position. We compare in figure 2c these detuning values for the σ − , π and σ + transitions starting from the ground state |J = 8, m J = −8 . We observe that, when increasing the detuning ∆, the π and σ + transitions become off-resonant, while the local detuning of the σ − transition tends to a finite value. The detuning of the latter transition is denoted in the following as This predominance of the σ − transition leads to optical pumping of the electronic spin towards the absolute ground state, as confirmed by a direct measurement of the spin composition of the atomic sample with a Stern-Gerlach separation of spin levels. To achieve this, we release the atoms from the MOT, apply a vertical magnetic field gradient of about 30 G/cm during ∼ 4 ms, and let the atoms expand for 20 ms before taking an absorption picture (see figures 3a,b,c). We measure the spin composition for several detuning values ∆, and we show two examples of spin population distributions in the insets of figure 3d. The mean spin projection J z inferred from these data is plotted as a function of the detuning ∆ in figure 3d. It reveals an almost full polarization for As the MOT position deviates from the magnetic field zero position, optical transitions of polarization σ − become more resonant than π and σ + transitions, leading to optical pumping into the absolute ground state. (c) Local detunings for the three optical transitions involving the absolute ground state, inferred from the laser detuning and the Zeeman shifts at the MOT position. We take into account an ambient magnetic field gradient along z, δG = −0.094(2) G/cm, measured independently. The local detuning of the σ − transition (circles) saturates to a fixed value for large detunings, showing that the MOT position adapts to keep the local detuning fixed. The π and σ + transitions (square and diamonds, respectively) do not exhibit this saturation.
large detunings ∆ −2π × 1 MHz, which is the first main result of our study.
In the spin-polarized regime, the theoretical description of the magneto-optical trap can be simplified. As the atomic gas is spin-polarized and the σ − component of the MOT light dominates over other polarizations, the atom electronic states can be restricted to a two-level system, with a ground state |J = 8, m J = −8 and an excited state |J = 9, m J = −9 [26]. The radiative force is then calculated by summing the contributions of the six MOT beams, projected on the σ − polarization. For the sake of simplicity we restrict here the discussion to the motion on the z axis, but extending the model to describe motion along x and y directions is straightforward (see Appendix B). The motion along z is governed by the radiative forces induced by the four MOT beams propagating in the plane x = 0 (see figure 1b). Taking only the σ − component of these beams into account and neglecting interference effect between the beams, the total force for an atom at rest on the z axis is obtained by summing the radiation pressure forces from the four beams in the x = 0 plane. It takes the simple form F rad = kΓ 2 s 1 + 2s + 4∆ 2 loc /Γ 2 e z at the MOT position [29] ‡. The MOT equilibrium position corresponds to the condition of the radiative force compensating gravity F rad + mg = 0, leading to where η 168 for the considered optical transition. The local detuning ∆ loc thus only depends on the laser intensity s and not on the bare detuning ∆, as (2) ‡ For a different MOT beam geometrical configuration, with propagation directions along ±e x , ±e y and ±e z , we would obtain the same expressions for the radiative force, as well as the optical pumping rates described in the following.
This expression accounts well for the experimental data presented in figure 2c: the saturation of the local detuning ∆ loc for large detunings corresponds to the spinpolarized regime, where equation (2) applies. For simplicity we do not take into account magnetic forces associated with the magnetic field gradient, as it leads to ∼10% corrections on the local detuning ∆ loc , which is below our experimental error bars.
To go further this simple approach, we also developed a model taking into account the populations in all Zeeman sublevels. The Zeeman populations in the ground and excited states are calculated as the stationary state of optical pumping rate equations, including the Zeeman shifts corresponding to a given position. We then calculate the radiative force by summing the contributions of all optical transitions, which allows determining the MOT position z c from the requirement of compensation of the radiative force and gravity. As shown in figure 3d, this model accounts well for the measured population distributions.

Equilibrium temperature
The main interest in using narrow optical transitions for magneto-optical trapping lies in the low equilibrium temperatures, typically in the 20 µK range. In order to investigate the effect of the spin composition of the gas on its equilibrium temperature, we investigated the variation of the temperature T -measured after time-of-flight § -with the laser detuning ∆ (see figure 4). Far from resonance, we observe that the temperature does not depend on ∆, in agreement with the picture of the spinpolarized regime discussed above (see figure 4a). The temperature slightly decreases for −1 MHz < ∆/(2π) < −0.3 MHz, i.e. when leaving the spin-polarized regime, before significantly raising closer to resonance.
We also investigated the influence of the laser intensity, by probing the temperature variation with the saturation parameter s (see figure 4b). The observed temperature raise upon increasing s can be intuitively understood from equation (2): the local detuning from resonance increases when raising s; we then expect a temperature increase according to the Doppler cooling theory (in the regime |∆ loc | > Γ/2 considered here).
A more quantitative understanding requires adapting the theory of Doppler cooling to the experiment geometry. Here we restrict the discussion to the atom dynamics along z, in the spin-polarized regime (see Appendix B for a generalization to 3D). The radiative force produced by the four MOT beams in the x = 0 plane can be calculated for an atom of velocity v along z, leading to the damping force for small velocities which coincides with the usual damping force formula for Doppler cooling for a pair of counter-propagating laser beams of saturation parameter s and detuning ∆ loc , up to § The magnetic field gradient is kept on during expansion in order to avoid eddy current effects. We checked that magnetic forces play a negligible role in the expansion dynamics for the flight durations used for this measurement. numerical factors related to the geometry of our laser configuration. In Appendix A, we discuss additional experiments on temperature equilibration dynamics, which can be explained qualitatively using the damping rate value (3). The momentum diffusion coefficient along z, denoted as D zz , is calculated taking into account the contribution of the six cooling laser beams [30], as where Π = 1/η is the population in the excited state. The temperature T is then obtained as the ratio k B T = D zz /(mα), leading to the expression Extending this analysis of the atom dynamics to the two other spatial directions x and y leads to slightly different equilibrium temperatures along these axes, but the expected difference is less than 10%, which is below our experimental resolution (see Appendix B). According to equation (4) , the MOT temperature is minimized down to 3.4 µK for s = 2/(η−2), which corresponds to a local detuning ∆ loc = Γ/2 according to equation (2). As η 1, this value is very close to the standard Doppler limit T D = Γ/(2k B ).
We investigated this behavior by measuring the gas temperature as a function of the saturation parameter s for large laser detunings. As shown in figure 4b, the measured temperatures are reasonably well reproduced by equation (4) for 0.3 ≤ s ≤ 10. The measured temperatures deviate from theory for small saturation parameter values s 0.1, which can be explained by the shaking of the atomic sample due to the noise of the cooling laser intensity. We calculate in Appendix C the temperature increase expected from the measured intensity noise spectrum (corresponding to r.m.s. fluctuations of the saturation parameter s of 5 × 10 −3 ). We obtain the dashed curve in figure 4b, which qualitatively reproduces the temperature raise observed for small saturation parameters. A minimal temperature of 15(1) µK is measured for s = 0.2. For large saturation parameters s 10, the gas does not remain spin-polarized, and the temperature is no longer accounted for by equation (4).
We also observe a raise of the MOT temperature when increasing the atom number, similarly to previous studies on alkali and alkaline-earth atoms [31,32,23,33]. We discuss this effect in the Appendix A.

Cloud sizes and atom density
The atom density in the MOT is an important parameter to consider for efficiently loading the atoms into an optical dipole trap. In this section we characterize the cloud sizes and atom densities achieved in our setup. We measured the horizontal and vertical r.m.s. sizes σ y and σ z , respectively, using a gaussian fit of the optical density -the absorption image being taken in situ (see figure 5a). We observe that, in the spinpolarized regime, the vertical size σ z weakly varies upon an increase of the laser detuning ∆, while the horizontal size σ y increases over a larger range.
We now explain how one can account for this behavior within the simple model developed above. In the spin-polarized regime, the analytic form of the radiative force allows expressing the equilibrium shape of the atomic sample in a simple manner. Close to the equilibrium position, the radiative force can be expanded linearly as F rad = −κ x x e x − κ y y e y − κ z z e z , where the spring constants are given by where δµ = µ − µ is the difference between the magnetic moments in the excited and ground electronic states, denoted as µ and µ, respectively. Note the simple expression (5) for κ x , in which the influence of the detuning ∆ and saturation parameter s only occurs via the equilibrium position z c . We remind that the magnetic field gradient is twice larger along x than along y, which explains the relation (6) between the spring constants κ x and κ y . The r.m.s. cloud sizes σ u (u = x, y, z) are then determined using the thermodynamic relations As both the temperature T and local detuning ∆ loc are constant in the spin-polarized regime, equations (7) and (8) predict a constant r.m.s. size σ z , consistent with the weak variation observed experimentally. The variation of the size σ y is also well captured by this model. A more precise analysis would require taking into account trap anharmonicities, which cannot be completely neglected given the non-gaussian cloud shape (see figure 1c). We also observe a variation of the cloud size when increasing the atom number N with fixed MOT parameters (see figure 5b). Such an effect is expected from the repulsive interaction between atoms being dressed by the MOT light, due to the radiation pressure of fluorescence light they exert on each other [34]. We plotted in figure 5b the cloud volume V as a function of the atom number N . The volume V is defined as V = (2π) 3/2 σ x σ y σ z , so that n peak = N/V represent the atom density at the trap center. In order to calculate the volume V , we use equation (6) to deduce σ x from the measured σ y value. The measurements are consistent with a volume V independent of N for low atom numbers ('temperature-limited' regime) and linearly varying with N for large atom numbers ('density-limited' regime) [34,35]. The data is fitted with the empirical formula where Θ is the Heaviside function. For the MOT parameters corresponding to figure 5b, we obtain V single atom = 3.7(1) mm 3 , N c = 3(1) 10 7 and α = 8(1) × 10 −9 mm 3 .
Far in the density-limited regime (N N c ), we expect the atoms to organize as an ellipsoid of uniform atom density n max corresponding to the maximum atom density that can be reached in the MOT [34,35]. In such a picture, the volume determined from the r.m.s. sizes varies linearly with the atom number, with α = 0.52 n −1 max (taking into account the non-gaussian atom distribution in this regime). From the fit (9) of our data, we infer for large atom numbers a maximum atom density n max = 7(1) 10 10 cm −3 , a value comparable to the ones typically reached with alkali atoms [35].

Atom losses due to light-assisted collisions
The variation of the cloud sizes shown in figure 5 indicates that the atom density could be maximized by setting the cooling laser light close to the optical resonance, so as to achieve the smallest cloud volume. However, we observe an increased rate of atom losses near resonance, eventually leading to a reduced atom density. In order to understand this behavior, we present in this section an experimental study of atom losses in the magnetooptical trap, and we interpret the measurements with a simple model based on molecular dynamics resulting from light-induced Van-der-Waals interactions [36,37,38,39].
The loss of atoms is quantitatively characterized by measuring the variation of the atom number N with the time t spent in the magneto-optical trap, in the absence of Zeeman slowing light (see figure 6a). The atom decay is fitted with the solution of an atom loss model taking into account one-body atom losses due to collisions with the residual gas and two-body losses, described by the equatioṅ In this equation, τ 1 = 12(2) s is the one-body lifetime due to collisions with background atoms (background pressure 4 × 10 −10 mbar), β is the two-body loss coefficient and n = n peak /(2 √ 2) = N/(2 √ 2V ) is the average atom density in the trap. An example of fit of the atom number decay is presented in figure 6a .
The figure 6b shows the variation of the loss coefficient β with the detuning ∆. We observe that the loss coefficient stays almost constant in the spin-polarized regime, with β = 2.6(5) × 10 −11 cm 3 /s, and it increases by a factor ∼ 20 when approaching resonance. We can also compare our measurements to the loss coefficient β = 3.7(4) × 10 −11 cm 3 /s reported in reference [15]. As shown in figure 6b, the two measurements are in good agreement after renormalizing the laser detuning to account for the different saturation parameters used in the two studies, so as to compare atom samples with identical spin composition ¶.  (12). The red square corresponds to the decay coefficient measured in reference [15], with the arrow indicating the required detuning renormalization (see text).
We now interpret our loss coefficient measurements using a simple theoretical model in which atom losses originate from light-induced resonant Van der Waals interactions. As shown in figure 7, when an atom is brought to an excited electronic state by absorbing one photon, it experiences strong Van der Waals forces from nearby atoms. For reddetuned laser light, an atom pair is preferentially promoted to attractive molecular potentials. Once excited in such a potential the pair rapidly shrinks and each atom may gain a large amount of kinetic energy. When the molecule spontaneously emits a photon, both atoms return to the electronic ground state with an additional kinetic energy that can be large enough for the atoms to escape the MOT.
We model this phenomenon using a simple description of molecular dynamics, inspired from [36,40,39]. The large electron spin of Dysprosium leads to an intricate structure of 2(2J + 1)(2J + 1) = 646 molecular potential curves that we calculated numerically. The complete description of this complex system is out of the scope of this paper. Fortunately, the main physical effects occurring in the experiment can be captured by a simplified model corresponding to a single effective molecular potential V mol (r) = −λ Γ/(kr) 3 , with a 1/e molecule lifetime (µΓ) −1 , where λ and µ are dimensionless numbers. The mean values of these parameters averaged over the 323 attractive molecular potentials areλ 0.68 andμ 1.05.
The calculation of the two-body loss rate within this model is detailed in Appendix D. We show that the laser excitation of atoms from the spin level |J, m J to |J , m J + q (q = −1, 0 or 1) contributes to the loss coefficient as where Γ E is the Euler Gamma function. We remind that Π m J is the atom fraction in the state |J, m J and ∆ (m J →m J +q) loc is the local detuning for the considered optical transition at the MOT position. The exponential factor corresponds to the probability that a molecular association event leads to the loss of the atom pair. The total loss coefficient β is then obtained by summing the contributions β m J ,q of all optical transitions.
From equation (11), we see that the atom losses associated with this mechanism are exponentially suppressed when considering broad optical transitions. This suppression plays a large role for alkalis, for which other loss mechanisms dominate, e.g. finestructure-changing collisions [36]. In the case considered here, the exponential factor has a moderate effect: for the data shown in figure 6b, it takes the value 0.6 for the transition |J, m J = −J → |J , m J = −J in the spin-polarized regime, assuming λ =λ and µ =μ.
In the spin-polarized regime, the predominance of the transition |J, m J = −J → |J , m J = −J leads to a simpler expression for the loss coefficient As the local detuning ∆ loc does not vary with ∆ in this regime, we expect a constant loss coefficient β, as observed for the data presented in figure 6b for ∆ −2π × 1 MHz. In figure 6b we show the prediction of the full model for two sets of values for the parameters λ and µ. Using the mean valuesλ andμ only provides a qualitative description of the measured loss coefficients. A better agreement is obtained using λ = 0.75 and µ = 0.5, possibly indicating the important role played by subradiant molecular states, which correspond to µ < 1.

Conclusions and perspectives
We presented a detailed experimental study of narrow-line magneto-optical trapping of Dysprosium, together with theoretical models supporting our measurements. We showed that the optimal operation of the MOT is obtained for large laser detunings, leading to a spontaneous spin polarization of the atomic sample and to minimal two-body atom losses.
This understanding allows us to prepare gases in ideal conditions for transferring them into an optical dipole trap. In such a non-dissipative trap, it is crucial to For the duration τ spent in the excited state, the atoms attract each other, and acquire each a kinetic energy E k before returning to the electronic ground state. The atoms are lost as soon as their kinetic energy exceeds a threshold energy E * , which is typically much larger than other involved energy scales (see Appendix D). The condition for atom losses is represented as a light gray area, corresponding to 2 E * = 1000 Γ.
produce atomic samples polarized in the electronic ground state in order to avoid dipolar relaxation [41,42]. In preliminary experiments, we were able to trap about 2×10 7 atoms into an optical dipole trap created by a single laser beam of wavelength λ = 1070 nm and optical power P = 40 W, focused to a waist of 35 µm. Optimal efficiency of the dipole trap loading is obtained by first preparing the MOT at the lowest achieved temperature of 15 µK (corresponding to parameters ∆ = −2π ×1.8 MHz, s = 0.2 and G = 1.7 G/cm), and then superimposing the dipole trap over the MOT center during 600 ms. The phase space density of 8 × 10 −5 reached after the dipole trap loading corresponds to a good starting point to reach quantum degeneracy via evaporative cooling [1,3]. Our study will be of direct interest for magneto-optical traps of other atomic species featuring both narrow optical transitions and a spinful electronic state, such as the other magnetic Lanthanides. Future work could investigate sub-Doppler cooling to very low temperatures in optical molasses. Contrary to sub-Doppler mechanisms observed with magnetic Lanthanides in broad-line magneto-optical traps [43,13,44], reaching temperatures below the Doppler limit of the 626 nm optical transition would require applying an optical molasses at low magnetic field [45,46]. UQUAM) and the Idex PSL Research University (ANR-10-IDEX-0001-02 PSL ). L. S. acknowledges the support from the European Union (H2020-MSCA-IF-2014 grant n°661433).

Appendix A. Additional temperature measurements
In this appendix we describe further temperature measurements related to the equilibration dynamics and to the influence of the atom density.
We studied the equilibration dynamics in the magneto-optical trap by measuring the time evolution of the temperature right after the trap parameters have been set to the 'compressed MOT' values. We show such an evolution in figure A1a, corresponding to MOT parameters s = 0.65, ∆ = −2π × 1.84 MHz and G = 1.71 G/cm. The temperature variation is fitted with an exponential decay of 1/e time constant τ = 29(11) ms, with a baseline of 23.6(9) µK. No significant atom loss is observed over the duration of equilibration. This measurement allows extracting a damping coefficient α = 1/(2τ ) = 17(6) s −1 , comparable but smaller than the value α = 47(2) s −1 given by the simplified model leading to equation (3).
We also investigated the raise of the MOT temperature when increasing the atom number [31,32,23,33]. In previous studies, such an effect was attributed to multiple scattering of photons within the atomic sample [32], leading to a temperature raising linearly with the peak atom density n peak , as T (n) = T single atom + γ n peak .
We investigated this behavior by measuring the gas temperature for various atom densities. The atom density was varied by loading different atom numbers 4 × 10 7 ≤ N ≤ 2 × 10 8 or using different gradient values 1.1 G/cm ≤ G ≤ 5.3 G/cm. Note that the highest atom density used for this study (n peak 1.1×10 11 cm −3 ) exceeds the maximum density n max discussed in the main text as we use here larger magnetic field gradients. As shown in figure A1b, our measurements are compatible with a linear variation of the temperature with density, with a slope γ = 8(1) × 10 −11 µK cm 3 . This value is comparable with the one obtained with Cs gray molasses [32].

Appendix B. MOT temperature in the spin-polarized regime
In this section we give a more detailed description of the MOT temperature calculation and extend the analysis to the atom dynamics in three spatial directions. We restrict the discussion to the spin-polarized regime.
In order to extract the damping coefficient, the radiative force can be expanded at the equilibrium position r c as where α was introduced in the main text, see equation (3). The anisotropy of the damping comes the specific geometry of our setup (see figure 1b). The momentum diffusion tensor D is calculated as the sum of the diffusion tensors D abs and D em , associated with the stochastic absorption and spontaneous emission events, respectively. Taking into account the geometry of our experiment (see figure  1b), we obtain the diffusion coefficients The temperature is then obtained according to k B T = D/(mα) and reveals weak anisotropy: Appendix C. Temperature increase due to the laser intensity noise In this section we give more details on the calculation of the temperature increase due to fluctuations of the cooling laser intensity, leading to a time-dependent saturation parameter s(t). The main heating effect comes from fluctuations of the trap center z c (t), which can be related to s(t) using the equilibrium condition (1). For simplicity we restrict the discussion to the atom dynamics along the z axis. The atom dynamics is described by Newton's equation (C.1) Here, F L (t) is the Langevin force associated with stochastic radiative processes, such that F L (t) = 0 and F L (t)F L (t ) = 2Dδ(t − t ), involving the diffusion coefficient introduced in the main text. By integrating equation (C.1), we calculate the r.m.s. fluctuations of the velocity, as where S(ω) is the spectral density of the saturation parameter noise and ω 0 = κ z /m. The equilibrium temperature is then obtained as T = m ż 2 /k B . Note that for small damping rates α and in the absence of Langevin forces, equation (C.2) is consistent with the heating rates expected from reference [47] for shaken conservative traps. The dashed line in figure 4b is calculated using equation (C.2) and the measured noise spectrum, shown in figure C1 for the lowest saturation parameter s = 0.065 explored in figure 4b. Note that, during the MOT loading and compression, the laser intensity is servo-locked to a PID controller, typically over the range 0.065 s ≤ 50.
Van der Waals forces induce fast atom dynamics, leading to atom losses once the molecule is de-excited after spontaneous emission. A precise modelling is a challenging task, as it requires calculating the 646 molecular potentials and the corresponding excitation amplitudes, taking into account the local magnetic field value and the orientation of atom pairs. We consider here a single effective molecular potential V mol (r) = −λ Γ/(kr) 3 , of 1/e lifetime (µΓ) −1 . For simplicity we consider a uniform atom density n = N/V . The atom loss can be described by the equationṄ = − n 2 2 dr 1 dr 2 2Γ asso (|r 1 − r 2 |)P loss (|r 1 − r 2 |), (D.1) where Γ asso (r) is the rate of molecule formation for a pair of atoms of relative distance r, and P loss (r) is the probability to lose this pair of atoms after de-excitation. The factor 1 2 avoids double counting and the factor 2 accounts for the fact that each loss event corresponds to the loss of two atoms. Equation (D.1) can be recast aṡ N = −βnN, β = dr Γ asso (r)P loss (r).
Consider a pair of atoms in the MOT separated by the distance r i . We assume the rate of molecular association, triggered by the absorption of a photon (q = −1, 0 or 1 referring to σ − , π or σ + polarizations, respectively) by an atom of spin |J, m J , to be given by the standard algebra of atom-light interaction, as Γ asso (r i ) = | J , m J + q|J, m J ; 1, q | 2 µΓ 2 2s 1 + 4[∆ loc − V mol (r i )/ ] 2 /(µΓ) 2 ,(D.2) where we ignore intensity saturation effects.
Once the molecule is formed, the atom pair evolves according to Newton's law (m/2)r = −∂ r V mol (r), which can be solved implicitly. We neglect the initial atom motion as it corresponds to a weak energy scale for this problem. The electronic excitation decays after a duration τ , corresponding to a distance r f (r i , τ ), such that The acquired kinetic energy per atom E k is obtained from energy conservation, as The atom pair is lost if the acquired kinetic energy exceeds a threshold energy E * . This energy can be calculated by solving numerically the equation of motion of an atom of kinetic energy E * , initially located at the equilibrium position r c , and subjected to the radiative forces. For the MOT parameters corresponding to the data of atom number represented in figure 6a, we estimate a capture velocity v * of about 0.6 m/s, corresponding to an energy E * 500 Γ.
The loss probability is then obtained as the probability for the atoms to acquire enough kinetic energy during the molecular dynamics: where µΓe −µΓτ is the density probability for the spontaneous emission to occur at time τ (Θ is the Heaviside function). We obtain As the threshold energy E * is much larger than other energy scales, we can safely replace the factor f [·] by f (0) = √ πΓ E (5/6)/[3 Γ E (4/3)]. The loss rate can then be obtained by calculating numerically the integral (D.1). We find that replacing the Lorentz absorption profile (D.2) by a strict resonance condition (i.e. the suitable Dirac δ function) only introduces minor differences for the numerical value of β. After this replacement the integral can be calculated analytically, leading to the formula (11) in the main text.