Quantum chaos in Feshbach resonances of the ErYb system

We investigate ultracold magnetic-field-assisted collisions in the so far unexplored ErYb system. The nonsphericity of the Er atom leads to weakly anisotropic interactions that provide the mechanism for Feshbach resonances to emerge. The resonances are moderately sparsely distributed with a density of $0.1\,{\rm G}^{-1}-0.3\,{\rm G}^{-1}$ and exhibit chaotic statistics characterized by a Brody parameter $\eta \approx 0.5-0.7$. The chaotic behaviour of Feshbach resonances is accompanied by strong mixing of magnetic and rotational quantum numbers in near-threshold bound states. We predict the existence of broad resonances at fields $<300\,{\rm G}$ that may be useful for the precise control of scattering properties and magnetoassociation of ErYb molecules. The high number of bosonic Er-Yb isotopic combinations gives many opportunities for mass scaling of interactions. Uniquely, two isotopic combinations have nearly identical reduced masses (differing by less than $10^{-5}$ relative) that we expect to have strikingly similar Feshbach resonance spectra, which would make it possible to experimentally measure their sensitivity to hypothetical variations of proton-to-electron mass ratio.


Introduction
One of the most important recent achievements in the field of ultracold quantum gases was the production of samples of highly magnetic atoms: chromium (Cr) [1], erbium (Er) [2] and dysprosium (Dy) [3][4][5][6]. At large distances these atoms interact mostly through an anisotropic dipole-dipole interaction which decays more slowly than the usual van der Waals atom-atom interaction. As a result, a sample of highly magnetic atoms can exhibit a rich variety of phenomena distinct from the behaviour of systems of alkali-metal or alkaline-earth metal atoms. Highly magnetic atoms promise to be a key building block in accomplishing the vision of making a versatile toolbox for quantum simulations of condensed-matter physics and manybody phenomena [7]. Some of these promises have already been fulfilled experimentallyfor example, the extended Bose-Hubbard model was recently realized experimentally using ultracold Er atoms confined to a 3D optical lattice [8]. Quantum degenerate Er was used to study the deformation of a Fermi surface by the dipole-dipole interaction [9] and the formation of quantum droplets [10,11]. Apart from the applications in quantum many-body physics [12,13] the Dy atom is notably known for its role in studies of physics beyond the Standard Model [14,15].
Magnetic Feshbach resonances are one of the most important tools for controlling the interactions of ultracold atoms [16]. In collisions of open-shell lanthanides with a non-zero orbital quantum number the spectrum is very dense: the anisotropic interaction potentials couple the magnetic Zeeman states with end-over-end rotational quantum number, which results in multiple bound states being coupled to the entrance s-wave channel. Such a pattern was initially predicted by Kotochigova and Petrov [17,18], and confirmed experimentally for Er [9], and later for Dy atoms [19][20][21][22]. The resonance pattern is not only dense, but also features the characteristics of quantum chaos, evidenced by correlations in the resonance spacings.
The positions of magnetic Feshbach resonances appear for field values at which the molecular bound states cross the zero-energy threshold. Quantum chaos emerges when the mean spacing between these bound states is comparable to the coupling strength between them. In this regime the couplings cause multiple avoided crossings to emerge which results in a spaghetti-like pattern of repelling energy levels and a highly congested spectrum of overlapping resonances. The resonances can not be interpreted individually, but their statistical behaviour can be analyzed instead (see pioneering work of Porter and coworkers eg. [23,24]) and the tools for such analysis include the Random Matrix Theory (RMT) developed by Wigner [25] and Dyson [26]. For ultracold gases the first evidence of chaos in a Feshbach resonance spectrum was reported in Er dimer by Frisch et al. [27]. Since then, the chaotic character of bound states and Feshbach resonances in ultracold collisions was also found for several other systems, for example, in mixtures of Yb atoms in 1 S 0 and 3 P 2 states at large magnetic fields [28], and molecular collisions of Li atom with CaH and CaF [29]. On the other hand, the resonance spacings in collisions of Er and Li atoms in a magnetic field [30], and if highly magnetic europium ( 7 S) with alkali metal atoms [31] have both been shown to follow the non-chaotic poissonian distribution.
Here we investigate the interactions and Feshbach spectra between weakly anisotropic lanthanides and heavy spin-singlet atoms using ErYb as a prime example. Ytterbium, alongside Sr, is the most widely used spin-singlet atom with applications for optical clocks, quantum gases and quantum simulation [32][33][34][35]. ErYb should also be very attractive due to the fantastic mass-scalability of both Yb [36,37] and Er (as shown here). Three of the bosonic Yb isotopes, 168 Yb [38] and 170 Yb [39] and 174 Yb [40] form stable Bose-Einstein condensates, and several Fermi and Bose-Fermi and Bose-Bose gases [41,42] have also been demonstrated experimentally. Er and Yb atoms could form many isotopic mixtures, including two fermion-fermion mixtures with nearly equal masses and very close polarizabilities (hence, trapping properties). While an Er-Yb quantum degenerate mixture has not been achieved yet, it is clearly within the reach of present state-of-the-art techniques. Our predictions could be tested by forming a dual Mott insulator of Er and Yb atoms. Keeping unity filling for Er would prevent losses due to Er-Er Feshbach resonances.
We unveil the impact of anisotropy on the chaotic behavior in this system. The interaction between these two atoms is by far simpler than in case of the lanthanide dimer or lanthanidealkali metal atom case and can be tackled by realistic ab initio calculations, which might help shed some light on much more complex systems involving lanthanides (mixed-lanthanide quantum gases, lanthanide+alkali metal atom). Given the progress in the field and possibility of mutual trapping of alike atoms the results presented in this paper should hold for most combinations, including (Ho, Dy, Tm)+(Sr, Yb, Hg).
The paper is planned as follows. In the next chapter we describe the interaction between the erbium and ytterbium atoms focusing on details of ab initio calculations and their implications. Next we describe the methodology of solving the close-coupling equations in the magnetic field. In Section 4 we analyse the mass-scaling behavior of the zero-field scattering lengths and show example ErYb magnetic Feshbach spectra. Next, we focus on the statistical properties of resonance positions showing evidence of quantum chaos. In Section 6 we show that the weak anisotropy present in the interaction between Er and Yb atoms is  (d) Anisotropic Legendre coefficients V k for the short-range interaction potential based on the broken symmetry X 3 Σ − state computed using the CCSD(T) method with added anisotropy. responsible for its chaotic behaviour. Finally we look into the sensitivities of magnetic Feshbach resonances to small changes of the proton-to-electron mass ratio that, in ErYb, could be measured experimentally thanks to the existence of two isotopic combinations with nearly equal reduced masses. Section 8 concludes the paper.

