Suppressing dipolar relaxation in thin layers of dysprosium atoms

The dipolar interaction can be attractive or repulsive, depending on the position and orientation of the dipoles. Constraining atoms to a plane with their magnetic moment aligned perpendicularly leads to a largely side-by-side repulsion and generates a dipolar barrier which prevents atoms from approaching each other. We show experimentally and theoretically how this can suppress dipolar relaxation, the dominant loss process in spin mixtures of highly magnetic atoms. Using dysprosium, we observe an order of magnitude reduction in the relaxation rate constant, and another factor of ten is within reach based on the models which we have validated with our experimental study. The loss suppression opens up many new possibilities for quantum simulations with spin mixtures of highly magnetic atoms.


Introduction
Experiments with ultracold atoms or molecules are often limited by unfavorable inelastic collision rates.Several methods have been developed to control collisions such as isolating atoms in deep lattices [1], reducing collisional channels via confinement [2], or by mitigating their effects through the enhancement of elastic collisions via Feshbach resonances [3].Polar molecules, in particular, have been shielded from chemical reactions at short range by using repulsive interactions between electric dipoles, either in two dimensions or via microwave dressing [4,5].

Basic principles
The dipole-dipole interaction is attractive in the case of a tip-to-tail orientation and repulsive for the side-by-side one.Constraining atoms to an xy plane, with a magnetic moment aligned perpendicularly along z, leads to a largely side-by-side repulsion and generates a dipolar barrier.The dipolar length 1 a dd = µ0 4π µ(10µ B ) 2 2 represents the strength of the interaction and the two-particle oscillator length a z = /µω z the extension of the cloud in the z direction.We denoted µ and ω z the reduced mass and the trap frequency respectively.A dipolar barrier appears when the dipolar length a dd > 0.34 a z [16], which we refer to as the quasi-2D regime.Thus, experiments with dysprosium require 10,000 times higher axial frequencies than polar molecules to compensate for the 100 times smaller dipolar length.Our experiments have reached this regime with a z = 20 nm and a dd = 10 nm.
Three parameters determine the loss rate in quasi-2D: the ratio a dd /a z set by the confinement, the temperature T , and the magnetic field B. The potential barrier increases with confinement, ultimately reaching the pure-2D limit as a z → 0 as shown Fig. 1a.As the temperature decreases, the wave function of an incoming pair is suppressed by the barrier over a longer range, thereby decreasing the chance of two atoms reaching close range.This shielding effect on the wave function is illustrated in Fig. 1b.As the magnetic field increases, the range where dipolar relaxation occurs is shortened and the shielding increases.Indeed, a higher magnetic field leads to a higher released energy, and correspondingly a more rapidly oscillating outgoing wave function (see red curve Fig. 1b).Since the dipolar potential falls off as 1/r 3 , the majority of the decay will come from the first oscillating lobe of the outgoing wave function, as seen in Fig. 1c.The range of dipolar relaxation, therefore, decreases as the magnetic field increases.This can also be explained in a semi-classical picture: the Franck-Condon principle predicts spin flips to occur at the classical turning point of the outgoing wave function [2,17], i.e. when the released Zeeman energy equals the energy of the centrifugal barrier.Correspondingly, higher magnetic fields cause spin-flips to occur at a shorter range, ultimately behind the barrier felt by the incoming atoms, where they are strongly suppressed.Therefore, shielding qualitatively changes the magnetic field dependence of the dipolar relaxation rate.For bosons, in both 3D [2,15,18] and unshielded 2D geometries, the relaxation rate increases with magnetic field.When accounting for the barrier in 2D, however, the signature of dipolar shielding appears: the relaxation rate decreases with magnetic field (see Supplementary Information).

