## Abstract

Ergodic quantum systems are often quite alike, whereas nonergodic, fractal systems are unique and display characteristic properties. We explore one of these fractal systems, weakly bound dysprosium lanthanide molecules, in an external magnetic field. As recently shown, colliding ultracold magnetic dysprosium atoms display a soft chaotic behavior with a small degree of disorder. We broaden this classification by investigating the generalized inverse participation ratio and fractal dimensions for large sets of molecular wave functions. Our exact close-coupling simulations reveal a dynamic phase transition from partially localized states to totally delocalized states and universality in its distribution by increasing the magnetic field strength to only a hundred Gauss (or 10 mT). Finally, we prove the existence of nonergodic delocalized phase in the system and explain the violation of ergodicity by strong coupling between near-threshold molecular states and the nearby continuum.

## INTRODUCTION

Ultracold atomic physics is now poised to enter a new regime, where far more complex atomic species can be cooled and studied. Magnetic lanthanide atoms with their large magnetic moment and large orbital angular momentum are extreme examples of these species. Ultracold gases of magnetic lanthanides provide the opportunity to examine strongly correlated matter, creating a platform to explore long-range dipolar magnetic interactions and exotic many-body phases, such as quantum ferrofluids, quantum liquid crystals, and supersolids. In particular, experimental advances in trapping and cooling Dy, Ho, Er, and Tm atoms (*1*–*8*) are paving the way toward these goals.

Most lanthanides have a partially filled submerged 4f-electron shell, which lies beneath a filled 6s^{2} shell. The 4f electrons are added with parallel spins, leading to both a large magnetic moment and a large total and orbital angular momentum. These “unquenched” angular momenta generate a rich atomic and molecular structure as well as collective phenomena of these systems (*8*, *9*).

In previous studies (*8*–*11*), we explored the scattering and interactions between neutral magnetic lanthanide atoms in the ultracold regime with temperatures well below 1 mK. Over the past few years, we developed a framework for understanding the complex anisotropic interactions between these dipolar atoms and their chaotic behavior. In particular, our theoretical model, in combination with advanced numerical treatments, bridges the conceptual gap between “simple” alkali-metal and alkaline-earth atoms and complex lanthanides. It allows us to explain the origin of the dense spectra and chaotic character found in the statistics of the observed Er and Dy collisional resonances as due to the anisotropy of both the short- and long-range interactions between the atoms.

Here, we apply modern theoretical approaches from quantum theories of chaos to characterize the observed chaotic distribution of Fano-Feshbach resonances in weak magnetic fields. We turn to an analysis of near-threshold bound states and, in particular, the eigenstate wave functions of our molecular Hamiltonian matrices, which include both localized and continuum basis states or configurations. By introducing generalized moments of wave function mixing coefficients, we are able to identify extended eigenfunctions as well as a number of localized eigenstates in our system. Localized states are characterized by a wave function with a small percentage of mixing. Our extended eigenfunctions exhibit strong, disordered mixing that is partially chaotic and approaches a unique type of universal behavior. The degree of mixing is controlled by the external magnetic field *B*. We compute the fractal statistics of the extended states and attempt to answer the question of whether our quantum system has ergodic properties where all accessible configurations are equally available.

Our model goes beyond the canonical theory of spectral statistics. It focuses on a discrete set of bound rovibrational levels coupled to scattering continua, where states exist for any energy. The level statistics are only chaotic in the energy region of the bound states. In practice, we treat continua with a discretization procedure by introducing a large spatial/radial cutoff. As a result, we include a large number of continuum states to converge the eigenenergies and eigenvectors of the selected bound states. This approximation works well for many scattering problems.

