Emergence of chaotic scattering in ultracold Er and Dy

We show that for ultracold magnetic lanthanide atoms chaotic scattering emerges due to a combination of anisotropic interaction potentials and Zeeman coupling under an external magnetic field. This scattering is studied in a collaborative experimental and theoretical effort for both dysprosium and erbium. We present extensive atom-loss measurements of their dense magnetic Feshbach resonance spectra, analyze their statistical properties, and compare to predictions from a random-matrix-theory inspired model. Furthermore, theoretical coupled-channels simulations of the anisotropic molecular Hamiltonian at zero magnetic field show that weakly-bound, near threshold diatomic levels form overlapping, uncoupled chaotic series that when combined are randomly distributed. The Zeeman interaction shifts and couples these levels, leading to a Feshbach spectrum of zero-energy bound states with nearest-neighbor spacings that changes from randomly to chaotically distributed for increasing magnetic field. Finally, we show that the extreme temperature sensitivity of a small, but sizeable fraction of the resonances in the Dy and Er atom-loss spectra is due to resonant non-zero partial-wave collisions. Our threshold analysis for these resonances indicates a large collision-energy dependence of the three-body recombination rate.


I. INTRODUCTION
Anisotropic interactions are a central and modern tool for engineering quantum few-and many-body processes [1]. A prominent example of such an interaction is the long-range dipole-dipole interaction (DDI) acting for instance between polar molecules [2], Rydberg atoms [3], or magnetic atoms [4]. Over the years, fascinating quantum effects of the anisotropy have been observed, such as the d-wave collapse of a dipolar Bose-Einstein condensate [5], the deformation of the Fermi surface [6], and the control of stereodynamics in dipolar collisions [7]. Moreover, the DDI is expected to give rise to a plethora of few-and many-body phenomena, which still await observation, such as universal few-body physics [8,9], rotonic features [10,11], two-dimensional stable solitons [12], and the supersolid phase [13].
Recently, atomic species in the lanthanide family became available to the field of ultracold quantum gases. The interaction between magnetic lanthanide atoms, such as Er [14,15] and Dy [16,17], is highly anisotropic. This is not only due to the long-range DDI, originating from their large magnetic moment, but also to the shorter-ranged van-der-Waals interaction [18], which exhibits anisotropic contributions arising from the large orbital angular momentum of their valence electrons.
For magnetic lanthanides, which also include the successfully laser-cooled elements Ho [19] and Tm [20], the orbital anisotropy is a consequence of a partially filled submerged 4f electron shell that underlies a closed outer 6s shell. This leads to an electronic ground state with a total atomic angular momentum  with j 1. Consequently, in collisions between such atoms there exist (j + 1) 2 non-degenerate (gerade) molecular Born-Oppenheimer (BO) potentials and a correspondingly large manifold of collision channels with associated molecular bound states. This is in sharp contrast to the one or two BO potentials encountered in alkalineearth and alkali-metal atom collisions. In addition, the anisotropy or orientation dependence of the BO potentials strongly mixes collision channels with large relative orbital angular momentum between the atoms even for our ultracold collisions with a = 0, s-wave initial channel. The complexity of the molecular forces are reflected in a dense spectrum of Fano-Feshbach resonances as a function of magnetic field B, as recently observed in Er [14,21] and Dy [22]. In Er a statistical analysis of the spacings between resonances has shown correlations that revealed chaotic scattering. The data set of the initial Dy experiments was too small to extract statisticallysignificant correlations.
Chaotic behavior is manifest in a variety of complex systems ranging from atomic to nuclear and solid-state physics. In atomic physics chaos was originally studied with Rydberg states of H and He in a magnetic field [23]. Later on a variety of more complex atoms and ions in highly-excited states have shown signatures of chaotic spectral distributions [24]. The origin of chaos in these systems was traced back to a strong mixing of many-electron excited states by the Coulomb interaction [25]. A chaotic level distribution is also common in a variety of solid state systems ranging from those with strong many-body interactions to the motion of particles in irregular potentials [26,27]. Experiments in nuclear physics [28,29] have also produced substantial evidence for chaotic neutron resonance spectrum fluctuations, which agree with predictions of random matrix theory (RMT). Similar agreement was found from numerical simulations based on nuclear shell models [30,31]. Moreover, Refs. [32,33] suggested that chaos is a generic property of nuclei with multiple degrees of freedom (i.e. multiple active shells), which become completely mixed.
This article describes a joint effort to understand ultracold scattering and Fano-Feshbach spectra of strongly magnetic Er and Dy atoms. In particular, we report on the measurement and statistical analysis of Fano-Feshbach spectra for Dy and Er between B = 0 G to 70 G at gas temperatures T below and around 1 µK. Here 1 G=0.1 mT. We observe that both elements have similar chaotic scattering. We present a random-matrix-theory (RMT) inspired model to gain insight into their statistical properties as well as theoretical evidence based on coupled-channels calculations with a microscopic Hamiltonian that chaotic scattering requires both strong molecular anisotropy and Zeeman mixing to fully develop. Limitations of the RMT are also discussed. Finally, we present experimental data and a comparison to a resonant trimer model to show that our increase in resonance density with temperature is a consequence of the strong collision-energy dependence of transitions from entrance d-wave channels of three free atoms to resonant trimer states.