Experiment
Here we study these principles experimentally.We load ∼ 8×10 4 spin-polarized 162 Dy atoms in the excited |J = 8, m J = 8 Zeeman level (see Methods for details) in an optical lattice and get a stack of about 45 thin pancakes ('crêpes').The crêpes reach an a z /2 = 10 nm root-mean-square (RMS) width and a 5.7 µm radius.The peak density is 2.9 × 10 9 cm −2 .The experiment is performed at T ≈ 1.6 µK, above the BEC transition temperature (300 nK), to prevent convolving our results with changes in the two-particle correlation function [19,20].The quantization axis is set by an external magnetic field along the z direction.The lattice beam is blue detuned, with its radial repulsion compensated by a coaxial red-detuned optical dipole trap, as shown in Fig. 2a.Axial trap frequencies are limited to ω z /2π = 260 kHz by the maximum laser power of the compensation beam.
By measuring the atom losses we determine the inelastic decay coefficient, β 3D , as defined by the differential equation for the 3D density n: We obtain densities from the measured atom number, temperature and trap frequencies, and average over the stack of crêpes (also see the Methods section).We sometimes refer to the 2D loss rate β 2D in cm 2 /s, which uses the 2D ) Effective radial potential between two atoms from equation ( 6) for: no confinement (light blue, 3D), quasi-2D with ωz/2π = 300 kHz (steel blue) and pure-2D (dark blue).The incoming energy is given by the temperature T = 1 µK.(b.) Wave function solutions of equation ( 6) with n = 0 and initial orbital momentum m i = 0 for the three confinement strengths described above, and in red the spinflipped outgoing wave function to m f = 2 and B = 500 mG.The effect of shielding of the outgoing wave function is negligible for these parameters.(c.) Integrand of Fermi's golden rule (equation ( 2) and see equation (S21) in Supplementary Information).Each curve is the product of the respective wave function in (b.), the outgoing wave function, and the double spin-flip operator from equation ( 4) integrated with the harmonic oscillator wave functions in the z-direction.The shielding we implement here corresponds to the difference between the light blue and steel blue curves.The minimum attainable decay rate for this incoming energy corresponds to the dark blue curve.See Supplementary Information for insights on the behavior of the integrand.
density instead.It is related to β 3D through the axial harmonic confinement via β 2D = β 3D /(a z √ π).Our experimental results are shown in Figs.2b-c.We also compare the theoretical shielded decay rate (solid blue) with the one we would expect in the same crêpe geometry if there was no elastic dipolar potential to repel the atoms (dashed blue).In contrast to the loss rate in a 3D geometry (red), which increases with √ B (see Supplementary Information for comments about this scaling), we observe the signature of shielding in Fig. 2b: a much weaker dependence on magnetic fields (solid blue).
1064 red beam 741 blue lattice Fig. 2 Experiment scheme and results.(a.) Trap geometry.A blue-detuned 741 nm retroreflected beam repels the atoms to create a 1D lattice.The finite contrast of the lattice and the zero-point motion of the atoms in the ground state create a repulsive transverse potential, which is compensated by a 1064 nm red-detuned beam to create an adjustable transverse harmonic confinement.(b. and c.) Experimentally measured β 3D in a large volume trap (red) and in a thin layer (blue).The lines are theory curves obtained by using Fermi's golden rule (see Supplementary Information for derivations).The red curve shows the decay rate in 3D [15,18], the dashed blue curve is for non-shielded atoms in a lattice [2].The solid blue line takes into account the shielding induced by the elastic dipole-dipole interaction.All theoretical curves are thermally averaged over the incoming momenta.The shaded blue region corresponds to the inclusion of van der Waals contact interactions (see Supplementary Information).(b.) Measurement of β 3D as a function of magnetic field.The axial trap frequency is ωz/2π = 185 kHz which corresponds to az/2 = 13 nm.(c.) Measurement of β 3D in a constant magnetic field of 200 mG while varying the trap frequency.The uncertainties are set by the atom number stability, cloud temperature measurement and trap frequency measurements (see Methods section).
We operate in the quasi-2D regime which differs from the pure-2D one in several aspects.Compared to pure-2D, the finite axial extent of the quasi-2D geometry softens the radial barrier, reducing the barrier height to energies comparable to typical temperatures in the experiment.Furthermore, for Zeeman energies that are larger than the axial trapping frequency, new collisional channels open, with a portion of the released energy converted into axial excitation and the remainder into radial motion.As a result, the relaxation for these processes is shifted to larger distances, thereby weakening the shielding.The first channel opening is visible in Fig. 2b around 100 mG as well as in Fig. 3a-b.The aforementioned factors lead to a relaxation rate that does not decrease with magnetic field, as it would in the pure-2D case, but instead shows a weaker increase compared to the case without a dipolar barrier (dashed blue).Fig. 2c shows the loss rate coefficient as a function of axial confinement.The loss rate decreases with confinement due to enhanced shielding by the dipolar repulsion and closing axially excited states channels.
We have reduced the loss rate coefficient to approximately 1×10 −12 cm 3 /s.Over a large range of magnetic fields in a lattice, we achieved more than an order of magnitude reduction in the dipolar relaxation rate coefficient compared to the unshielded case.The agreement between the numerical calculations and the experiment enables extrapolation beyond the current limitation of the experiment: very favorable loss rate coefficients of 2 × 10 −13 cm 3 /s can be achieved at 200 mG with an axial confinement of 500 kHz at 1 µK.This matches the lowest rate obtained with fermions through Pauli suppression in reference [15].Under such conditions, axial excitations are energetically forbidden and the 2D decay rate is less than a factor of 3 above the pure-2D limit.By lowering the temperature to 100 nK, the relaxation rate would be suppressed by an additional factor of three and reach the 10 −14 cm 3 /s regime.To further understand how these numbers are computed, we describe our theoretical model in the following paragraphs.

Theoretical model
Dipolar relaxation rates can be calculated from Fermi's golden rule.The decay rate Γ of 2 particles is given by where ρ(E) is the final density of states at energy E. The incoming wave function is an excited Zeeman state with transverse momentum k i in the lowest harmonic oscillator state, n i = 0.The outgoing wave function is a lower Zeeman state with momentum k f in the harmonic oscillator state n f .The loss rate coefficient β 2D is related to Γ through β 2D = πL 2 Γ, with L being the radius of the transverse box used to normalize the wave functions.The atoms are coupled by the magnetic dipole-dipole interaction: where r is the interatomic separation (with corresponding unit vector u r ).
The magnetic field points along z.
with z = z/r and r+ = (x + iy)/r.Equation (4) shows the three effects of the dipolar interaction: an elastic scattering process, a single spin-flip proportional to z, and a double spin-flip which implicitly depends on z through r.In the two-dimensional limit where z = 0, the single spin-flip term vanishes and the elastic term is a purely repulsive 1/ρ 3 potential (where ρ = x 2 + y 2 ).This potential has an analytic solution at zero temperature (k i = 0) [21], while other cases have to be solved numerically.
We assume a quasi-2D geometry where we ignore the effect of V dd,0 on the z motion, which is then factorized and described by harmonic oscillator wave functions (see Methods for a discussion on this approximation).The elastic portion of the operator in equation ( 4) is averaged over the z direction.This leads to an effective repulsive potential (see Fig. 1a) in the one-dimensional radial Schrödinger equation: Here, the state |n is the n th harmonic oscillator's state along z.We focus on incoming states with zero projection of orbital angular momentum, m i = 0, as this channel dominates for any reasonable magnetic field (see Supplementary Information).We solve the Schrödinger equation for the radial wave function using numerical techniques, and use it to perturbatively calculate the dipolar relaxation rate with Fermi's golden rule (2).In Fig. 1 we show how dipolar repulsion (Fig. 1a) modifies the incoming wave function (Fig. 1b) and reduces the integral of the transition matrix element (Fig. 1c).
Without axial excitation, only double spin flips to the final spin state |j 2 = |7, 7 and orbital state m f = 2 are allowed.At sufficiently high magnetic field the energy released during the collision can exceed ω z , thereby opening up new collisional channels resulting in axial excitations.Energy conservation requires The single spin-flip channel (∆j = 1) requires odd ∆n due to the odd symmetry of the z term in equation ( 4), whereas double spin flips (∆j = 2) require even ∆n.Newly opened channels increase the decay rate, as shown in Fig. 3a-b.Furthermore, as previously explained, they also decrease the shielding factor, as visible in the small notch in Fig. 3c.
Remaining in the ground state of the harmonic oscillator is therefore necessary for obtaining extremely low relaxation rates, but that requires working at low enough fields.Unfortunately, the relaxation rates we measure at very low fields deviate from the theoretical values in Fig. 2b, most likely because of imperfect circular polarization of the lattice and compensating beams.The mixture of σ + and σ − light induces Raman couplings between |m J = +8 and other even |m J states, thereby opening additional relaxation channels via spin exchange [18].With a > 95% circular polarization purity, we find agreement between experimental decay rates and calculated dipolar relaxation rates for fields > 100 mG, where the Raman coupling is suppressed by Zeeman detuning.(c-d) The suppression factor defined as the ratio of β 2D obtained from shielded and free wave functions, both at fixed temperature T 1 (c), and at fixed magnetic field B 1 (d).For each of the graphs we present curves for ω 1 /2π = 300 kHz (dotted), ω 2 /2π = 1.8 MHz (dashed dotted), the pure-2D case (solid line) as well as the analytical approximation (grey dotted) from equation (S32) detailed in Supplementary Information.