Other types of ultracold atomic and molecular physics experiments are also starting to influence the developments in the quantum chaos and disorder theories. This is mainly due to the fact that these experiments have an unprecedented control and tunability of system parameters. Disorder potentials, for example, can be used to study Anderson localization (*12*). For ultracold atoms and molecules, these potentials can be created with speckle patterns, obtained by passing laser light through glass plates that introduce random phase profiles (*13*, *14*) or by preparing randomly located dipoles in optical lattices, which are periodic potentials for atoms and molecules (*15*, *16*). In the former case, the strength of the disorder potential is determined by the laser intensity. In addition, atom-atom interactions can be set at a known value with collisional Fano-Feshbach resonances (*17*), which is then seen to modify the localization (*18*). For randomly distributed dipoles, the eigenenergy spectra are controlled by the lattice filling fraction.

Disordered quantum systems exhibit a rich diversity of statistical properties and types of extended and localized eigenstates. Within an early statistical classification scheme, they were often characterized by two distinct regimes where the nearest-neighbor spacing distribution (NNSD) of eigenvalues follows either an integrable Poissonian or a chaotic Wigner-Dyson statistics (*19*). At the same time, many dynamical systems that are partially chaotic exist (*20*, *21*). In this intermediate regime, the best fit to NNSD is obtained using Brody distribution (*22*). Recent studies of Anderson localization on hierarchical lattices, such as random regular graphs, point to the existence of a phase transition that separates extended-ergodic (EE) states from nonergodic, multifractal extended (NEE) states (*16*, *23*–*25*). The correspondence between classical and quantum chaos in ultracold dipolar collisions is examined by Yang *et al.* (*26*).

Here, we broaden the class of partially chaotic systems accessible for theoretical and experimental study. We present the search for fractal dimensions in disordered scattering and near-threshold bound states of ultracold highly anisotropic bosonic and magnetic ^{164}Dy lanthanide atoms in an external magnetic field. Their high-density Fano-Feshbach resonance spectrum was recently seen in a number of experiments with dysprosium (*8*, *27*, *28*) and is without precedent in other ultracold quantum gases. For comparison, the resonance density in alkali-metal atom collisions is two orders of magnitude lower.

We concentrate on the hundreds of loosely bound dysprosium molecular states, with binding energies less than 0.5 GHz when expressed in units of the Planck constant *h*, as we increase the strength of an external magnetic field *B* from 0 to 250 G, where 1 G equals 0.1 mT. Unique is the notion of locality and the role of scattering thresholds. In the Anderson Hamiltonian with random onsite energies, this is defined through the location of the particle or spin. Eigenstates are superpositions of being at different locations. In our system, the zero–*B*-field molecular eigenstates are the natural equivalent (local) states. These *B* = 0 states have energies either below or above the thresholds. Finite *B*-field eigenstates are then superpositions of zero-field eigenstates. This is schematically illustrated in Fig. 1. The *B* = 0 Hamiltonian matrix (Fig. 1A) is diagonal with values that are randomly distributed. For nonzero magnetic field, off-diagonal matrix elements appear because of the Zeeman interaction and are shown as yellow boxes in Fig. 1B. Our Hamiltonian matrix has a few hundred bands.

By focusing on the mixing coefficients of molecular wave functions, we find a transition between two NEE regimes as we increase the applied magnetic field. We also isolate localized states within our chaotic level structure. The transition is further elucidated by computing the spectrum of multifractality using the concepts of De Luca *et al*. (*23*) and correlated with statistical descriptions of the nearest-neighbor energy spacing distributions of the near-threshold quantum states. Our previous analyses (*8*, *9*) have already shown that the spacing distribution of these resonances closely, but not completely, follows a Wigner-Dyson distribution for the larger magnetic field strengths.

## RESULTS