Interaction of Er and Yb atoms
The erbium atom in its 3 H 6 ground state is slightly anisotropic due to two missing electrons in the 4 f shell screened by the doubly occupied 6s 2 shell (4 f 12 6s 2 ). Since Yb has a spherical 1 S ground state, the ground state ErYb molecule, shown schematically Fig. 1(a), inherits its quantum numbers from the erbium atom. The projection Λ of the erbium atom's large orbital angular momentum (l = 5) onto the interatomic axis gives rise to eleven states labeled by Λ = −l . . . l. The (Λ = 0) 3 Σ − state is unique while the remaining five |Λ| > 0 states are doubly degenerate.
Electronic structure calculations for such a heavy system as the ErYb molecule are challenging. Our general strategy follows Ref. [30] and involves the use of two complementary ab initio methods. We first obtain the differences between complete active space self consistent field (CASSCF) total energies of all 11 states [43,44]. Since the CASSCF method itself does not reproduce the dynamic correlation properly, which results in a lack of the attractive dispersive interaction, we then add the calculated CASSCF anisotropies to the broken symmetry spin-restricted CCSD(T) lowest potential [45,46]. At long ranges we use van der Waals coefficients for isotropic part of the potential, C 6,0 , and the leading anisotropy, C 6,2 .
For all of our electronic structure calculations we used the Molpro quantum chemistry code [47]. Quasi-relativistic effective core potentials ECP28MWB and ECP60MDF for the first 28 and 60 electrons in Er and Yb atoms, respectively [48,49], were used to replace the core electrons. The ECP28MWB gaussian basis set for Er included the [14s13p10d8 f 1g] functions with uncontracted s and p functions, and an additional g = 4.0 exponent. For Yb, we used the [7s7p6d8 f 1g] basis set in which the f and g functions were added from the Er atom. The basis set was also augmented by midbond functions to improve the accuracy in the potential well region. Finally, the description of both atoms was modified by adding one shell of even-tempered diffuse functions per orbital angular momentum.
In their early studies of interactions of the lanthanides, Buchachenko and coworkers [50][51][52] demonstrated that the interaction has a peculiar character: the atomic anisotropy is suppressed due to the shielding by the doubly occupied outermost electronic shell. Systems like Tm+He and Er( 3 H 6 )+Li( 2 S 1/2 ) have been previously studied using CASSCF [30,43,44,51]. This method preserves the symmetry associated with the projection Λ of the orbital angular momentum l on the interatomic axis; therefore, it allows to compute potentials corresponding to pure |Λ| states. In our study of the Er( 3 H 6 )+Yb( 1 S 0 ) system, we have performed CASSCF calculations with an active space including 4 f 6s6p orbitals for the Er atom, and 6s6p orbitals for the Yb atom. We carefully monitored the orbital symmetry of the optimized states by controlling the projection of orbital angular momentum on the molecular axis. In Fig. 1(c) we show the differences between |Λ| > 0 potentials with respect to the state with the lowest energy, 3 Σ − state. The potential anisotropies corresponding to |Λ| > 0 states range from a few to several dozens of cm −1 and increase with |Λ|.
The CASSCF interaction potentials do not reproduce the dynamic correlation correctly which results in an absence of dispersion forces and, as a consequence, incorrect long range asymptotics. To fix this we performed spin-restricted CCSD(T) calculations, which for similar systems provides few-percent accuracy [53]. However, we were only able to converge the CCSD(T) calculations for four states (instead of all eleven), and we were not able to control the |Λ| quantum number in them. Hence we used the interaction energy for the lowest state obtained with CCSD(T), shown in Fig. 1(b), combined with the energy differences obtained with CASSCF to model the atomic interaction.
As a sanity check, we have also performed broken-symmetry MP2 calculations based on the restricted open-shell Hartree-Fock determinant, as well as CASPT2 calculations [54]. The spectroscopic parameters of potential energy curves obtained with CCSD(T), MP2 and CASSCF methods are shown in Table 1. The spectroscopic constants predicted by all three methods are quite similar. The dissociation energy for ErYb obtained with broken-symmetry CCSD(T) is 639 cm −1 . This value can be compared with experimental estimates of the dissociation energy for Yb 2 of 738 cm −1 [55] and RbYb of 785 cm −1 [56]. In the latter case, ab initio calculations based on the ECP60MDF core potential as reported in Ref. [56] predicted a well depth of 704 cm −1 , hence we expect our result could be underestimated by about 10%.
At long ranges the interaction between Er and Yb atoms is dominated by the van der Waals interaction potential given by V(R, θ) = −C 6,0 R −6 −C 6,2 R −6 P 2 (cos θ). Higher anisotropy Table 1. Spectroscopic constants for the lowest interaction potential of the Er( 3 H)+Yb( 1 S) system including equilibrium distance R e , well depth D e , rotational constant B rot and harmonic frequency ω 0 . The latter two correspond to the 166 Er+ 174 Yb isotopic pair composed of the two most abundant isotopes. terms decay much faster with distance (at least R −8 ) and are thus ignored in this work. We have estimated the C 6,0 coefficient by means of a combination rule proposed by Tang [57] where C 6,0 can be obtained from static polarizabilities of atoms α Er , α Yb and homonuclear coefficients C 6,0,Er 2 and C 6,0,Yb 2 : . (1) The anisotropic C 6,2 coefficient can be evaluated by scaling the isotropic C 6,0 by the ratio of the erbium atom's tensor and scalar polarizabilities: The necessary atomic polarizabilities and homonuclear C 6,0 factors can be found in Refs. [36,[58][59][60][61]. For ErYb we find C 6,0 = 1835 E h a 6 0 and C 6,2 = 33.62 E h a 6 0 (the Hartree energy E h ≈ 4.359744650(54) × 10 −18 J and Bohr radius a 0 ≈ 5.2917721067(12) × 10 −11 m are the atomic units of energy and distance).