Discussion and Outlook
We have shown that confinement in thin layers not only reduces the number of available collisional channels, but additionally provides dipolar shielding, thereby strongly suppressing dipolar relaxation between atoms.In principle, arbitrarily low loss rates and infinite shielding factors are possible at very low temperatures.Strong magnetic fields are also predicted to reduce the shielded collision rate to arbitrary low values if strong axial confinement suppresses the opening of collision channels.As we have discussed above, rather straightforward improvements in axial confinement, purity of polarization and temperature should result in rate coefficients in the 10 −14 cm 3 /s regime.
Our simulations and experiments show that there is already substantial shielding at thermal energies comparable to the barrier height.Lowering the temperature well below the barrier eventually results in exponential suppression [22].For our experimental parameters, going from 1 µK to 100 nK would increase the suppression by a factor of three.
In this work, we have discussed the interplay between the elastic and inelastic aspects of dipolar interactions.Both scale with the dipolar length, which could be 10,000 times larger for polar molecules.Yet the large total angular momentum J = 8 works in favors of dysprosium over molecules, as the elastic part of the dipole-dipole potential scales as J 4 in a stretched state, while the relaxation rate scales as J 3 for single spin-flips and J 2 for double spin-flips.
An important point of comparison is the elastic scattering rate.At 1 Gauss in a trap with a 2 MHz axial frequency, the inelastic 2D cross-section would be 20 nm without shielding.Shielding drops this number to 0.3 nm, while the semi-classical dipolar elastic collisional cross section is σ SC = 180 nm [21].
Shielding is necessary to obtain a ratio of good to bad collisions in excess of 100.
Dipolar shielding has previously been observed in polar molecules with fermionic statistics [4], for which the shielding is qualitatively different.Since identical fermions already have an isotropic p-wave barrier, adding moderate dipolar interactions in a confined geometry will first strengthen this barrier in the radial direction but also weaken it in the axial one.As a result, the inelastic collision rate will first decrease with the dipole moment and then increase [23].This cannot be seen with bosons.For both particle types, the inelastic collision rate will eventually decrease when entering more deeply into the 2D regime, as we have explored in this work.Our technique would be crucial to study spin mixtures of bulk gases of bosonic dipolar species.
In conclusion, we have demonstrated a way to realize long-lived spin mixtures in dense bosonic lanthanide clouds, opening up new possibilities for quantum simulation experiments in two dimensions.With such technique, dysprosium can be used to study quantum materials with dipolar interactions in regimes different from those currently possible for polar molecules [24] and Rydberg atoms [25].Stable spin mixtures are important for implementing spinorbit coupling and artificial gauge potentials via Raman coupling of spin states [26,27].By suppressing dipolar relaxation, one can take advantage of the ground state orbital angular momentum of lanthanides to avoid the substantial photon scattering rates of the Raman beams for alkali atoms [28].

Sample preparation
We prepare spin-polarized samples of ∼ 8 × 10 4 162 Dy atoms in the |J = 8, m J = −8 state in an optical dipole trap just above the transition temperature.The samples are obtained after evaporative cooling in a crossed optical-dipole trap (ODT) which is loaded from the narrow-line magnetooptical trap described in reference [S29].Working with a thermal gas makes it easier to determine dipolar relaxation rate coefficients without accounting for a varying condensate fraction.
The highest spin state |m J = +8 is populated via adiabatic rapid passage using an RF sweep in a magnetic field of 3.5 G along the z direction.A stack of quasi-2D layers, which we refer to as crêpes due to their extreme aspect ratio, is created using a 1D optical lattice formed by retroreflecting a 741 nm laser beam along the z axis.The lattice beam is blue-detuned from the 741 nm transition (which has a linewidth Γ/2π = 1.8 kHz) by several GHz, thus providing frequency-controllable tight axial confinement.A coaxial vertical optical dipole trap is used to compensate for the transverse repulsion resulting from the blue-detuned lattice.The lattice and the vertical dipole trap are turned on using exponential ramps with a 50 ms time constant to adiabatically load the atoms into the lowest vibrational level of the 2D layers.During the first 40 ms of the lattice ramp, the magnetic field is rapidly reduced to 40 mG to minimize the dipolar relaxation losses.The magnetic field is then ramped up to its final value during the last 10 ms of the lattice loading ramp, after which the decay of the sample due to inelastic collisions is measured.

Zeroing the magnetic field
Achieving control of low magnetic fields is critical for minimizing dipolar suppression by preventing higher outgoing vibrational channels from opening.We have devised a method to zero the magnetic field that relies on the large disparity of Clebsch-Gordan coefficients for dysprosium.When an atom's magnetic moment is aligned along the propagation of a circularly polarized imaging beam, the amount of scattered light strongly differs whether the magnetic dipole moment is oriented parallel or anti-parallel to the propagation of the imaging beam.By using absorption imaging for various external magnetic fields, as shown in Fig. S1, one can observe when the dipole moment has flipped, which determines the zero of the external magnetic field.
More specifically, in a spin-polarized (m J = −8) sample of bosonic dysprosium, the Clebsch-Gordan coefficients for σ − , π and σ + transitions are 1, 1/9 and 1/153 respectively.We perform absorption imaging of a spin-polarized sample with left-circularly polarized (σ L ) light along the magnetic field quantization axis z.We work with low enough light intensity and imaging time to prevent optical pumping.At large positive magnetic field bias, the atoms see σ − light with a corresponding Clebsch-Gordan coefficient of 1, resulting in a large atom count.At large negative magnetic field bias, the atoms see σ + light with a corresponding Clebsch-Gordan coefficient of 1/153 leading to a low atom count.The lower the transverse magnetic field, the sharper is the transition when the longitudinal field is varied.In this way, the zero settings for all components of the magnetic field are determined.