The Hamiltonian of our diatomic ^{164}Dy quantum system in a homogeneous magnetic field with strength *B* was discussed by Maier *et al*. (*8*), and only a brief account of the features relevant in this discussion is presented here. The atoms have angular momentum with *i* = *a* or *b*, *j*_{i} = 8, and projection *m*_{i} along the magnetic field. For *B* = 0, the total molecular angular momentum and parity (that is, even versus odd relative orbital angular momentum or partial wave) are conserved. For even parity, the Hamiltonian describes coupled rovibrational motion in 81 distinct gerade anisotropic potentials that depend on the interatomic separation and atomic spins. A schematic of the molecular potentials, scattering thresholds, and a rough estimate of the level density in the relevant energy region is shown in Fig. 2. A sensible choice of basis functions is the weakly bound and above-threshold continuum *B* = 0 eigenstates |ν*JM*〉 with eigenvalue . For nonzero field and fixed *M*, unit-normalized eigenfunctions are then(1)with expansion coefficients *c*_{νJM}. The coefficients are determined by diagonalizing the Hamiltonian, with diagonal matrix elements coupled (and shifted) by the matrix elements of the Zeeman Hamiltonian. Crucially, the coupling strengths are proportional to *B* and only nonzero when |Δ*J*| = 0, 1. Hence, the Hamiltonian leads to a banded matrix that is schematically depicted in Fig. 1. Figure 3 shows the exact eigenenergies of our Hamiltonian for four typical *B*-field strengths. The numerical approach and convergence properties of the Hamiltonian are described in Materials and Methods.

We analyze the fractal properties of eigenfunctions of near-threshold bound states of ^{164}Dy_{2} that can lead to the magnetic Fano-Feshbach resonances, when their binding energy goes to zero for increasing *B*. Our analysis follows the spectral statistics and the spectrum of fractal dimensions (SFDs) as used by De Luca *et al.* (*23*). We define the generalized moment or partition function as well as a grand potential for each eigenfunction |Ψ; *M*〉 as(2)respectively, where is a generalization of the inverse of a temperature, and for simplicity, we suppressed the eigenstate labels Ψ and *M* and replaced the indices ν*JM* by the single index *i*. Here, log_{N} is the base-*N* logarithm. Unit normalization of the eigenstates ensures that *I*(*q* = 1) = 1 and τ(*q* = 1) = 0. The function *I*(*q*) can also be interpreted as a generalization of the inverse participation ratio (IPR) to which it coincides for *q* = 2.