Quantum scattering calculations
While Yb, in its 1 S 0 ground state, is a structureless atom, Er has a highly magnetic 3 H (l = 5, s = 1) ground state which splits into three fine sublevels j = 4, 5, 6 with the 3 H 6 state being the lowest in energy. Bosonic isotopes of both Er and Yb also lack hyperfine structure. Since the fine splitting of more than 6000 cm −1 between the 3 H 6 ground state and the closest excited state, 3 H 5 , is an order of magnitude larger than the potential depths, we carry out our scattering calculations for just the lowest 3 H 6 + 1 S 0 asymptote. Our basis spans the Hilbert space of consists of, respectively, kinetic and rotational energy, the atomic interaction potentialsV and Zeeman interaction of the isolated Er atom with the external magnetic field B,Ĥ Z = µ B Bg jĵz . The Landé g-factor g j = 1.166; R and µ denote the interatomic distance and reduced mass.
It is standard procedure for the interaction potential operatorV to be expanded in Legendre polynomials [62][63][64]: where k = 0, 2, ...2l. The matrix elements of the angular part P k (cos θ) in the j-coupled basis set are where m k = M L − M L , symbols in parentheses and curly braces are 3j and 6j symbols [65]. The potential matrix changes the m l quantum number such that m l + M L is conserved and the mixing between different m j exists for each element with nonzero k.
To be more precise, k > 0 directly couples the channels with ∆m j = ±1, ±2, . . . , ±k at the cost of M L quantum number. The higher anisotropies k = 2, 4 . . . were obtained using the following transformation [66,67] Fig. 1(d) shows the anisotropic Legendre coefficients for erbium-ytterbium interaction for k = 2 . . . 8. All calculations with the collisional Hamiltonian were done with the close-coupling (coupled channel) code implemented by Hutson and LeSueur in bound (for the calculations of bound states in the magnetic field), field (for converging the value of magnetic field at which the bound state crosses the threshold) [68], and molscat [69] programs (for scattering length). We used the hybrid log-derivative-Airy method of Alexander and Manopoulos [70] with a fixed step size of 0.002 a 0 for R between 6.2 a 0 and 30 a 0 and a variable step size for R between 30 a 0 and 800 a 0 . The collision energy was E/k B = 10 −8 K.
The anisotropic interaction potentialV couples the initial s-wave (L = 0) scattering channel to many closed channels that correspond to higher partial waves. A practical calculation, however, can only include a finite number of partial waves, and the resulting spectrum of Feshbach resonances will depend on the selected basis set. Generally, as the number of included partial waves is increased, the total number of Feshbach resonances also grows. The width, however, of resonances originating from higher partial waves quickly becomes vanishingly small, as the coupling of the s-wave to higher partial waves becomes weaker, and the density of states of closed channels at the threshold decreases as L quantum number increases. We have found that the number of resonances with an appreciable width -larger than 0.01 G -saturates quickly and for the remainder of the paper we calculated the properties of resonances which include L up to 20: a comparison with calculations performed for maximum L = 50 for the magnetic fields up to 100 G shows that the calculations with basis set L = 50 produces 10 more resonances (48 compared to 38) but all new resonances have a width smaller than 0.01 G. Scattering lengths for fully anisotropic interaction. The shaded area corresponds to a full cycle of scattering lengths. The light grey vertical lines correspond to reduced masses of ErYb isotopic combinations. The red and blue lines correspond to the respective "off-resonance" and "on-resonance" cases discussed later in text.