Lattice light choice
The need for deep optical lattices requires a tightly focused lattice beam, which causes undesirably strong radial confinement if one uses a red-detuned beam.By choosing a blue-detuned lattice we avoid adiabatic compression of the cloud in the transverse direction and the substantial corresponding increase in temperature when ramping up the optical lattice.The choice of a bluedetuned lattice also exposes the atoms to lower light intensities and reduces the unwanted Raman transitions due to imperfect circular polarization.However, the radial deconfinement created by the lattice needs to be compensated, which we achieve with a red-detuned optical dipole trap that enables independent control of the axial and transverse trap frequencies (see Fig. S2 left).
The lattice was created by near-resonant 741 nm light from a Ti:Sapph laser which can deliver about 300 mW of light to the atoms, after fiber coupling and intensity stabilization.

Trap geometry
Atoms are loaded into an optical dipole trap consisting of three 1064 nm laser beams: two beams with 40 µm beam waists crossed at 8 • in the horizontal plane, and a beam with a 64 µm waist propagating along the (vertical) z direction.During the dipolar relaxation experiment, the horizontal beams are switched off, and the vertical beam serves to compensate for the deconfinement of the blue-detuned lattice.The lattice beam is focused down to a waist of 50 µm.It is typically detuned by 14.25 to 2.25 GHz to the blue side of the narrow 1.8 kHz transition [S30].The transverse antitrapping potential is compensated using 8 W in the vertical trapping beam.We verified with in-situ images (obtained with detuned imaging light due to the high optical densities) that the blue-detuned lattice is correctly compensated without displacement of the cloud.
The RMS extension of the cloud along the lattice direction before loading is σ ODT 4.7 µm.Given the layer separation of λ/2 371 nm, around 4 √ πσ ODT /λ = 45 crêpes are loaded with initially 3 × 10 4 atoms and a central density of n 0 = 2.9 × 10 9 cm −2 .The density distribution in the i th pancake is described by (see Fig. S2) with

Lifetime analysis
The decay of the cloud can be described via equation (1) for the 3D densities ) or by using a 2D equation The densities in each pancake are related by with σ z = 4µωz = a z /2 10 nm.When integrating equation ( S2) and equating it to (S3), we obtain In the main paper, we are using β 3D to characterize the decay.We will omit the 2D subscript for the densities in the rest of the manuscript.Equation ( S3) -in a local-density approximation -needs to be integrated over the cloud volume to relate to the observed quantity N , the number of atoms: The effective volume V eff is determined as follows.After integration of equation ( S1), one gets Identifying V eff in equation (S6) gives We note that n = N 2 3/2 V eff , where each √ 2 factor comes from the Gaussian averaging along one axis.To take into account the moderate heating during the experiment, we perform a linear fit of the temperature T (t) = T 0 + v T t which is used to scale the effective volume V eff (t).The solution of the differential equation ( S6) that we fit is from which we determine β 2D .The atom number N (t) is measured as a function of hold time (typically tens of ms) using time of flight imaging.
Here we have assumed that every dipolar relaxation event leads to the loss of both atoms.This is justified since the effective trap depth of a few microkelvins is negligible compared to the kinetic energy gained by the spin-flip for magnetic fields larger than a few tens of milligauss.The experiment is sufficiently fast (tens of ms) such that photon scattering in the lattice, background collisions and residual evaporation are not important.

Error bars
The uncertainties represented by the errorbars in the plots come from the statistical error due to the curve fitting, as well as our best estimate in the uncertainties of ω ODT , T ODT , ω ⊥ , ω z and T lattice .ω ODT is measured in the ODT by observing the oscillation of the cloud after suddenly applying a magnetic force.ω ⊥ is measured by modulating the intensity of the vertical trapping beam and observing the parametric heating resonance.Axial trap frequencies ω z (which go up to several hundreds of kHz) are also measured via parametric heating.Due to the bandwidth of the drive electronics, this was done only in shallow lattices and extrapolated to deeper lattices.Temperatures are observed in time of flight.

Theoretical decay rate
We summarize here the main equations to produce the theoretical predictions in Figs. 1, 2 and 3. Detailed derivations are given in Supplementary Information.
We recall equation ( 4) which is used to compute the potential used in equation ( 6): ) This is the equation that we solve numerically for both the incoming (m = 0, n = 0) and outgoing (m = 2, n even or m = 1, n odd) wave functions.The code we developed combines grids of multiple step-sizes to account for the need to appropriately average the potential along z, describe the shortrange shielding at small ρ and normalize correctly the wave functions at large distances.Given the temperature, magnetic fields, z trapping frequencies and desired precision, the code determines an appropriate grid, and computes the incoming wave function and the harmonic oscillator states on this specific grid.
It then distributes those results on multiple cores, computing the outgoing wave function for each of the different decay channels and the respective integral of Fermi's golden rule.This method enables the code to produce the plots presented in this paper on a simple laptop in a reasonable time.
The normalization condition reads: L 0 dρ φ n,m (ρ) 2 = 1 for a cylinder of radius L. Our model accounts for the modification of both incoming and outgoing wave functions by the dipolar interaction.The free radial wave function solution with momentum k is φ The 2D loss rate coefficient for the channel 2 with χ n being the n th harmonic oscillator's state wave function: The total rate is then the sum over all channels: .
and relates to the 3D rate as β 3D = √ πa z β 2D .The rate is eventually averaged over the thermal distribution of incoming momenta (see Supplementary Information) for the Fig. 2 b-c, and computed at the mean momentum for all of the other figures.