A. Measurement
The experimental study of Fano-Feshbach resonances in Er and Dy is based on high-resolution trap-loss spectroscopy on spin-polarized thermal samples. Ultracold bosonic 164 Dy samples are created by direct loading from a narrow-line magneto-optical trap (MOT), operating on the 626 nm cycling transition, into a single-beam optical dipole trap (ODT) [34]. By moving the last focusing lens of the ODT, the atoms are transported from the MOT chamber to the science cell. The ODT is created with a 100 W fiber laser at a wavelength of 1070 nm. We achieve a transport efficiency close to unity. This fiber laser, however, causes atom loss due to its longitudinal multimode structure [35]. Therefore, we transfer the atoms into a second single-beam ODT, created by a 55 W solid-state laser at a wavelength of 1064 nm. Finally, forced evaporative cooling in a crossed ODT leads to a sample of 10 5 atoms in the energetically-lowest Zeeman sublevel m j = −8 at T = 600 nK.
High-resolution trap-loss spectroscopy is performed on a spin-polarized bosonic 168 Er sample at T = 1400 nK and compares this spectrum with that obtained at a four times lower temperature measured in previous work both for 168 Er as well as for fermionic 167 Er [21]. The experimental procedures for creating bosonic and fermionic samples are described in Ref. [14] and [15], respectively. Bosons (fermions) are prepared in the lowest Zeeman sublevel m j = −6 (m f = −19/2). Erbium samples are trapped in a crossed ODT and contain about 10 5 atoms.