Zero-field scattering lengths and Feshbach spectra
Due to the large reduced mass of the Er+Yb system, our isotropic potential V 0 supports a large number of bound states, between 64 for the lightest and 66 for the heaviest isotopic combinations. In absence of anisotropy, the mass variation of the zero-field s-wave scattering length [ Figure 2(a)] follows the well known periodic behaviour, a(B = 0) ≈ a 1 − tan (φ − 3π/8) , where φ = −1 ∞ R t −2µV(r )dr is the WKB phase integral and a = 2 −3/2 [Γ(3/4)/Γ(5/4)](2µC 6 / 2 ) 1/4 ≈ 73 a 0 is the "mean scattering length" [71]. The range of reduced masses covered by all isotopic combinations (vertical lines in Figure 2(a)) covers approximately one and a half cycle of a. Figure 2(b) shows the scattering length once the anisotropic interaction is switched on. The anisotropic V k terms couple the s-wave to channels corresponding to higher partial waves causing additional (Feshbach) resonances to emerge. As result we obtain a complex, yet periodic structure. The resonance pattern, as we will show, is in itself chaotic. The periodicity, however, may stem from the near-threshold physics being governed mainly by the long range van der Waals interaction, and the WKB phase, φ. The length of this cycle, of 2.6 u, aligns perfectly with the length of the cycle in V 0 alone.
Since we cannot predict the scattering lengths ab initio, we instead study the Feshbach spectra for two representative cases, which we shall dub "off-resonance" and "on-resonance". The "off-resonance" case, shown in Fig. 2 in red, corresponds to a zero-field s-wave scattering length close toā, which should be representative for most isotopes. For the much fewer isotopes that lie on-resonance, we also calculate a spectrum for a reduced mass that for our model coincides with a resonance (Fig. 2, blue).   Figure 3 shows the magnetic Feshbach spectra, resonance positions and near-threshold bound states as a function of the magnetic field B. For clarity we only show the range of B ∈ [0, 300] G, although we have calculated the spectra up to 1000 G. The resonance spectrum is dense, with many overlapping resonances. The resonance density is significantly higher at small fields: the full 1000 G scan yields 168 and 173 resonances for the "offresonance" and "on-resonance" cases, respectively, but as much as 31 and 38 are contained within the first 100 G. While the overall densities of the spectra between the two samples are similar, no individual features from one spectrum can be discerned in another. The complexity of the Feshbach spectrum originates from the chaotic behavior of near-threshold bound states characterised by multiple avoided crossings caused by the anisotropic interaction.
Recently, Frye et al. [72] published a procedure for extracting parameters, namely the positions B res , widths ∆B, and background scattering lengths a bg , of Feshbach resonances that conform to the Breit-Wigner profile, We have found it applicable to the majority of the resonances in our spectra. Histograms of the resulting widths (in log scale) are shown in Figures 4(a) and (b). In both cases the vast majority of resonances is narrower than 1 G and about 10% is narrower than 1 mG. The respective median widths for the "off-resonance" and "on-resonance" datasets are nearly identical: ∆B = 0.241 G and ∆B = 0.246 G. The distribution quartiles are qualitatively similar: the middle 50 % of resonances are contained within the ranges of 0.026, 0.800 G and 0.126, 0.659 G, respectively. Narrower resonances originate from higher partial waves than broad resonances. We illustrate this by evaluating the differences in magnetic moments, ∆m eff j g j µ B , between the initial scattering channel and the bound states the resonances originate from. Since the interaction potential only couples channels with matching m j + m l , ∆m j determines ∆m l , which in turn provides indirect information on the partial wave composition of the bound state. Bound states are, in principle, mixtures of many partial waves. While we can not determine the actual composition of a bound state in terms of pure m j states, we can determine its magnetic moment ∆m eff j g j µ B = ∂E/∂B. In practice, we evaluate this derivative by probing the resonance position upon small changes to threshold energy E. Fig. 4(b) shows log resonance widths as a function of ∆m eff j . The correlation coefficients ρ = −0.81 and ρ = −0.78 for the respective "off-resonance" and "on-resonance" datasets show clear correlation between (log 10 ) resonance widths and ∆m eff j .