Pure-2D limit
In pure-2D the double spin-flip potential is V dd,2 (ρ) ∝ 1/ρ 3 .If we ignore the shielding, in the low-temperature limit, we find that the 2D decay rate is If we incorporate shielding we find that which goes to zero at zero temperature.Under certains assumptions detailed in Supplementary Information and noting x 2 the first zero of the Bessel function J 2 (x), one can find that the decay rate in a certain field range behaves as ) which vanishes at high magnetic fields.

Discussion of various approximations
Unitarity limit.The perturbative results will get modified when the decay rate approaches the unitary limit.However, in our range of parameters, the decay rates are much smaller than the unitary limit.A β 3D of high 10 −12 cm 3 /s corresponds to a β 2D = σ k/µ in the low 10 −6 cm2 /s.This gives a σk 0.1 4 which puts us safely in the non-unitary regime 2 .
Wave function substitution approximation.The system is perturbed by two part of the dipolar potential: one which is diagonal in the spin states basis and therefore elastic, the other part is non-diagonal and causes transitions.Usually, the weakest part of the Hamiltonian should be treated perturbatively.It is the case here since V inelastic dd / V elastic dd 2 1/J = 1/8 as shown in equation ( S7).This is why we evaluate the decay rate on the shielded wave function.Following previous treatments [S2, S15, S18], we have not checked the importance of higher-order terms in the perturbation theory.
Effective potential approximation.To compute the wave functions we assumed an effective potential obtained by averaging V dd in the n th state of the harmonic oscillator.This is the diabatic limit of a coupled channels calculation.There exists a fully adiabatic method to compute the molecular potential of two interacting dipoles in a quasi-2D geometry [S16].It would mix the harmonic oscillator states but we found it would only affect the wave function at short distances, which is important only at high magnetic fields.In our experiment, kf remains on the order of 1, and restricting the z-motion to the pre-existing harmonic oscillator states is acceptable.
Fermi's golden rule approximation.The use of Fermi's golden rule with the original density distribution is valid only if the decay rate is smaller than the other time constants of the system.The relaxation rate Γ = β 3D n 10 2 s −1 is indeed smaller than the collision rate which is around 10 3 s −1 , or the trap frequencies of 200 Hz.This assumption is therefore fulfilled.The system will stay in (quasi-) equilibrium when the loss rates are smaller than the trapping frequencies and smaller than the rate of elastic collisions which provides thermalization.
Neglecting molecular potentials.The background s-wave scattering length of dysprosium a = 5.9 nm can modify the wave functions.A previous paper [S2] studied extensively its influence on chromium.However, the decay rate we observed in a large volume 3D trap (red curve in Fig. 2) agrees better with the theory which does not take the scattering length into account.Another dysprosium experiment [S15] found a similar result in an even wider range of fields.However, it is possible that the short-range molecular potential plays a role in the 2D results, and could possibly explain why we obtain rates a few times smaller than the theory predicts (see shaded areas in Fig. 2).Indeed, a sizeable contribution to the loss comes from interatomic distances smaller than the van der Waals length a vdW = 4.3 nm and the scattering length a = 5.9 nm.Since the real wave function rapidly oscillates at short-range, the contribution to the overlap matrix element should vanish.

Supplementary Information
We present here our theoretical model to compute the expected rates shown in Fig. 2. The calculations use the formalism of Fermi's golden rule.For completeness, we also derive and discuss the Born approximation previously sketched in [S2] and used extensively in [S31].We comment on their equivalence and then discuss the pure-2D limit.Both Fermi's golden rule and the Born approximation are perturbative approaches to quantum scattering.The former decomposes incoming and outgoing waves into many different channels, reducing the calculation effectively to 1D (i.e. the radial coordinate).The latter computes the perturbative impact of the dipolar relaxation Hamiltonian on an incoming plane wave, leaving the calculation in 2D.Both results are equivalent.It is easier to generalize Fermi's golden rule for our approach, which uses initial wave functions modified by the dipolar elastic potential.The Born approximation will be presented only for un-modified plane waves, whereas Fermi's golden rule formalism will be applied to both modified and un-modified wave functions.

Dipolar units
For the following derivations, it is convenient to use dipolar units, indicated with a tilde.Dimensionless lengths are set with the dipolar length, such that ãz = a z /a dd , wave functions become φ = a 1/2 dd φ, momenta become ki = k i a dd and energies are measured in units of dipolar energy E dd = 2 2µa 2 dd , so Ṽdd = V dd /E dd .The equation (4) becomes and equation ( 6) for an arbitrary harmonic oscillator channel n becomes The normalization condition reads: L 0 dρ φn,m (ρ) 2 = 1 for a cylinder of radius L = La dd .Our model accounts for the modification of both incoming and outgoing wave functions by the dipolar interaction.The free radial wave function solution with momentum k is φ(free) n,m (ρ) = π k ρ L J m ( k ρ), which does not depend on n.

Fermi's golden rule derivation
Here we derive the expression for the 3D loss rate coefficient for the channel |j 0 → |j f , |0 → |n f : χn are the harmonic oscillator wave functions in the z direction.ãz is the harmonic oscillator length in units of the dipolar length.

Fermi's golden rule
Fermi's golden rule conveniently expresses the decay rate from an initial state to a continuum having a certain density of states.It is natural to compute Fermi's golden rule with free incoming and outgoing plane waves, which we do at first.We then expand it into cylindrical waves, which we finally substitute for their shielded version solutions of equation ( S13).The starting point is using equation ( 2) to compute the decay rate Γ j f ,n f plane of a pair of polarized atoms coming in a plane wave in the ground state of the harmonic oscillator, and outgoing in another asymptotic plane wave in the oscillator state n f with a spin state |j f .The particles are assumed to be contained in a cylinder of length L and by a harmonic oscillator potential in the z direction.
with the total wave function |Ψ j,n = | k ⊗ |n ⊗ |j and

Density of states
In a two-dimensional box, the volumic density of states is µ 2π 2 .In a particular direction of angle dθ k , the density of states is ρ p (k f ) = µL 2 4π 2 dθ k .