The functions *I*(*q*) and τ(*q*) classify eigenstates of the disordered Hamiltonian. For example, a localized eigenstate with a single dominant *c*_{i} has *I*(*q*) ≈ 1 and τ(*q*) ≈ 0 independent of *q* and *N*. On the other hand, near-equal population with |*c*_{i}|^{2} ≈ 1/*N* for all basis states *i* leads to *I*(*q*) ≈ *N*^{1 − q} and τ(*q*) ≈ *q* − 1. Hence, the grand potential τ(*q*) does not depend on the number of basis functions. The set of eigenstates |Ψ; *M*〉 within a small eigenenergy interval that still contains a fairly large number of eigenstates is ergodic when the “local” (that is, fixed *i*) average 〈|*c*_{i}|^{2q}〉 over such a subset of eigenstates approaches 〈*I*(*q*)〉/*N* when *N* → ∞. The system is thus said to uniformly explore all basis functions with relatively weak fluctuations. On the other hand, a departure from these relations and, in particular, τ(*q*) ≠ *q* − 1 is a signature of nonergodic eigenstates.

Figure 4 (A to D) shows the grand potential τ(*q*) for Dy_{2} eigenstates |Ψ; *M*〉 with eigenenergy in at four finite *B*-field values and GHz. All curves are monotonic and concave, meet at *q* = 0 and 1, and, for *q* > 1, approach a linear functional form with a slope that is less than one. The curves can be compared to τ(*q*) = 0 for a localized state (or any of our eigenstates at *B* = 0) and τ(*q*) = *q* − 1 for ergodic GOE states. Here, GOE is an abbreviation for the Gaussian orthogonal ensemble of *N*-dimensional Hamiltonians, whose matrix elements are randomly chosen from a Gaussian distribution. Averages over eigenstates within a small energy window for the whole ensemble of these GOE Hamiltonians satisfy τ(*q*) → *q* − 1 for *N* → ∞. In the context of metal-to-insulator transitions (*29*), 〈τ(*q*)〉 has been used to distinguish insulator phases with 〈τ(*q*)〉 = 0 from metallic states with 〈τ(*q*)〉 ≠ 0. Figure 4 (A to D) also shows that τ(*q*) depends not only on near-threshold eigenstates but also on magnetic field. At *B* = 10 G, with already significant Zeeman couplings between zero-field basis states, a relatively large number of localized states with a slope close to zero for *q* > 1 remains. For stronger *B* fields, mixing increases and more eigenstates have a finite slope. At *B* = 250 G, only one state remains localized, whereas τ(*q*) for all other states seem to collapse onto a single universal curve with a critical slope that is less than one for *q* > 1. We identify this as the unique universal characteristic of disorder in our system.

We further quantify these observations by studying the fractal dimension *D*_{q} = 〈〈τ(*q*) − τ(1)〉〉/(*q* − 1). Here, 〈〈 ⋯ 〉〉 indicates a double average over eigenfunctions with energy in and over realization of the molecular Hamiltonian. In particular, we changed the strength of the strongest anisotropic potential contribution, as identified by Maier *et al*. (*8*), in the zero-field Hamiltonian over ± 5 % from its physical strength assuming a uniform probability distribution. This range of strengths is sufficiently large that the *B* = 0 bound state spectrum changes significantly but, at the same time, is small enough to not change its mean spacing and random distribution. (A ± 5 % variation also reflects the uncertainty in this strength.) As a consequence, the accuracy of our statistical analysis improves without changing our conclusions.

Figure 4E shows the computed *D*_{q} as a function of *B* for *q* = 1 and 2. For *q* → 1, *D*_{q} corresponds to the derivative of averaged τ(*q*) with respect to *q* and, in fact, is an effective entropy(3)of the probabilities |*c*_{i}|^{2}. By construction, the entropy is 0 ≤ *S* ≤ 1 and only equals one when for all *i*. By construction, the domain of the entropy is 0 ≤ *S* ≤ 1 and *S* only equals one when for all *i*. The choice of *D*_{q = 2} is a compromise. It gives a reasonable estimate of the large *q* trends and, at the same time, describes the dimensionality near the IPR point.

Two observations can be made about Fig. 4E. First, we note that the fractal dimension *D*_{q} depends on *q*. This is a characteristic of a multifractal system. Second, at fixed *q*, the entropy *S* and *D*_{2} linearly increase with magnetic field starting from zero, that is, the system becomes chaotic and then saturates at a finite value near 0.5 much smaller than one. The transition between the two behaviors occurs between 50 and 100 G. Figure 4E also shows the square root of the variance of the fractal dimension as a function of *B*. A smaller variance indicates increased universality with a larger fraction of the eigenfunctions with a similar τ(*q*). It is smallest for *B* > 200 G.

A measure of the multifractality of *D*_{q} or, equivalently, the *q* dependence of 〈〈*τ*(*q*)〉〉 is given by the SFD or the Legendre transform of 〈〈τ(*q*)〉〉 (*29*). It is defined as and implies *d*〈〈*τ*(*q*)〉〉/*dq* = α. Figure 5A shows *f*(α) as functions of α for six representative *B* values ranging from 1 to 250 G. Each curve is defined between (*α*_{−}, *α*_{+}) with 0 < *α*_{−} < *α*_{+}. They are positive and resemble half of an ellipse with a maximum value *f*_{max} ≈ 1 at α = α_{max}, the slope of *τ*(*q*) at *q* = 0. A small value of *f*(α) indicates that corresponding tangent lines of *τ*(*q*) extrapolate to the origin. Finally, *α*_{−} is the slope of *τ*(*q*) when *q* → + ∞.

We characterize the behavior of *f*(α) as multifractal with a trend toward the ergodic limit for increasing magnetic field strength, where the ergodic limit corresponds to a delta function at α = 1. For increasing *B*, the SFD sharpens up with a maximum that shifts to smaller α and approaches α_{max} = 2.3. The approach is shown in Fig. 5B. The SFD reaches a universal function for *B* > 50 G, leading to the conclusion that weak multifractality persists over the whole range of magnetic fields accessible in our numerical analyses. In contrast to the results of De Luca *et al*. (*23*) for a Bethe lattice, we do not obtain a triangular-shaped SFD even for small magnetic field strengths. The shape of our curves is always flat-topped, typical for weak multifractality. The lack of ergodicity in the collisional dynamics of Dy atoms will have significant implications for the observed spectra and invites further investigation.

The onset to universality has its origin in our Hamiltonian matrix structure. It contains the *B* = 0 G eigenenergies on the diagonal, with mixing and shifts induced by the Zeeman Hamiltonian with each of its matrix elements proportional to *B*. Naively, we can then expect that the expansion coefficients of eigenstates |Ψ, *M*〉 become independent of *B* when the energy scale of the Zeeman interaction is larger than that of the *B* = 0 G eigenenergies. However, in practice, the energy scales for the two contributions are hard to estimate because of the presence of the (discretized) continua near our threshold bound states with energy between . The coupling strength between the bound and spatially extended continuum states has a different character than that between bound *B* = 0 G basis states. We have only been able to empirically determine this *B* value (that is, between 50 and 100 G) from the spectrum of the fractal dimensions.

We finish by connecting the statistical properties of eigenstate wave functions with those of the NNSD of their eigenenergies. This spectral statistics of quantum systems is often characterized by either a Poisson distribution *P*_{P}(Δ) = *e*^{−Δ} or a Wigner-Dyson distribution (*30*), where Δ = *s*/〈*s*〉 is the dimensionless scaled spacing between levels. The Brody distribution is used to describe quantum systems that experience a transition from Poisson to Wigner-Dyson statistics (*21*). Here, we construct the NNSD as a histogram from eigenenergy spacings within the energy interval and over several realizations of the molecular Hamiltonian, as described in the previous subsection. We then determine the Brody parameter η that provides the best fit of the Brody distribution (*22*) to the histograms, where *b* is a known constant such that for *k* = 0 and 1.