Resonance statistics and evidence of quantum chaos
Our calculated magnetic Feshbach spectra exhibit a chaotic behaviour. To show this, we quantify the resonance statistics using two standard measures of chaos: the Brody parameter η, which probes the propensity of adjacent resonances to repel each other [73], and resonance number variance, a measure of long range resonance correlations [74]. In our analysis of the chaotic behaviour of resonance spacings we take into account the slow variation of resonance density by following the standard "unfolding" procedure [28,75]. Figures 5(a) and (b) show a staircase function S (B) = i Θ(B i ), where Θ(x) is the Heaviside step function and B i are the resonance positions for the "off-resonance" and "on-resonance" datasets, respectively. S (B) can be understood to give the total number of resonances for fields from zero up to B. We fit a smoothing spline, ξ(B) to the staircase function. The function ξ(B) is, by design, a continuous function of B and can be viewed as a "generalized" resonance number with respect to B. Its values at resonance positions provide the unfolded resonance positions, ξ i = ξ(B i ). The derivative of ξ with respect to the magnetic field, ∂ξ/∂B gives an estimate of the local resonance density, shown in Fig. 5(e). In Er+Yb the resonance density for low fields B < 100 G is slightly above 0.3 G −1 , or approximately one resonance every 3 G. With higher fields the Feshbach spectrum becomes sparser, with ∂ξ/∂B ≈ 0.1 G −1 , or one resonance per 10 G, at about B ≈ 1000 G. We find this behaviour to be the same for both the "off-resonance" and "on-resonance" cases.  P WD (s; η) = (πs/2) exp(−πs 2 /4), where the probability density of two resonances overlapping is zero. Real systems can be described by the Brody distribution, P B (s) = c η (1 + η)s η exp(−c η s η+1 ), where c η = Γ ((η + 2)/(η + 1)) η+1 . The variable parameter η makes it possible to smoothly interpolate between the idealized Poisson (η = 0) and Wigner-Dyson (η = 1) distributions. Following Green et al. [28] we calculate the Brody parameter η through max-likelihood estimation: η can be found by maximizing the log-likelihood l(η) = i ln P B (s i ; η), while the standard deviation u(η) = ∂ 2 l/∂η 2 −1/2 . In practice, both η and u(η) can be easily found by fitting with a low order polynomial to a scan of l(η) between η = 0 and η = 1. For the "off-resonance" case we find η = 0.71 (11) and η = 0.54 (10) for the "on-resonance" spectrum. These two values of η are statistically consistent and characterize ErYb as a moderately chaotic system.
The number variance, Σ 2 (∆ξ), is defined as the variance of the number of resonances within the interval [ξ i , ξ i + ∆ξ], calculated over all resonances ξ i . For a Poisson distribution of resonance spacings, Σ 2 P (∆ξ) = ∆ξ. For the Wigner-Dyson distribution the resonance spacings are much more even, and the number variance for large ∆ξ can be approximated with Σ 2 WD (∆ξ) = 2π −2 ln(2π∆ξ) + γ + 1 − π 2 /8 , where γ ≈ 0.577 . . . is the Euler constant. Figure 5(f) shows the number variance calculated for the on-and off-resonance cases concerned. We find that Σ 2 for the off-resonance case closely matches the Wigner-Dyson limit indicating strong level repulsion. The on-resonance result deviates slightly from the ideal Wigner-Dyson behaviour, but still indicates a strong signature of quantum chaos.