Plane wave expansion
To fully use the symmetries of the dipolar potential, one can expand the plane wave into spherical waves: which have the following position representation: ) Summing over all possible outgoing directions, the rate is

Dipolar interaction
To simplify the problem we can look at the selection rules of the dipolar potential.The equation ( 3) can be written which gives equation ( 4) with r+ = ρ r e iθ and z = z/r.Hence V dd, 0 is independent of θ, V dd, 1 proportional to e iθ and V dd, 2 to e 2iθ .The e imθ in equation ( S15) makes the matrix elements k f , m f | Vdd, j f |k i , m i of the sum (S16) non-zero only if m f = m i + j f .Furthermore, as z is anti-symmetric, and both χ 0 and r+ are symmetric, V dd, 1 can only promote to odd n f states and V dd, 2 to even ones.This gives

Symmetrization
Since the atoms are bosons, the wave functions need to be symmetrized.All the spin states are already symmetrized.Each incoming and outgoing . It transforms the sum (S19) by multiplying the incoming and outgoing terms by √ 2 each and summing on even m i for bosons and odd for fermions.The density of states of the outgoing channels is divided by 2 to avoid double counting.Furthermore, due to the independence of the braket on θ f , the only terms in the sum giving a non-zero contribution after integrating over θ f are the terms diagonal in m i coming from the modulus.This gives: ) s-wave scattering Given our parameter range, we only keep the m i = 0 channel, i.e. the swave channel.This reflects that the relevant range of the dipolar potential is much smaller than the incoming De Broglie wavelength.It comes from the fact that k f k i for most of the magnetic fields (see in the "Comparing Born approximation and Fermi's golden rule" section of the Supplementary Information for a discussion when the magnetic energy is comparable to the temperature).The outgoing wave function starts to oscillate at a distance ∼ |m f |/k f ∼ 1/k f which cuts out the integration at this distance (see Fig. S6).Before that point, the free wave function φ m (ρ) rises like a Bessel function in √ ρ(kρ) |m| , and the dipolar potential goes as 1/ρ 3 .The integral goes then like , such that all the terms for which m i = 0 are greatly suppressed as soon as the magnetic field energy is greater than the temperature (1µK, which is about 10 mG).Therefore

Wave function substitution
We have expanded the plane wave in equation ( S14) into cylindrical wave functions, solutions of the Schrödinger equation for a free particle, i.e. equation (S13) without the dipolar interaction term.We now replace these free cylindrical wave functions φ m by the solutions of the full Schrödinger equation with the dipolar interaction term φ n,m , which now depend on the harmonic oscillator channel n.The initial state is now the shielded wave function, and Fermi's golden rule describes its dipolar decay.This substitution is valid as the plane wave expansion (S14) still holds at large distance since the dipolar potential in 1/ρ 3 decays faster than the centrifugal 1/ρ 2 potential.

Calculation of β 2D
The rate Γ sym is the probability per unit of time of the two bosons decaying when placed in a harmonic oscillator state n = 0 in a cylindrical box of radius L. We are interested in the decay rate β 2D defined through dn2D dt = −β 2D n 2 2D , which for a homogeneous gas gets integrated into dN dt = −β 2D N 2 πL 2 .For N atoms homogeneously spread in the area, the differential equation sums the decay rate on all the possible pair combinations N (N − 1)/2.For each event two atoms get lost, so The φ n,j wave functions are either free wave functions like the light blue curve in Fig. 1b or modified wave functions through equation ( 6) such as the blue and navy curves on the same figure.

Calculation of β 3D
The final equation for β 3D follows from equation (S5), so S21) Note that the prefactor E dd a 3 dd / is useful to check the units but inconvenient to look at the scaling with the dipolar interaction.When the wave functions are not modified by the dipole-dipole potential they read: φn,m (ρ) = π k ρ L J m ( k ρ) z , which gives for the 2D rate: which explicitly shows the a 2 dd ∝ (10µ B ) 4 scaling.

Total decay rate
The total rate is then the sum over all channels: .

Momentum averaging
The rate obtained depends on the incoming momentum ki .The gas being thermal with many occupied states in the transverse direction, we integrate over momentum to obtain the average decay rate The results presented in Fig. 2b-c are the average momentum rates.But since this is rather computationally heavy and blurs the channel opening, we only present results computed at the mean momentum ki = κ √ π/2 in the rest of the paper.For instance, the incoming energy of the wave functions in Fig. 1a-b

Van der Waals interaction
Short-range van der Waals interactions exist on top of the dipolar ones.To get a sense of the sensitivity of our model to that contact interaction we also used simulated wave functions with a hard-core potential at a s = 5.9 nm, the effective background s-wave scattering length of dysprosium [S32], while keeping the dipolar potential elsewhere.We put a node in the incoming and outgoing radial wave functions φ at this position, and integrated from this distance outward.This produced the lower bound of the shaded area in Fig. 2.

Definitions
The Born approximation is the standard way to describe scattering by a potential.The method in 2D has been described in [S2], but it is lacking of an explicit final formula, which we would like to present here.We look for eigenstates of the full Hamiltonian: ) and Vdd defined through equation (S17).We define the Green operator Ĝ0 as (E − Ĥ0 ± iη) Ĝ0 = 1 and will eventually take the limit η → 0. We omit the surface normalization coefficient in the wave function such that: 3 It follows that:

Born approximation
Following the definitions, |Ψ = Ĝ0 Vdd |Ψ is one solution of the eigenvalue equation (S22).If furthermore we have a solution |Ψ 0 to the equation (E − Ĥ0 ) |Ψ 0 = 0 then |Ψ = |Ψ 0 + Ĝ0 Vdd |Ψ is also a solution.The first Born approximation consists in computing the first-order part of the solution: Let us take |Ψ 0 = | k i ⊗ |n = 0 ⊗ |j 0 which has the energy: . Note that we will later symmetrize the incoming state by replacing , so it is good to keep track of the e i ki• ρ terms that will eventually become We decompose the problem into channels (n f , j f ).We eventually want the position representation of the scattered wave function but the Green operator has a simpler expression in momentum space.When projecting |Ψ on j f | n f | k| the first term in equation ( S24) disappears for (n f , j f ) = (n i = 0, j 0 ).The matrix element to compute is thus When acted on the left side, the Green operator passes through the states which are the eigenstates of Ĥ0 from equation (S23), giving S25) with ∆m = 1, 2 for channels j 1 , j 2 respectively.By setting: which uses the notation of equation ( S17) for the potential.
Calculation of Ψ n f ,j f ( ρ) In position representation we have Computing the integral I requires a few steps in Mathematica: with K 0 being the Bessel K function with a complex argument and H (1) 0 the Hankel function of the first kind.We only keep the + solution from the square root in the Hankel function to have an outgoing flux and perform the far field expansion To further simplify we first introduce the Fourier transform H n f of the harmonic oscillator wave functions product χ * n f (z)χ 0 (z): and then define q = k f − k i + q z u z and r 1 = ρ 1 + z u z which gives differential scattering cross section into an angle dθ f is So the total scattering cross section is The loss coefficient is β = σv = σ ki µ .It still depends on the direction of k i .It is then averaged in all possible incoming directions: Comparing Born approximation and Fermi's golden rule  The regions shaded in gray are calculated with Fermi's golden rule from equation (S20) with free wave functions, including multiple m i incoming partial waves.The orange line is the Born approximation.The insert zooms in on small magnetic fields.
The Born approximation and Fermi's golden rule give the same results once all incoming m i partial waves are taken into account in the sum in equation ( S19).This is expected since they are both calculated in first-order perturbation theory.It is not obvious from the final expressions but can be seen in Fig. S3.At high magnetic field, only the m i = 0 channel is contributing since the rate scales as k f (k i /k f ) |mi| .However, at low enough magnetic field, the incoming and outgoing momenta are comparable and the other channels become as important, especially the m i = −2 → 0.

Pure-2D limit
We mentioned in the main text that the decay rate in the pure-2D case is indefinitely suppressed when the magnetic field increases.The reason is that as the field increases, so does the outgoing momentum, and the relaxation becomes shorter-ranged.The inward region is shielded by the dipolar potential which repels the atoms.We offer here a derivation of the shielded rate in pure-2D.We first derive the un-shielded case.Then we use a known [S21] zerotemperature shielded wave function and patch it at long-distance with a free low momentum wave function to compute the shielded rate.This will exhibit the spectacular behavior of a decreasing relaxation rate with an increasing magnetic field.We also show further suppression by going to lower temperature.We have not achieved this regime in our experimental setup, but we believe it is within reach with minor improvements, and think this theoretical treatment can provide insights into the system.
One starts with the equation ( S21) to write the 2D rate in the case of infinite transverse confinement ( χ0 (z) = δ(z)).The only decay channel is j f = 2 and

Free wave functions case
If we ignore the shielding, the wave functions take a simple form φm (ρ) = π k ρ L J m ( k ρ).In the low-temperature limit, the incoming Bessel function goes like J 0 ( ki ρ) 1.The integral

Shielded wave functions case
The result is radically different if we take into account the dipolar shielding potential.We will show two effects: the rate decreases both with increasing magnetic field and decreasing temperature.Lower temperature means greater shielding, and higher magnetic field shortens the range of the interaction which increases the shielding as well.To keep the model simple, we will only take the incoming wave function to be shielded since the outgoing one already experiences a centrifugal barrier.In pure-2D, the differential equation for the incoming wave function follows from equation (6): for which an analytical result is known from [S21] but only at zero temperature ( ki = 0) and involves modified Bessel functions.At finite temperature and in the absence of the dipolar interaction term, the solution is known with regular Bessel functions.We therefore use the modified Bessel function solution at short-range until ρ0 such that 2/ρ 3 0 = k2 i .We then patch this function with the finite temperature free solution with a phase shift δ.The incoming wave function then reads where K 0 is the modified Bessel function and J 0 and Y 0 the Bessel functions of the first and second kind.α is a normalization coefficient obtained by equating the wave functions and their derivative at the patching location ρ = ρ0 in a box of length L, which gives We now explore independently two limits: low temperature, and high magnetic field.

Low temperature limit
One can find the low temperature behavior of the incoming wave function.We expect the decay rate to be suppressed at low-temperature as the shielding increases.The r coefficient can be rewritten: , which can be used to expand δ through equation (S31) 4 : It gives the following expression for the incoming wave function at low temperature: This is reflected in the decay rate and β pure 2D shielded ∝ 1/ log ki 2 , which goes to zero at low temperature.This is very different from the free wave function case which had a finite limit even at zero temperature.The increase in the shielding factor with colder temperatures is shown in Fig. 3d.

Moderate magnetic field approximation
Here we derive an approximate analytical formula for the integral in equation (S28) with shielded wave functions.There is no simple expression for the integral of a Bessel K function multiplied by a Bessel J or Y and the dipolar potential5 .We have done full numerical calculations (see below).However, we can find an approximate analytical result for moderate magnetic fields by approximating the Bessel functions by ones having known integrals.Fig. S4 explains the different approximations we do to compute the integral.For high-enough magnetic field compared to the temperature6 , the outgoing Bessel wave function oscillates and has its first zero before the patching point ρ0 (see Fig. S4).We can then cut off the integral at this value as the remaining part is damped by the 1/ρ 3 potential.This is valid only if the incoming wave 4 Note that equating this equation with the form of the phase shift in 2D for a hardcore potential tan δ π 2(ln( ki ã/2)+γ) allows to recover the universal dipolar scattering result from [S21] that ã = 2e 2γ .function is not increasing too much after this zero, as it would compensate for the decrease due to the potential, which would result in the second lobe of the oscillation contributing more than the first one.The incoming wave function increases exponentially up to a certain distance ρi which sets an upper bound on the magnetic field our model tolerates7 .Defining D = x 2 / kf < ρ0 with x 2 5.13 the first zero of J 2 (x), the decay rate from equation (S28) is  The left part of the integrand J 2 ( kf ρ)/ρ 2 corresponding to the outgoing wave function, the potential and the surface element is plotted in solid blue.The outgoing wave function is cut off at the first zero in x 2 / kf and then Taylor expanded at this point inward.The right part of the integrand corresponding to the incoming wave function is plotted in solid red.It is patched in ρ0 to determine the normalization factor, and eventually expanded into a simple exponential function (dashed red).The total patched wave function (red and then green) is compared to the fully simulated wave function (dashed grey).Our model requires ρi < x 2 / kf < ρ0 .
The patching is only used to determine the normalization factor of φf .As we mentioned, we further simplify the functional form of the wave function as illustrated in Fig. S4.We use a Taylor expansion in ρ = x 2 / kf of the Bessel J function J 2 ( kf ρ) ρ2 gives overall  cos(δ)J 0 ( ki ρ0 ) − sin(δ)Y 0 ( ki ρ0 ) The agreement of this analytical expression with the numerical integration is presented in Fig. S5.