Figure 6 shows our Brody parameter η as a function of *B* between 0 and 250 G. Examples of our fitting of the Brody parameter to our histograms for two *B* fields, 25 and 160 G, are presented as insets. We observe that the Brody parameter changes monotonically from η ≈ 0.1 for *B* ≥ 1 G to an η just below 1 for *B* > 100 G. Hence, our system undergoes a transition from an integrable system with Poissonian distributed levels to a “hard-chaotic” system that satisfies the Wigner-Dyson distribution. Crucially, we observe two physical regimes: For *B* < 50 G, the Brody parameter increases linearly, whereas for larger field strengths, it saturates to an asymptotic value and becomes universal. The critical or transition *B* value is ~60 G, as defined in the figure, and closely follows the universality as observed in the fractal dimension and multifractality of Figs. 4 and 5.

## DISCUSSION

Here, we obtained a detailed understanding of the chaotic quantum dynamics of two weakly bound bosonic lanthanide ^{164}Dy atoms in a magnetic field. The near-threshold chaotic molecular structure is a consequence of the magnetic Zeeman interaction–induced mixing of 81 gerade electronic potentials of Dy_{2} that dissociate to two ground-state Dy atoms (*31*). The remarkably complex nature of the electronic structure is due to the large total and orbital angular momentum of the dysprosium atom.

Our work was motivated by successful experiments in producing quantum degenerate gases of dysprosium atoms (*3*, *8*) as well as the measurement of its dense magnetic Fano-Feshbach spectrum. The latter is directly related to the collisional properties of atoms in an optical trap. In a previous study of chaos in this system (*8*), we showed that the NNSD of Fano-Feshbach resonances was closer to the Wigner-Dyson distribution than to the Poisson distribution.