Role of anisotropy
The mechanism behind magnetic Feshbach resonances in ErYb is the anisotropy of atomic interactions that couples different even partial waves. To gain further insight into how chaos emerges in the Feshbach spectra we introduce an anisotropy scaling parameter λ, replace the anisotropic potentials V k=2,4,... by λV k=2,4,... and repeat our statistical analysis as λ is varied between zero and two. We should also point out that our knowledge of the anisotropy is of pure ab initio origin and thus we also test the robustness of our predictions with respect to small inaccuracies in the total anisotropy. Figures 6(a-f) show the previously discussed properties of the ErYb system as a function of λ for the "off-resonance case".
We first look into zero-field near-threshold bound state energies. At zero anisotropy the vibrational levels have well defined rotational quantum numbers L [solid lines in Figure 6(a)]. As anisotropy is gradually introduced, the initial degeneracy of the M L sublevels is lifted and bound states split into L + 1 components which correspond to states with different |M L | quantum numbers. Interestingly, also quasibound states with L > 0 become bound as the anisotropy emerges. At λ 0.8 signatures of avoided crossings start to emerge and the system becomes manifestly chaotic for λ 1.1. The unscaled ErYb system corresponds to λ = 1 and this corroborates our previous conclusion that ErYb is a 'moderately chaotic' system.
The statistical measures of chaos for our calculated Feshbach spectra also reveal the decisive impact of anisotropy on resonance correlations. Figure 6(b) shows the Brody parameter η. At λ close to zero, the Brody parameter η ≈ 0.1 − 0.3 indicating a near complete lack of level repulsion. With λ ≈ 1 the values of η are close to 0.7, again showing moderate chaos. For larger λ the values of the Brody parameter appear to saturate near η ≈ 0.8 − 0.9, consistent with strong level repulsion. Similarly, the number variances, shown in Figure 6(c) for intervals ∆ξ = 1, 2, 3, all exhibit a downward trend as λ is ramped.
We have also performed calculations of near-threshold bound states with the anisotropy limited to just the V 2 term [dotted lines in Fig. 6(a)]. Interestingly, the pattern of bound states in such case remains very close to one obtained with the full anisotropy (blue solid lines). Thus, it is clear that the role of interactions between closed channels is far more important than the direct coupling of L = 0 scattering state with corresponding outgoing channels via V 4 and higher anisotropies. The dominating role of V 2 potential can be used in future for introducing a model potential that could reproduce the experimental spectrum of resonances using a small number of parameters. A similar situation can likely be expected in interactions of other systems with ultracold lanthanides: the most important information on the interaction could be hidden in the dominant anisotropic term, and the complex structure would originate from the coupling between closed channels.
The anisotropic potential V k mixes many channels with M L quantum number and m j which makes the bound states a composite of m j quantum states. In the absence of anisotropy the differential magnetic moment with respect to the initial state ∆m eff j (as defined in Sec. 4) takes even integer values between 2 and 12 [ Fig. 6(d)]. As the anisotropy is ramped up ∆m eff j becomes non-integer and starts from 0. Surprisingly, the mixing of m j even for a small anisotropy is quite strong: for λ = 0.1 the ∆m eff j bound states start to deviate from integer values Since it is possible to experimentally measure the magnetic moments of the sub-threshold bound states by radio-frequency spectroscopy measurements [76], the statistical distribution of ∆m eff j could be investigated. The experimental results could be later confronted with theoretical predictions to test the anisotropy obtained with ab initio calculations.