High field limit
The complexity of f makes difficult to grasp the behavior of the decay rate β.It is possible to do an expansion at high-magnetic fields of f which greatly simplifies the equation.However this expansion becomes valid only for fields of several thousands of Gauss, which invalidates the initial assumptions that x 2 / kf is larger than the limit ρi of the exponentially suppressed region of the incoming wave function.Nonetheless it describes well the overall behavior of the curve, and brings some insight into the suppression of the decay rate.In the large field limit, when setting the numerical and normalization prefactor κ to be:  or β pure-2D shielded ∝ B 1/8 exp −ξB 1/4 , which is a radically different behavior from the free wave function case where the rate is increasing linearly with the magnetic field.Fig. S5 shows this exponential suppression.The suppression factor can be made arbitrarily large.The overall shape of the curve is correctly reproduced by our analytical formulas (S32) and (S33) even at high fields, although our approximation does not fully capture the precise amplitude of nor the zeros in the decay rate that appear when x 2 / kf becomes < ρi .

Classical turning point
We have explained in the main text how the Franck-Condon principle predicts spin-flips to occur at the classical turning point of the outgoing wave function.The low temperature shielded situation we just presented provides a counterexample where the incoming wave function also has a classical turning point that needs to be taken into account.When the field is sufficiently high, the outgoing wave function oscillates multiple times in the suppressed region of the incoming wave function.Therefore the integrand gets contributions from multiple oscillations, not just the first one.

Fig. 1
Fig.1Principle of dipolar shielding (a.) Effective radial potential between two atoms from equation (6) for: no confinement (light blue, 3D), quasi-2D with ωz/2π = 300 kHz (steel blue) and pure-2D (dark blue).The incoming energy is given by the temperature T = 1 µK.(b.) Wave function solutions of equation (6) with n = 0 and initial orbital momentum m i = 0 for the three confinement strengths described above, and in red the spinflipped outgoing wave function to m f = 2 and B = 500 mG.The effect of shielding of the outgoing wave function is negligible for these parameters.(c.) Integrand of Fermi's golden rule (equation (2) and see equation(S21)  in Supplementary Information).Each curve is the product of the respective wave function in (b.), the outgoing wave function, and the double spin-flip operator from equation (4) integrated with the harmonic oscillator wave functions in the z-direction.The shielding we implement here corresponds to the difference between the light blue and steel blue curves.The minimum attainable decay rate for this incoming energy corresponds to the dark blue curve.See Supplementary Information for insights on the behavior of the integrand.

Fig. 3
Fig. 3 Theoretical loss rate coefficients.(a-b) Channel-by-channel decomposition of the dipolar relaxation rates of the m i = 0 incoming state (valid when µ B B k B T , see Supplementary Information) in a 300 kHz trap, both for free wave functions (a) and shielded ones (b).The blue and red colors correspond to single and double spin flips, respectively.The different shades correspond to different harmonic oscillator states as they open up with increasing magnetic field.(c-d)The suppression factor defined as the ratio of β 2D obtained from shielded and free wave functions, both at fixed temperature T 1 (c), and at fixed magnetic field B 1 (d).For each of the graphs we present curves for ω 1 /2π = 300 kHz (dotted), ω 2 /2π = 1.8 MHz (dashed dotted), the pure-2D case (solid line) as well as the analytical approximation (grey dotted) from equation (S32) detailed in Supplementary Information.

Fig
Fig. S1 Determining the zero of the magnetic field.For spin-polarized m J = −8 dysprosium atoms and left-circularly polarized imaging light, the drastic difference in Clebsch-Gordan coefficients for σ + and σ − transitions produces a step-like change in imaging signal as the magnetic field traverses through zero.

Fig
Fig. S2 Trap geometry and relevant length scales.Left: Reproduction of Fig. 2a of the main text.Right: The spatial density of the cloud in the longitudinal direction is characterized by the axial RMS width σz = az/2, the lattice spacing λ/2 and the initial width of the loaded thermal cloud σ ODT .

Fig. S3
Fig.S3Comparison between Fermi's golden rule and the Born approximation.The regions shaded in gray are calculated with Fermi's golden rule from equation (S20) with free wave functions, including multiple m i incoming partial waves.The orange line is the Born approximation.The insert zooms in on small magnetic fields.

Fig. S4
Fig. S4 Illustration of the approximations used to compute the shielded decayThe left part of the integrand J 2 ( kf ρ)/ρ 2 corresponding to the outgoing wave function, the potential and the surface element is plotted in solid blue.The outgoing wave function is cut off at the first zero in x 2 / kf and then Taylor expanded at this point inward.The right part of the integrand corresponding to the incoming wave function is plotted in solid red.It is patched in ρ0 to determine the normalization factor, and eventually expanded into a simple exponential function (dashed red).The total patched wave function (red and then green) is compared to the fully simulated wave function (dashed grey).Our model requires ρi < x 2 / kf < ρ0 .

2 
Fig.S5Dipolar decay rates in two dimensions.Shown are simulated and analytical decay rate coefficients for two magnetic field ranges.The blue curves represent the rate with simulated wave functions taking into account the dipolar repulsion for both the incoming and outgoing channels (dashed) or only on the incoming one (blue).The red curves present analytical results where the shielding is only accounted in the incoming wave function.The solid line is equation (S32) and the dashed one is its high-field limit, equation (S33).
Atoms in the initial spin state |j 0 =