Research ArticlePHYSICS

Fractal universality in near-threshold magnetic lanthanide dimers

See allHide authors and affiliations

Science Advances  16 Feb 2018:
Vol. 4, no. 2, eaap8308
DOI: 10.1126/sciadv.aap8308


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.


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 (18) are paving the way toward these goals.

Most lanthanides have a partially filled submerged 4f-electron shell, which lies beneath a filled 6s2 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 (811), 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, 2325). 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 164Dy 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.

Fig. 1 Structure of our banded Hamiltonian for lanthanide dimers in a magnetic field with strength B.

Gray boxes correspond to matrix elements with a zero value. Others correspond to elements with a nonzero value. (A) For B = 0 G, the Hamiltonian is diagonal. Boxes with the same color are for states with the same total molecular angular momentum Embedded Image. (B) For finite B, the Hamiltonian is banded. The nonzero off-diagonal matrix elements (yellow boxes) are proportional to B.

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.


The Hamiltonian of our diatomic 164Dy 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 Embedded Image with i = a or b, ji = 8, and projection mi along the magnetic field. For B = 0, the total molecular angular momentum Embedded Image and parity (that is, even versus odd relative orbital angular momentum Embedded Image 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 Embedded Image. For nonzero field and fixed M, unit-normalized eigenfunctions are thenEmbedded Image(1)with expansion coefficients cνJM. The coefficients are determined by diagonalizing the Hamiltonian, with diagonal matrix elements Embedded Image 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.

Fig. 2 Asymptotic thresholds in the 164Dy + 164Dy system.

Schematic of (adiabatic) molecular potentials for B = 50 G as functions of separation R. For graphical purposes, partial waves are restricted to even ℓ up to 10. For R → ∞, the potentials dissociate to scattering thresholds labeled by ma + mb = −16,−15,−14,⋯ starting from the lowest limit. Bound states used in the statistical analysis lie within the energy range Embedded Image and 0 GHz just below the energetically lowest threshold.

Fig. 3 Near-threshold eigenstate energies at four typical B-field values as indicated below each of the panels.

We analyze the fractal properties of eigenfunctions of near-threshold bound states of 164Dy2 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〉 asEmbedded Image(2)respectively, where Embedded Image 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, logN 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 ci has I(q) ≈ 1 and τ(q) ≈ 0 independent of q and N. On the other hand, near-equal population with |ci|2 ≈ 1/N for all basis states i leads to I(q) ≈ N1 − 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 〈|ci|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 Dy2 eigenstates |Ψ; M〉 with eigenenergy in Embedded Image at four finite B-field values and Embedded Image 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.

Fig. 4 Grand potential and fractal dimension of the weakly bound Dy2 wave functions.

(A to D) τ(q) as a function of q for all weakly bound eigenstates (dark red curves) within an Embedded Image MHz energy range below the threshold at B = 10, 100, 200, and 250 G, respectively. The corresponding eigenenergies for these four B fields are shown in Fig. 3. The two dashed straight lines correspond to τ(q) of the ergodic GOE and localized limits, respectively. (Although not shown, the curves remain concave for q < 0. (E) Fractal dimension Dq(B), defined in the text, as functions of B for q = 1 and 2. The error bars indicate the spread of fractal dimension. Solid lines are for guiding the eye. For q = 1 and 2, the parameter d0 = 0.57 and 0.49, respectively.

We further quantify these observations by studying the fractal dimension Dq = 〈〈τ(q) − τ(1)〉〉/(q − 1). Here, 〈〈 ⋯ 〉〉 indicates a double average over eigenfunctions with energy in Embedded Image 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 Dq as a function of B for q = 1 and 2. For q → 1, Dq corresponds to the derivative of averaged τ(q) with respect to q and, in fact, is an effective entropyEmbedded Image(3)of the probabilities |ci|2. By construction, the entropy is 0 ≤ S ≤ 1 and only equals one when Embedded Image for all i. By construction, the domain of the entropy is 0 ≤ S ≤ 1 and S only equals one when Embedded Image for all i. The choice of Dq = 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 Dq depends on q. This is a characteristic of a multifractal system. Second, at fixed q, the entropy S and D2 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 Dq or, equivalently, the q dependence of 〈〈τ(q)〉〉 is given by the SFD or the Legendre transform of 〈〈τ(q)〉〉 (29). It is defined as Embedded Image 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 fmax ≈ 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 → + ∞.

Fig. 5 Multifractality of Dy2 wave functions.

(A) Spectrum of the fractal dimension f(α) for Dy2 eigenstates within a 500-MHz energy window below the threshold as a function of α for B = 1, 10, 50, 100, 175, and 250 G. The dashed vertical line indicates the ergodic limit, where f(α) is a delta function located at α = 1. Filled circles for selected B locate the maximum of f (α). (B) Location αmax of the maximum of f(α) as a function of B.

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 Embedded Image. 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 PP(Δ) = e−Δ or a Wigner-Dyson distribution Embedded Image (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 Embedded Image 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 Embedded Image (22) to the histograms, where b is a known constant such that Embedded Image 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.

Fig. 6 Brody parameter η extracted from the NNSD as a function of B for weakly bound states with binding energy less than Embedded Image GHz.

The dashed gray lines extrapolate the linear small and large B-field behavior. Their crossing defines a transition strength. Insets show the NNS probability distributions (blue filled circles) for B = 25 and 160 G as a function of the normalized energy spacing s/〈s〉, with mean spacing 〈s〉 obtained for each B. Solid blue lines are Brody distributions fit to the data. Gray shaded areas and red curves indicate the Wigner-Dyson (η = 1) and Poisson (η = 0) distributions, respectively.


Here, we obtained a detailed understanding of the chaotic quantum dynamics of two weakly bound bosonic lanthanide 164Dy 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 Dy2 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 Dy2 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 Dq, has a limiting value of ≈0.5 for large B, where complete ergodicity corresponds to Dq = 1.


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 ≤ Jmax, and use B = 0 G eigenstates with Embedded Image, where Emin < 0 < Emax. 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 Jmax = 36, Emin/h = − 120 GHz, Emax/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 Embedded Image. 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 Embedded Image with Embedded Image and Embedded Image. The energy of these eigenstates has converged with respect to the basis size. We chose Embedded Image GHz to correspond to the nominal Dy Zeeman energy, μBB, at the largest B-field strengths studied. Here, μB is the Bohr magneton and Embedded Image MHz/G. There are approximately 150 basis functions, with energy Embedded Image 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 Jmax 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 [Emin, Emax] 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.


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.

Stay Connected to Science Advances

Navigate This Article