Sensitivity of Feshbach spectra to variations of reduced mass
There is a number of pairs of ErYb isotopic combinations that have nearly identical reduced masses. In particular the smallest difference is between 168 Er 170 Yb and 170 Er 168 Yb, where the latter pair has a reduced mass larger by just δµ ≈ 5.6 × 10 −4 u (corresponding to a relative change δµ/µ ≈ 6.6 × 10 −6 ). Given the minuscule change in the reduced mass, and the deterministic nature of the Schrödinger equations, one can expect the Feshbach spectra to be nearly identical for these two isotopic combinations. Indeed, the calculated spectra, shown in Fig. 7(a), are very similar, with resonance shifts much smaller than the spacings between resonances. This makes it possible to unambiguosly match the resonances between the two spectra despite the chaotic nature of the resonance spacing. The resonance shifts, shown in Fig. 7(b) vary roughly as δB res ∼ √ B res and reach at most 0.3 G for resonance positions close to 1000 G. Resonances at larger fields have higher sensitivities to mass variations because they correspond to more deeply bound states, whose positions in turn are more sensitive [77,78]. The large variation in the shifts between adjacent resonances can be attributed to the repulsion between molecular bound states that is central to quantum chaos. We find no obvious correlations between the shifts and resonance widths ∆B, as shown in Fig. 7(c).
Feshbach resonances have been proposed as a means for searching for temporal variations of the proton-to-electron mass ratio [77,79]. In a molecular context, and in absence of hyperfine interactions, the sensitivity to this fundamental constant is mainly due to its direct influence on the reduced mass µ. Our calculated shifts δB res can thus be used to (trivially) evaluate the relative sensitivity of resonance positions to variations in the reduced mass and, by proxy, the proton-to-electron mass ratio. The same could be done in experiment: a comparison of experimental Feshbach   [38,39] and any bosonic Er isotope can in principle be condensed because the Er-Er interaction can be tuned magnetically. In ErYb K µ (B res ) are largest for small fields, where for B < 100 G we can expect many resonances with K µ (B res ) of a few hundred. Given that magnetic fields can be controlled experimentally up to about 10 −6 , the positions of Feshbach resonances could constrain the variations of proton-to-electron mass ratio up to about 10 −8 -10 −9 at best. Much more useful constraints could potentially be achieved by monitoring the scattering length near a Feshbach resonance [77]. Assuming a Breit-Wigner profile [Eq. (8)] the relative sensitivity of the scattering length due to shift in resonance position § is For realistic conditions of B res = 100 G, ∆B = 10 −4 G, K µ (B res ) = 100 and assuming a reasonable detuning from resonance, B − B res = 2∆B, we obtain K µ (a(B)) on the order of 10 8 . If the scattering length can be measured to within 10 −3 -10 −6 then variations of the proton-to-electron mass ratio can be constrained down to 10 −11 -10 −14 .