B. Feshbach spectroscopy
Feshbach spectroscopy is performed in a similar manner for the two species. The magnetic field is ramped up over a few milliseconds to a magnetic field value B, where the atoms are held in the ODT for 500 ms for Dy, 400 ms for 168 Er, and 100 ms for 167 Er. During this time inelastic three-body recombination causes atom loss from the ODT. At resonance, the recombination process is enhanced because of the coupling between the atomicthreshold state and a molecular state leading to a resonant increase of the atom loss. We identify the field locations of maximum loss as the positions of Fano-Feshbach resonances [36]. The atom number is probed by standard time-of-flight absorption imaging at low magnetic field. We record atom-loss features for magnetic field values between 0 G and 70 G in steps of a few mG. Figure 1( The understanding of the richness of the scattering in Er and Dy requires the development of sophisticated microscopic coupled-channels scattering models. We defer such analysis until later in this paper and first analyze our data following the statistical approach based on the random-matrix theory (RMT) advocated by Ref. [21]. In particular, we study the correlations between resonance locations via the nearest-neighbor spacing (NNS) distribution and setup a RMT-like model, which accounts for the structure of our B-dependent microscopic Hamiltonian to get intuition about these NNS. In our description of the coupled-channels calculations limitations of such a RMT-like model are discussed.  0  2  4  6  8  10  12  14  16  18  20  22  24  26  28  30  32  34   10   -1   10   0   36  38  40  42  44  46  48  50  52  54  56  58  60  62  64  66  68

C. Statistical analyses
Our statistical analysis starts with the construction of the staircase function, which is a step-like function that counts the number of resonances below magnetic-field value B [37]. Figure 1(b) shows the staircase function for Dy and Er. For both species the function is well-fit by a linear curve forced to pass through the origin. Its slopeρ corresponds to the density of resonances. Deviations below and above the fit occur for small and large B, respectively. The fitted resonance densities are given in the caption. Remarkably, the density of resonances of 168 Er at T = 1.4 µK is 25 % higher than the one observed at 350 nK. The discussion of the origin of this sensitivity is postponed until Section V. The density ρ for bosonic Dy is 50% larger than for bosonic Er. This is caused by the larger  of Dy and, thus, its larger number of allowed collision channels. The much larger density ρ of 25.6 G −1 for the fermionic 167 Er is due to its additional hyperfine structure.
Fluctuations in the number of resonances within a magnetic field interval ∆B is a second measure of the statistical properties of the spacings between resonances. Formally, it is defined as the dimensionless number vari-  Figure 1(c) compares Σ 2 for our Dy and Er data as a function of N . The fluctuations for both species monotonically increase with ∆B but are substantially less than the shot noise limit. While this behavior was previously demonstrated for Er [21], the present results provide the first evidence of correlation in Dy and indicate similarity between the species.
These correlations between resonance locations are further studied using the nearest-neighbor spacings distribution P (s) where s = δB ρ and δB is the field spacing between two adjacent resonances in the spectrum.
, an empirical function that interpolates between P P (s) and P WD (s) for η = 0 and 1, respectively, and b is a normalization constant [38]. The values for η reported in the caption, indicate intermediate or mixed behavior of the data. We present the magnetic-field resolved Brody parameter η(B) in Figs. 2(c) and (d) obtained from a fit to the NNS distribution of resonances located in moving intervals [B − ∆B/2, B + ∆B/2] with ∆B = 20 G. It has a non-negligible 1σ uncertainty equally limited by the quality of the fit and the number of Feshbach resonances in an interval or bin. The latter uncertainty is reflected in the bin-to-bin variation of η(B). For Dy we observe that η increase linearly with field for small B, which saturates at a value of ≈ 0.5 for B > 30 G. For Er the Brody parameter fluctuates around 0.5. Interestingly, the Er data at our two temperatures have a similar behavior, indicating that the larger density of resonances at higher T does not impact the degree of correlation between their spacings.

III. RMT ENSEMBLE MODEL
Random matrix theory is based on the powerful notion that the statistics of eigenvalues and eigenfunctions of a complex system can be studied by replacing the microscopic Hamiltonian by an ensemble of random Hamiltonians. In this spirit, we construct a RMT-inspired model for weakly-bound molecular dimer states to test the distribution of Fano-Feshbach resonances.
Our RMT model is based on the statistics of eigenvalues of the N × N real, symmetric matrix H RMT = H 0 + H Z , where matrices H 0 and H Z represent the B = 0 Hamiltonian and the Zeeman interaction of the two atoms, respectively. Without loss of generality we can assume that H Z is a diagonal matrix with matrix elements given by mgµ B B, where m is an integer between −2j and 2j, corresponding to the sum of the projection quantum numbers of the atomic angular momenta, g is the atomic Landé factor, and µ B is the Bohr magneton. The Zeeman interaction does not depend on the rotational state of the molecule and, thus, entries in H Z correspond to states with a definite value for and its projection. H 0 is then the B = 0 Hamiltonian expressed in this basis. It is also convenient to define H 0 = H d + H cpl , where diagonal matrix H d contains the diagonal matrix elements of H 0 and H cpl is the matrix of all its off-diagonal elements. The eigenvalues of H d can then be interpreted as the energies of ro-vibrational levels of the isotropic contribution of the molecular BO potentials, while H cpl describes mixing due to the anisotropic contributions of these potentials.
We generate members of our ensemble of H RMT by   [39].
We apply the RMT model to the case of 168 Er. The relevant species-specific quantities are j = 6, g = 1.16, and d is chosen to roughly reproduce the observed density of Fano-Feshbach resonances of 168 Er and is set to shows an example of a molecular spectrum, the eigenvalues of H RMT obtained with our RMT model as a function of B with η d = 0 and ν cpl /h = 2 MHz. We observe that as B increases weakly-bound molecular states avoid each other multiple times before reaching the twoatom threshold creating a Feshbach resonance. When we turn off H cpl the levels cross. Similar B-field dependencies of the eigenvalues occur for η d > 0.
We investigate the effect of the parameters ν cpl and η d on the NNS distribution of the Fano-Feshbach resonances as well as that of the B = 0 molecular levels. Figure  3(b) shows the NNS distribution of Feshbach resonances, obtained by averaging over 15 realizations of H RMT , for four values of ν cpl and η d = 0. For negligible ν cpl the distribution follows P P (s) and approaches P WD (s) when the anisotropic coupling strength ν cpl is large compared to d . In fact, we find that a larger d requires a larger ν cpl to develop correlations.

Figures 3(c) and (d) show
Brody parameters fit to NNS distributions as functions of ν cpl and η d . Panel (c) shows η for the B = 0 molecular binding energies. For ν cpl = 0 the Brody parameter is simply η d as expected from the distribution of the diagonal H d , while for larger interaction anisotropy ν cpl the parameter η ≈ 0.9, close to a Wigner-Dyson distribution, independent of η d . Panel (d) shows η extracted from the RMT Feshbach resonance locations as a function of ν cpl . It suggests that the correlation in the NNS of the resonances is caused by ν cpl , whereas it appears fairly independent on η d . More precisely, the Brody parameter fitted to these distributions rapidly increases from η ≈ 0 to η ≈ 0.8 for ν cpl < ∼ d and tends to one for larger ν cpl . We conclude from the RMT model that the correlations between the locations of the Fano-Feshbach resonance are essentially due to the avoided crossings between weakly-bound molecular states at finite B and only weakly dependent on the energy distribution at B = 0. In fact, these correlations increase for increasing ν cpl .

A. Realistic setup
A quantitative understanding of the origin of the chaotic resonance distribution requires coupled-channels and bound-state calculations with physically realistic angular-momentum couplings and interaction potentials. We do so here based on the time-reversal symmetric Hamiltonian for the relative motion of Dy and Er described in Refs. [21,40]. It contains the Zeeman Hamiltonian, the molecular vibration and rotation, and the molecular interactions with isotropic (orientationindependent) and anisotropic (orientation-dependent) contributions,V i (R) andV a ( R), respectively, where R describes the separation R and orientation of the atom pairR. The potential has eight tensor operators coupling the two atomic and relative orbital angular momenta,  1 ,  2 , and . For B = 0 the total angular momentum J =  1 +  2 + is conserved. For B > 0 G only the projection M of J along B is conserved. The zero-of-energy of the Hamiltonian is the energy of an atom pair in the absolute lowest Zeeman sublevel m jα = −j α .
The potentialsV i (R) andV a ( R) contain short-ranged exchange, medium-ranged van der Waals as well as longrange magnetic dipole-dipole interactions. We use the isotropic van der Waals coefficient C 6 = 1723 E h a 6 0 and anisotropic coefficients spread over ∆C 6 = 174 E h a 6 0 for Er [21]. For Dy we have improved the value of van der Waals coefficients of Ref. [40] by including additional experimental and theoretical transition frequencies and oscillator strengths [41][42][43][44] and now use C 6 = 2003 E h a 6 0 and spread ∆C 6 = 188 E h a 6 0 . In particular, the anisotropic spread for Dy has significantly increased.
Here, E h = 4.360 × 10 −18 J is the Hartree energy and a 0 = 0.05297 nm is the Bohr radius.

B. Bound-state calculations
In Ref. [21] we performed initial coupled-channels calculations of the scattering between ultracold Er atoms and predicted that tens of partial waves should have been included as the strength of the anisotropic contribution is large. We, however, were unable to reach numerical convergence with respect to the number of coupled equations.
Here, we circumvent this limitation by performing multichannel bound-state calculations, in which we use B = 0 eigenstates as a basis for those at B > 0 G. For B = 0, where J is a good quantum number, at most 49 and 81 Bose-symmetrized and parity-conserving channels are coupled for Er and Dy, respectively. The B = 0 coupled Schrödinger equations are discretized on the interval R ∈ [0, R max ] assuming zero boundary conditions and solved as a matrix eigenvalue problem [45][46][47][48]. For each J only eigenstates with energies between [E 0 , E 1 ] surrounding the zero of energy are computed and stored. The bound states for B > 0 G are solutions of the matrix eigenvalue problem that includes all computed B = 0 solutions with |M | ≤ J ≤ J max and their coupling due to the Zeeman interaction. Selection rules of the Zeeman interaction ensure that there only exists direct coupling between J and J zero-field eigenstates with J −J = 0, ±1. For both species R max = 1000 a 0 , E 0 /h = −3 GHz and E 1 /h = 0.9 GHz ensuring that Feshbach resonance locations below 70 G are converged.
In this section on the microscopic calculations we focus on analysing the spectra at our coldest temperatures, where the initial collision channel has s-wave ( = 0) character. Hence, we only need to consider even-channels with total projection quantum number M = −12 and −16 for 168 Er and 164 Dy, respectively and inclusion of zero-field solutions up to J max = 36 for Dy and 39 for Er is sufficient to reproduce the experimental resonance densities. In Sec. V we will discuss higher temperature collisions between Er atoms, where d-wave ( = 2) entrance channels must be considered and, hence, spectra at other M values (i.e. M between −14 and −10 for 168 Er) contribute.

C. Interaction anisotropies
We first look into the role of interaction anisotropies on the level distribution of the most-weakly bound molecular energy levels at zero magnetic field. There are two dominant components to the anisotropy, the dispersion V ∆C6 ( R) and magnetic dipole-dipole V MDD ( R) contribution. To distinguish the contributions of these two terms, we define: with variable strength λ ∆C6 and λ MDD . We systematically increase the strengths λ ∆C6 and λ MDD from zero, where we recover the full physical strength for λ MDD = λ ∆C6 = 1.
For completeness we note that the dominant tensor operator for the anisotropic dispersion contribution is with strength c a < 0 found with the methodology discussed in the subsection IV A. Weaker contributions indicated by dots are included in our calculations. Moreover, where µ 0 is the magnetic constant. and dispersive (λ MDD = 0, varying λ ∆C6 ) anisotropic interaction, respectively. For λ MDD = λ ∆C6 = 0 the binding energies are regularly structured with many near degeneracies. In fact, the corresponding states are rovibrational levels of the isotropic centrifugal potentialŝ V i (R) and labeled by . In our 3 GHz energy window an swave channel has at most three bound states, while even > 0 channels with their centrifugal barriers have fewer [21,49]. For small λ ∆C6 and λ MDD the degeneracy is lifted and levels shift linearly. The linear dependence for increasing strength of the dipole-dipole is approximately valid up to the physical value of λ MDD = 1. Hence, the dipole-dipole interaction does not lead to our chaotic level distributions. In fact, Fig. 4(a) shows that at λ MDD = 1 and λ ∆C6 = 0 the NND distribution is Poissonian.
On the other hand, for a relatively small anisotropic dispersion strength λ ∆C6 ≈ 0.1 levels start to avoid each other. Starting from λ ∆C6 ≈ 0.5 most avoided crossings are noticeable on the 3 GHz scale of the figure. At the nominal λ ∆C6 = 1, where there are 56 levels with −3 GHz < E/h < 0 GHz, a significant fraction of the levels have undergone multiple avoided crossings and can not be described by a single dominant partial wave. The level spacing is chaotic as confirmed by the NND distribution for λ ∆C6 = 1 and λ MDD = 0 in Fig. 4(c). We have computed the weakly-bound J = 16 levels for λ MDD = λ ∆C6 = 1. Visually the level distribution is much the same as the one shown in Fig. 4 Fig. 4(b) and (d). For increasing dipole-dipole strength λ MDD and no anisotropic dispersion (λ ∆C6 = 0) the parameter is always zero indicating the prevalence of small level spacings. On the other hand, in the absence of the dipole-dipole interaction, increasing λ ∆C6 leads to an increasing η. It evolves from η = 0.2 for λ ∆C6 < ∼ 0.5 to η = 0.7 for λ ∆C6 = 1, indicating a depopulation of small energy spacings. Note that our systems does not reach a Wigner-Dyson distribution, which corresponds to η = 1.
We compare in Fig. 4(f) two NNS distributions of B = 0 weakly-bound states of 164 Dy 2 obtained for the full anisotropic interaction (λ MDD = λ ∆C6 = 1). Both distri- The individual-J NNS distribution is non-Poissonian as levels with the same J repel each other. The combined-J distribution, however, follows a Poisson distribution indicating that energies of bound states with different J are uncorrelated. In other words, even though the Hamiltonian, i.e. the set of coupling operators between  1 ,  2 and , is the same, differences in the matrix elements and thus coupling strengths between channels lead to uncorrelated eigen energies.

D. Atom scattering in a magnetic field
The study of the B = 0 G multi-channel bound states has shown that interaction anisotropies mix channels with the same J, while states with different J are uncorrelated. The Zeeman interaction mixes these molecular levels and leads to the Fano-Feshbach spectrum. cluded, respectively. The figure shows that the Dy level density is higher than that for Er. This simply follows from the larger atomic angular momentum of Dy, leading to a larger number of channels with the same J − |M |. We also observe that for both species the level structure in the 0 G to 10 G, small field region is qualitatively different from that in the larger field region. For small B the avoided crossings are substantially narrower than for larger B. Moreover, at small field the levels cluster, while at larger field they are more uniformly distributed. These changes are a consequence of the linearly-increasing Zeeman coupling between vibrational levels with different Js as a function of B. Figure 6(a) shows effective length a s (B) as a function of B. It diverges at every resonance location and is closely related to the scattering length of a zero-energy collision. Our calculations can not be directly used to define the scattering length as we use a hard-wall potential for R ≥ R max . This wall leads to a discrete set of states with positive energy and using the lowest of these E s (B) we can define the effective length a s (B) shown in the figure by solving for E s (B) =h 2 π 2 /[2µ r (R max − a s (B)) 2 ] with µ r = m/2 and atomic mass m [50].
It is of interest to briefly discuss the convergence properties of our calculations. The data in Figs. 6(a) and 5(c) and (d) are based on computations with channels with J up to J max = 39. Figure 6(b) shows the 168 Er 2 Feshbach resonance densityρ as a function of J max . The resonance density increase linearly from ≈ 0.5 1/G at J max = 12 but then is seen to "saturate" for larger J max . At J max = 39 the experimental density is reproduced. In addition, Fig. 6(c) shows the field location of resonances between 50 G and 55 G as a function of J max . The res- onance locations change significantly for J max < 22, but then rapidly converge. This implies strong mixing among bound states with those J. On the other hand, the location of resonances that appear for J ≥ 22 are almost immediately converged indicating weak mixing to smaller J states.

E. Comparison of experiment and coupled-channels model
In Figs. 2(a) and (b) we show the NNS distribution of converged Feshbach-resonance locations based on our multi-channel data between B = 0 G and 70 G for 164 Dy 2 with J max = 36 and 168 Er 2 with J max = 39, respectively. For both species the distribution clearly deviates from a Poisson distribution, consistent with the experimental distributions that are also shown. The fitted experimental and coupled-channel Brody parameters agree within their error bars.
The anisotropy parameters λ ∆C6 and λ MDD in the coupled-channels calculations and the parameter ν cpl in the RMT play analogous roles in the Hamiltonian and in the emergence of chaotic level distributions, even though no explicit quantitative connection exists. This role is most manifest in the Brody parameters of the B = 0 G bound states and that of the Feshbach resonance spectra for the two models. For 168 Er the corresponding Brody parameters from the coupled-channels calculations are ≈ 0.01 and 0.68 at the physical λ MDD = λ ∆C6 = 1, respectively. Within the RMT model the small η value for the B = 0 G level distribution requires weak coupling ν cpl d and η d ≈ 0. In contrast the Brody parameter for the Feshbach resonance spectrum requires ν cpl ≈ d and points at limitations of the current RMT model. Similar conclusions hold for bosonic Dy. Future advanced RMT models might circumvent these limitations by incorporating overlapping, uncoupled chaotic series as is found from our B = 0 G coupled-channels calculations.
We plot the B-field-resolved Brody parameter η(B) of the theoretical coupled-channels data in Figs. 2(c) and (d). A comparison with the experimental η(B) shows excellent agreement for 164 Dy, while the agreement for 168 Er is less satisfactorily. A possible explanation for the discrepancies in 168 Er is the larger bin-to-bin fluctuations as bins contain fewer resonances than for 164 Dy.
For 164 Dy the theoretical field-resolved Brody parameter in Fig. 2(c) linearly increases from zero for small B fields and saturates at η(B) ≈ 0.5 for fields larger than 35 G where the size or width of the avoided crossings between weakly-bound states is larger. For 168 Er in Fig. 2(d) we find a much more rapid increase of η(B) at small fields. This is followed by a plateau at η(B) ≈ 0.5 between B = 20 G and 50 G, after which η(B) → 0.9 with an uncertainty of 0.2 close to a Wigner distribution. The initial rise of η(B) for both atomic species is a consequence of weakly-bound vibrational levels, uncoupled and randomly distributed when B = 0 G, that start to repel each other as the Zeeman interaction increases in strength for increasing B. The plateau at η(B) ≈ 0.5 and the sudden increase of η(B) to one for 168 Er have no simple explanation and are determined by the notfully-explored complex interplay between the Zeeman and anisotropic inter-atomic interactions. It does, however, indicate that Wigner's assumptions on ensembles of Hamiltonians do not hold for fields below 50 G.

V. TEMPERATURE DEPENDENCE OF THE RESONANCE DENSITY
We now describe the origin of the strong temperature dependence of some of the resonances in our atom-loss spectra and thus explain the accompanying increase of the resonance density. Here, atom loss is solely due to three-body recombination, where three ultra-cold atoms collide to form a di-atomic molecule and an atom that both are lost from the atom trap. Figure 7(a) shows atom-loss spectra for one such resonance for 168 Er at four temperatures below 2 µK. Atom loss, indeed, is larger for larger temperatures, but we also observe a broadening of the B-field width and a shift of the maximum loss position to larger B fields. Resonances with a weak temperature dependence show none of these behaviors.
We show with an intuitive resonant "trimer" model that a strongly temperature-dependent resonance is due to scattering processes with entrance d-wave channels even though the two-body d-wave centrifugal barrier, V b /k B = 250 µK, is a hundred times larger than our highest temperature, where k B is the Boltzmann constant. The difference in the power-law Wigner-threshold behavior of the recombination rate with collision energy for s-and d-wave entrance-channel collisions can explain our observations. Three-body recombination has been extensively studied in the context of Efimov physics [51][52][53][54]. We follow Refs. [55][56][57] and start from a coupled-channels description in the (mass scaled) hyperradius ρ, which describes the size of the three-atomic system, and basis functions in the five other hyperspherical coordinates that are ρdependent eigen states of the squared "grand-angular momentum operator". Similar to the coupled-channels description for two atoms, there are entrance, open, and closed channels. The collision starts in one of the entrance channels with atoms in the energetically-lowest Zeeman state and relative three-body kinetic energy E 3 , the dimer plus atom are the open channels, and bound states in closed channels can lead to resonances. These closed channels dissociate to three free-atom states with at least one atom in a Zeeman level with higher internal energy. The bound states are resonant "trimer" states giving us our name for the model. It should, however, be realized that their origin lies in bound states of pairs of atoms and that the resonant state is better thought of as a pair bound state that hops from pair to pair. We define E 3 ≡h 3 k 2 3 /(2µ 3 ) ≡ µ 3 v 2 3 /2 with the three-body reduced mass µ 3 = m/ √ 3, where k 3 and v 3 are the relative wave vector and velocity, respectively.
The potentials in the entrance channels have longrange repulsive centrifugal potentials, governed by the asymptotic behavior of the grand-angular momentum operator and depend on the relative orbital angular momentum N of the three atoms. In fact, the centrifugal potentials areh 2 (λ + 3/2)(λ + 5/2)/(2µ 3 ρ 2 ) with nonnegative integer quantum number λ [55]. For N = 0 the least repulsive potential has λ = 0, while that for N = 2 has λ = 2.
For an isolated trimer resonance in a closed channel coupled to both entrance and other open channels we can apply the resonance theories by Fano and Feshbach and derive that the recombination rate coefficient at collision energy E 3 and entrance channel with quantum number λ is given by is a resonant expression for the square of a dimensionless S-matrix element, where B 0 is the trimer resonance location and µ is the magnetic moment of the resonant trimer relative to that of the entrance channel. The definition for |S(E 3 , B)| 2 also contains the entrance-channel energy width Γ(E 3 ) = A λ E λ+2 3 to the trimer resonance with a characteristic power-law energy dependence that reflects the threshold behavior of the scattering solutions in the centrifugal potentials. The energy width Γ br determines the decay or breakup rate of the resonance into the fast atom and dimer pair and is independent of E 3 . Finally, Γ tot (E 3 ) = Γ(E 3 ) + Γ br . For simplicity we have assumed that non-resonant, direct recombination from the entrance to open channels is weak. We also note that for N = 0 and λ = 0, L 3 (E 3 , B) approaches a finite constant for E 3 → 0 as expected.
In our experiments we have thermal samples of Er and we require the thermally-averaged rate coefficient and normalization Z = ∞ 0 E 2 dE e −E/kT = 2(kT ) 3 . In order to increase signal to noise we have allowed a significant fraction of atoms to be lost (See Fig. 7(c)), which assuming a homogeneous sample can be modeled by the rate equation dn(t)/dt = −3L 3 (T, B)n 3 (t) for atom density n(t) [53] with solution where N (t h , B) is the remaining atom number after hold time t h , N 0 is the initial atom number, and n 0 is the initial density. This non-linear time evolution adds additional broadening to the lines. Figures 7(b) and (c) show our model event rates L 3 (T, B) as a function of B for N = 0, λ = 0 and N = 2, λ = 2 , respectively. Curves are for the same four temperatures as in Fig. 7(a). A comparison of panels (b) and (c) shows a striking difference. The strongest features in panel (b) are for the smallest temperatures, while those in panel (c) are for the largest temperatures. This behavior naturally follows from an approximation of the integrant in L 3 (T, B) under the conditions kT Γ br Γ(E) [58]. In this limit the Lorentzian is sharply peaked around E 3 = µ(B − B 0 ) for B > B 0 and after some algebra it follows that L 3 (T, B) as a function of B has a maximum value proportional to (kT ) λ−1 located at B = B 0 + (λ + 2)kT /µ. Consequently, for λ = 0 and 2 the maximum loss rate decreases and increases with T , respectively. Even for less restrictive parameter values as used in Fig. 7 this trend remains.
Our experimental data has a temperature trend as in panel (c). In fact, Fig. 7(a) compares our experimental loss data with model N (t h , B) for N = 2, λ = 2 using the same parameters as in Fig. 7(c) and requiring a ≈ 50% maximum atom loss as in the experiment. It is worth noting that, from our theoretical calculations, the magnetic field width of L 3 (T, B) is noticeably smaller than that for N (t h , B) indicating that the finite hold time does indeed lead to broadening. The agreement of the experimental data and the prediction of our model for the losses is satisfactory for all four temperatures given the limitations and approximations within our modeling. We conclude that our strongly T -dependent resonances correspond to d-wave or more precisely N = 2 entrance channel collisions. Note that we haven't observed any resonances with temperature dependence similar to Fig. 7(b) in our spectra. In the case of resonances with three-body swave entrance channel, which would correspond to such a dependence, we infer that the loss spectra is saturated. This will be subject for future investigations.
As a corollary, this implies that for two colliding atoms as described in Sec. IV temperature-dependent resonances are due to collisions with entrance d-waves for which there are multiple allowed values of the total angular projection quantum number M . Here M = -14 to -10 for bosonic Er and M = -18 to -14 for bosonic Dy. Numerical computations, not presented here, show that their zero-field bound-states and thus resonance locations are again uncorrelated and random.

VI. CONCLUSION
In summary, we have experimentally and theoretically studied the resonant scattering of ultracold Er and Dy atoms in a magnetic field. We have shown that chaotic scattering as witnessed by chaotic nearest-neighbor spacings between Feshbach resonance locations emerges due to the anisotropy in the molecular dispersion.
Our study has also revealed several unique features of colliding magnetic lanthanides that have not been observed in any other ultracold atomic system. These lanthanides are characterized by their exceptionally large electron orbital angular momentum, which lead to large anisotropic dispersion interactions between these atoms. Our theoretical estimate shows that in both Er and Dy collisions the ratio of anisotropic to isotropic dispersion interaction ∆C 6 /C 6 is about 10%. This anisotropy leads to significant splittings among the 48 and 81 gerade short-range potentials that dissociate to the ground-state atomic limits of Er and Dy, respectively. We have shown that each potential has its own ro-vibrational structure, which by coriolis forces and the Zeeman interaction interacts with that of other potentials, creating a dense distribution of levels near the threshold and initiating chaos. In fact, we find a very large number of partial waves contributing to the creation of Fano-Feshbach resonances.
On the other hand, if we just consider the anisotropy from the magnetic dipole-dipole interaction alone, our coupled-channel calculations indicate that chaos in the level distribution does not appear. The strength of the dipole-dipole interaction is too small. In addition, we have shown that the NNS distributions for Dy and Er are very similar, as can be expected from their similar ∆C 6 /C 6 ratio. The difference in their magnetic moment only plays a small role. This further confirms that chaos is due to the anisotropic dispersion interaction.
The distribution of Feshbach resonances of ultra-cold ground-state alkali-metal, alkaline-earth, Yb, and Cr atoms, as experimental studies have shown, are not chaotic. This is because these atoms have a zero electron orbital angular momentum and, hence, only an isotropic dispersion interaction. Eventhough alkali-metal and Cr atoms have a non-zero magnetic moment of 1µ B and 6µ B , respectively, these moments do not lead to chaos. We would expect that other magnetic lanthanides and actinides with non-zero orbital angular momentum will exhibit chaotic properties in their collisions. In addition, collisions between mixed species, such as magnetic lanthanides and alkali-metals, like K+Dy or Na+Er, might be susceptible to chaos. A first theoretical analysis for Li+Er [59], however, estimates a small 2% dispersion anisotropy and no chaos is predicted.
Another interesting property of magnetic lanthanide gases is the extreme sensitivity of the atom-loss spectra and, in essence, three-body recombination to the temperature. This phenomenon was first observed in Ref. [22] for loss spectra of Dy. The number of Dy resonances increases by 50% when the temperature is increased from 420 nK to 800 nK. Here, we observe a 25% increase in the Er resonance density when the temperature rises from 250 nK to 1400 nK. We have shown by a comparison of resonance profiles taken at several temperatures and predictions of a theoretical model of three-body recombination via the formation of a trimer, or more precisely of a shared pair bound state, that the origin of the temperature-dependent resonances lies in the "partial wave" of the three-atom entrance channel. Entrance channels with zero and non-zero total orbital angular momentum N lead to line shapes with a different temperature behavior. Those with N = 0 or "s-wave" entrance channels have sharply decreasing recombination rates with temperature, whereas those with N = 2 or "d-wave" entrance channels have an increasing recombination rate. Temperature sensitive resonances can only be explained by "d-wave" collisions. It is worth noting that for alkali-metal-atom collisions a number of entrance-channel p-wave resonances have been observed (See for example Ref. [60] for cesium). Analysis of the temperature dependent rate coefficient, however, was not performed.