Here, we identified the key parameter that governs chaos and randomness of the vibrational diatomic wave functions and how it changes with increasing magnetic field to a system with a higher level of chaoticity. We analyzed the fractal properties of the Dy_{2} eigenfunctions in the near-threshold region, which allowed us to characterize the degree of localization and ergodicity for each eigenstate. We were surprised to find that some states in the chaotic system are localized, whereas others are nonlocalized, making our system multifractal. Moreover, by increasing the magnetic field strength, we are able to change our system from partially chaotic to completely chaotic with a remarkable degree of universality. At the same time, we realized that the degree of ergodicity, characterized by fractal dimensions *D*_{q}, has a limiting value of ≈0.5 for large *B*, where complete ergodicity corresponds to *D*_{q} = 1.

## MATERIALS AND METHODS

Our Hamiltonian contains bound states as well as a continuum spectrum separated by a threshold at energy *E* = 0 corresponding to two ground-state Dy atoms at rest and in their energetically lowest magnetic Zeeman level. Following the study of Maier *et al*. (*8*), it suffices to study even parity states with *M* = − 16, limit the sum over *J* to ≤ *J*_{max}, and use *B* = 0 G eigenstates with , where *E*_{min} < 0 < *E*_{max}. The continuum spectrum is “discretized” by limiting the extent of the vibrational motion to a large radius chosen such that weakly bound *B* = 0 states are unaffected. The number of basis functions *N* is then finite. We use *J*_{max} = 36, *E*_{min}/*h* = − 120 GHz, *E*_{max}/*h* = + 60 GHz, and a largest radius of 52.9 nm. Thus, Eq. 1 has *N* = 46,300 coefficients. Approximately 7% of the matrix elements of our banded matrix are nonzero, and ≈15% of the basis states have . The value of matrix elements between continuum states is a few orders of magnitude smaller than those between bound states.

Our statistical analysis looks at eigenstates |Ψ; *M*〉 with energy in the small energy interval with and . The energy of these eigenstates has converged with respect to the basis size. We chose GHz to correspond to the nominal Dy Zeeman energy, μ_{B}*B*, at the largest *B*-field strengths studied. Here, μ_{B} is the Bohr magneton and MHz/G. There are approximately 150 basis functions, with energy corresponding to a mean spacing of 〈*s*〉/*h* = 3.3 MHz. Their NNSD is Poissonian (*8*). Energy shifts induced by the diagonal Zeeman matrix elements are already larger than this mean spacing when *B* = 1 G, a field that is smaller compared to our typical field strength.

We emphasized that the matrix size of the lanthanide Hamiltonian does not play the same role as it typically does in the conventional random matrix Hamiltonians (*32*). For random matrices, the eigenenergies do not converge to fixed values with increasing *N*. The statistical properties of eigenpairs as a function of *N* are often an integral part of a study and rely on the notion that basis functions are equivalent (that is, adding spins in a Heisenberg model or lattice sites in an atomic Hubbard model). This is not true for our Hamiltonian because basis states have a distinct physical interpretation. We can increase *J*_{max} and thus *N*; however, there are diminishing returns from these states, because both the eigenenergies converge and the additional number of zero-field states within the energy range [*E*_{min}, *E*_{max}] decreases.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank B. L. Altshuler for an insightful discussion.

**Funding:**Work at the Temple University is supported by the U.S. Air Force Office of Scientific Research and the Army Research Office (grants FA9550-14-1-0321, W911NF-17-1-0563, and W911NF-14-1-0378). Work at the Joint Quantum Institute is supported by the U.S. NSF (grant PHY-1506343).

**Author contributions:**S.K. and E.T. performed theoretical analysis and wrote the manuscript. C.M. and M.L. performed numerical simulations, contributed to writing the manuscript, and claim responsibility for the figures.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper. Additional data related to this paper can be requested from the authors.

- Copyright © 2018 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).