Conclusion
In conclusion, we have investigated the ultracold collisions of bosonic Er and Yb atoms by means of accurate ab initio calculations. Due to the slightly anisotropic Er atom, the ground state of ErYb is split into eleven asymptotically degenerate states with Λ = 0, ±1, . . . , ±5 with similar well depths, the deepest of which is 639 cm −1 . The potentials are split by at most 40 cm −1 at equilibrium distance. The same small anisotropy (C 6,2 /C 6,0 ≈ 1%) is the mechanism behind the Feshbach resonances: the anisotropy couples different rotational quantum numbers. This coupling is also responsible for the emergence of chaos in the positions of both near-threshold bound states and magnetic Feshbach resonances. The resonance density is approximately 0.3 G −1 for small fields (< 100 G) and decreases to about 0.1 G −1 at 1000 G. The resonance spacings exhibit features of quantum chaos as evidenced by a Brody parameter η = 0.5 − 0.7 and a number variance behaviour consistent with the predictions of the Wigner-Dyson theory. Despite the chaotic behavior we reliably observe several broad resonances for fields < 300 G that may be useful for scattering length control or molecular magnetoassociation. In our calculations using scaled anisotropy we have demonstrated the decisive impact of anisotropy on the emergence of quantum chaos in the Feshbach spectra.
The ErYb system offers 20 bosonic isotopic combinations whose reduced masses range from 82.43 u (for 162 Er 168 Yb) to 85.92 u (for 170 Er 176 Yb) and cover a full cycle of scattering lengths. Uniquely, the pairs 168 Er 170 Yb and 170 Er 168 Yb have reduced masses that match to better than 10 −5 relative. The resonance spectra for these pairs are so similar that it is possible to identify the resonances between them. In experiment this could be used to determine the sensitivity of these resonances to small variations of the proton-to-electron mass ratio with implications for fundamental physics [77][78][79]. Experimentally, our predictions could be tested most cleanly by forming a dual Mott insulator of Yb and Er atoms. Separately, Mott insulators have already been demonstrated for both species. Keeping unity filling for the Er gas could help avoid additional losses due to the much denser Er-Er Feshbach spectrum.
We expect our results for bosonic ErYb to qualitatively hold for other systems of highly magnetic lanthanides (like Dy, Ho, or Tm) combined with spin-singlet atoms (like Sr, Yb, or Hg). Since atoms in the 1 S 0 are spherically symmetric, the anisotropy in a given combination of atomic species will only depend on the selected highly magnetic atom. Some variation between the different combinations may be expected due to the different reduced masses which in turn influence the rotational couplings. A combination of species involving the comparatively light Sr atom could form a mixture with a stronger mass imbalance that could be of interest in context of Efimov physics.
Lastly, ErYb offers a number of fermion-fermion and fermion-boson combinations which may be of interest for the quantum simulation and many-body physics communities. The behaviour of many-fermion ensembles is central to most important areas of physics, and in vastly different contexts: whether the design of new superconducting materials or understanding the structure of neutron stars. Finding ways to design new types of fermionic mixtures with various types of interactions (dipolar for Er and isotropic for Yb) can help build new types of experimental setups ready for new challenges in many-body quantum physics. Before the advent of quantum gases composed of lanthanides the research on fermionic mixtures was limited only to alkali quantum gases: potassium-40 and lithium-6, or their mixtures [80][81][82]. The Er+Yb system can make 2 Fermi-Fermi mixtures: 167 Er combined with either the 171 Yb or 173 Yb isotope. The 167 Er ground state has f = j + s + i = 19/2. The first excited hyperfine state f = 17/2 is separated by a 2.7 GHz gap [83]. Thus, the first crossing of the two lowest thresholds of 167 Er+ 171 Yb and 167 Er+ 173 Yb systems alike corresponds to a magnetic field of 2.7GHz/g j µ B =1953 G. Close to the crossing the density of bound states supported by f = 17/2 can be expected to be large leading to a densification of Feshbach resonance spectra for fields slightly below this value. An additional difference from the bosonic spectra is the possible Zeeman splitting of Feshbach resonances into multiplets due to the nuclear magnetic moment of fermionic Yb, similarly as in RbSr or CsYb molecules [84,85]. The detailed analysis of interactions in fermion-fermion mixture, however, is beyond the scope of this paper, and we will study it in the near future.