Universality of Efimov states in highly mass-imbalanced cold-atom mixtures with van der Waals and dipole interactions

We consider a three-body system of two heavy particles (mass $M$) which interact resonantly with a light particle (mass $m$). The heavy particles are assumed to be either identical spin-polarized bosons or fermions, to model ^166,168,170Er and ¹⁶⁷Er atoms, respectively. As the light particle and heavy particles interact with a large $s$-wave scattering length $a_{\mathrm{HL}}$, a condition for the Efimov states to appear, we model the light-heavy interaction by the zero-range contact interaction [48]. On the other hand, the interaction between the heavy particles is assumed to be the sum of van der Waals and dipole interactions. This amounts to neglecting the van der Waals interaction length scale between the light and heavy particles, which should be a reasonable assumption because the van der Waals length between the heavy atoms is larger than that between the light-heavy atoms. Throughout this paper, we focus on the unitary limit $1/a_{\mathrm{HL}}=0$, and model the light-heavy interaction via the unitary Bethe-Peierls boundary condition [48]. We note that this treatment of the inter-species interaction is valid in describing broad Feshbach resonances, but less so for narrower ones [10]

In a highly mass-imbalanced system, the Born-Oppenheimer approximation provides a good description of the Efimov physics. Indeed, it successfully reproduces the Efimov scale factor, and the critical mass ratio for the Efimov effect $M/m=13.99\ldots$ [2] in good agreement with the exact value $M/m=13.60\ldots$ [35]. The Born-Oppenheimer approximation has also been confirmed to reproduce well the van der Waals universality of the Efimov states without the dipole interaction for mass ratios $14\lesssim M/m\lesssim 29$ [24]. It is therefore a good approximation for an Er-Li mixture, whose mass ratio is $M/m\approx 28$. The light particle’s Schrödinger equation is solved analytically and leads to the induced attractive interaction between the heavy particles $-\Omega^{2}/2mr^{2}$, where $\Omega=0.5671\ldots$ [3]. The Schrödinger equation between the heavy particles then reads

where $\theta$ is an angle measured from the $z$ axis taken to be parallel to the orientation of the dipoles (parallel to the external magnetic field). Here, $C_{6}$ and $C_{dd}$ are the van der Waals and dipole coefficients, from which we can define the van der Waals length $r_{\mathrm{vdw}}=\dfrac{1}{2}(MC_{6})^{1/4}$ and the dipole length $a_{dd}=MC_{dd}/3$ between the heavy atoms.

The binding energy of the Efimov states (i.e., three-body parameter) can be obtained by solving Eq. (1). It needs to be supplemented with a short-range boundary condition. For this purpose, we use a hard wall condition at $R_{\mathrm{min}}$: $\psi(r=R_{\mathrm{min}})=0$. We will see later that the results are insensitive to the way the boundary condition is introduced as long as $R_{\mathrm{min}}$ is small enough. As the dipole interaction couples different partial waves, the angular momentum is no longer a good quantum number, but the parity $\Pi$ and the $z$ component of the angular momentum $M_{z}$ are good quantum numbers. While various values of $M_{z}$ are allowed in general, we restrict ourselves to the $M_{z}=0$ channel in this paper. For the bosonic heavy atoms, the $M_{z}^{\Pi}=0^{+}$ state is allowed by the symmetry. This state involves the $L=0,2,4,...$ angular momentum channels, among which the $L^{\Pi}=0^{+}$ channel shows the Efimov effect and the others do not for the Er-Li mass ratio [2, 36, 49, 50, 3]. For the fermionic heavy atoms, the $M_{z}^{\Pi}=0^{-}$ channel is allowed by the symmetry. This state involves $L=1,3,5,...$ angular momentum channels, among which the $L^{\Pi}=1^{-}$ channel shows the Efimov effect and the others do not for the mass ratio in question [2, 36, 49, 50, 3]. Therefore, in the following, we study the $M_{z}^{\Pi}=0^{+}$ and $M_{z}^{\Pi}=0^{-}$ states for the bosonic and fermionic Er atoms, respectively.

For typical cold atomic species, $r_{\mathrm{vdw}}\gg a_{dd}$ so that the dipole interaction is negligible in Eq. (1). While this is not the case for Er atoms, it still provides us with a basis for understanding the universal behavior of the Efimov states with the dipole interactions. Therefore, in this section we analyze

with a short-range boundary condition $\psi(|\bm{r}|=R_{\mathrm{min}})=0$. While this equation was studied numerically in Ref. [25] and the van der Waals universality of the Efimov states was clarified for the mass-imbalanced system, we show below that Eq. (2) can be solved analytically and provide explicit formulas for the binding energies and wave functions of the Efimov states.

Without the dipole interaction, the angular momentum is a good quantum number, and Eq. (2) is essentially a single-channel Schrödinger equation with $1/r^{6}$ potential:

with $\psi(\bm{r})=\dfrac{u_{\ell}(r)}{r}Y_{\ell m}(\bm{r})$ and

Since we are interested in the Efimov states, we focus on the $\ell=0$ ($\ell=1$) channel where $s_{\ell}^{2}<0$ for the mass ratio of interest $M/m\approx 28$ for bosons (fermions), respectively [2, 36, 3, 49, 50]. Equation (3) can be solved analytically using Ref. [51], with a remark that the centrifugal term is attractive due to the Efimov attraction. In other words, the angular momentum in Ref. [51] should be generalized to an imaginary number. We can then apply the quantum defect theory (QDT) [52, 53] to obtain the analytical behavior of the Efimov states; using the QDT parameter $K^{c}$ introduced as

where $f_{\ell}^{c}$ and $g_{\ell}^{c}$ are two independent solutions of Eq. (3). The binding energies of the Efimov states can be obtained by imposing $u_{\ell}(r\rightarrow\infty)=0$, from which we find [52]

The parameters $\theta_{\ell}$, $\nu$, and $M_{\ell}$ can be obtained from the standard procedure of solving the van der Waals two-body problem (see Appendix. A). At low energy $|E|\ll 1/Mr^{2}_{\mathrm{vdw}}$, by performing the low-energy expansion in a similar manner to Refs. [53, 54], we find (see Appendix. B )

where $\xi_{\ell}=\arg\left[\Gamma\left(\frac{i|s_{\ell}|}{2}\right)\Gamma(1+i|s_{\ell}|)\right]$ and $n$ is an arbitrary integer satisfying $|E|\ll 1/Mr^{2}_{\mathrm{vdw}}$.

The QDT parameter $K^{c}$ is formally a matrix acting on multiple channels, but it simplifies as a scalar parameter for the single-channel scattering. Its value is related to the short-range phase, hence the short-range boundary condition. Indeed, at short distance, one can substitute the low-energy form of the wave function Eq. (28) because the energy is much smaller than the potential, and imposing $u_{\ell}(r=R_{\mathrm{min}})=0$, one finds

where $\Phi=2r_{\rm vdw}^{2}/R^{2}_{\mathrm{min}}$. In the second line, we have used the asymptotic form of the Bessel function (or equivalently Eq. (25)) valid for $R_{\mathrm{min}}\ll r_{\rm vdw}$. This relation between $R_{\mathrm{min}}$ and $K^{c}$, substituted into Eq. (7) gives an explicit analytical formula of the Efimov binding energies as a function of $R_{\mathrm{min}}$. The wave function of the Efimov states is also analytically obtained as ($x\equiv 2\left(r_{\rm vdw}/r\right)^{2}$)

where the log-periodic behavior characteristic of the Efimov states appears from the oscillation of the imaginary-index Bessel function.

While the $R_{\mathrm{min}}$ dependence of the binding energy and wave function may look like a non-universal short-range dependence of the Efimov states, we can recast it into universal forms by using the relation between $R_{\mathrm{min}}$ and the two-body universal scattering parameter between the heavy atoms; if we consider the two-body problem of the two heavy atoms

the only difference with Eq. (2) is the presence/absence of the $1/r^{2}$ attraction. At short distance $r\sim R_{\mathrm{min}}\ll r_{\mathrm{vdw}}$, this attraction is much smaller than the van der Waals force so that the two equations are essentially the same. More physically, the short-range phase of the three-body system is dominated by the potential between the heavy atoms. This assumption should be particularly true for a highly mass-imbalanced system, because due to larger atomic polarizability of the heavy atoms, their van der Waals force and higher-order forces are much stronger than those for the light atom. Therefore, the QDT parameter of the three-body system can be approximated by that of two heavy atoms $K^{c}_{\mathrm{3body}}\simeq K^{c}_{\mathrm{2body}}$, which amounts to adopting the same short-range cutoff $R_{\mathrm{min}}$ value for the two- and three-body systems. While this assumption is plausible for our model, it can be compromised by finite-range effects of the heavy-light atoms’ interaction. It may also be invalidated by the non-adiabatic effects beyond the Born-Oppenheimer approximation, which break down the single-channel description.

The van der Waals two-body problem Eq. (12) with the QDT condition Eq. (5) has been solved in Refs. [53, 54], and the following relations between the the QDT parameter and the two-body scattering parameters have been obtained:

where $a^{(\mathrm{HH})}$ and $v_{p}^{(\mathrm{HH})}$ are the $s$-wave scattering length and $p$-wave scattering volume between the heavy atoms, respectively. Noting that the bosonic (fermionic) heavy atoms interact dominantly with the $s$-wave scattering length ($p$-wave scattering volume) at low energy, and substituting Eqs. (13) and (14) into Eq. (7), the binding energies of the Efimov states are represented as

for the bosonic heavy atoms ($\ell=0$), and

for the fermionic heavy atoms ($\ell=1$), where $n$ is an arbitrary integer such that $|E|\ll 1/Mr_{\rm vdw}^{2}$.

As we will see in Sec. III.1, these analytical formulas excellently reproduce the numerical results for non-dipolar systems. These analytical formulas suggest that the binding energy of the Efimov states in the unitary limit, i.e. the three-body parameter, is insensitive to the short-range details of the system but universally determined by the van der Waals length and the $s$-wave scattering length ($p$-wave scattering volume) for the bosonic (fermionic) Efimov states. This van der Waals universality has been experimentally demonstrated for equal-mass [19, 20] and mass-imbalanced bosonic systems [29] close to a broad Feshbach resonance. The van der Waals universality has also been demonstrated theoretically with numerical calculations for equal-mass [23, 21, 22] and mass-imbalanced bosonic systems [24, 25]. In particular, the same Born-Oppenheimer equation (2) with the same short-range boundary condition has been employed in Ref. [25] and solved numerically to study the universality of the Efimov states. Our QDT analysis provides an analytical demonstration of the van der Waals universality of hetero-nuclear Efimov states; namely, the assumption of $K^{c}_{\mathrm{3body}}\simeq K^{c}_{\mathrm{2body}}$ is pivotal in universally characterizing the three-body parameter with the two-body scattering parameters between the heavy atoms. As we use the single-channel QDT, Eqs. (15) and (16) are invalid for describing the Efimov states near a narrow Feshbach resonance, which have been demonstrated both experimentally [26, 27, 28, 29] and theoretically [30, 31, 32, 33, 34] to deviate from the van der Waals universality; the multi-channel QDT [55, 56] is necessary for capturing their sensitivity to the channel coupling and nature of the internal spin states.

First, we show the results without the dipole interaction: in Fig. 1, we show the binding energy of the Efimov states (in the van der Waals energy unit $1/Mr_{\mathrm{vdw}}^{2}$) as the short-range boundary $R_{\mathrm{min}}$ is varied. As $R_{\mathrm{min}}$ gets smaller, the short-range van der Waals attraction increases, leading to the successive appearance of bound states and increase of the binding energy. The numerical data points almost perfectly agree with the analytical curves obtained with the QDT in Eq. (7), particularly at small energy. There are a few points which significantly deviate from the analytical curves. They are ascribed to higher angular momentum states of $\ell=2,4,...$ for the bosons and $\ell=3,5,...$ for the fermions, respectively. The thresholds at which the $\ell=2$ ($\ell=3$) bound states start to appear, denoted by the vertical dotted lines in Fig. 1 (a) and Fig. 1 (b), agree with most of the stray data points, supporting this interpretation. A few stray points which agree with neither are ascribed to higher angular momentum channels $\ell\geq 4$ ($\ell\geq 5$). This classification of the states has also been confirmed from the shape of the wave functions. Because the higher angular momentum states are localized to shorter distances than the $s$-wave ($p$-wave) state owing to the centrifugal barrier, they are more sensitive to $R_{\mathrm{min}}$, resulting in the steeper binding energy curves. The binding energy curves of different partial waves do not render an avoided crossing but simply cross because the angular momentum is a good quantum number in the absence of the dipole interaction.

In Fig. 2, we show the binding energies plotted against the $s$-wave scattering length $a^{(\mathrm{HH})}$ and $p$-wave scattering volume $v^{(\mathrm{HH})}$ between the heavy atoms, obtained from Fig. 1 and Eqs. (8) (13) (14). The numerical data with different values of $R_{\mathrm{min}}$ collapse into universal curves. Therefore, the three-body parameter of the Efimov states is universally characterized by ($a^{(\mathrm{HH})}$, $r_{\mathrm{vdw}}$) for the bosons and ($v_{p}^{(\mathrm{HH})}$, $r_{\mathrm{vdw}}$) for the fermions, respectively. In other words, they are insensitive to the short-range details, but are universally characterized by the two-body scattering parameters which describe the long-range asymptotic behavior of the two heavy atoms. The numerical data in Fig. 2 excellently agree with the curves of Eqs. (15) and (16), vindicating our QDT analytical formulas. Some points in Fig. 2 deviating from the curves are higher angular momentum states, whose energies change rapidly as $R_{\mathrm{min}}$ hence $a^{(\mathrm{HH})}$ and $v_{p}^{(\mathrm{HH})}$ are varied; although they significantly deviate from the analytical curves of Eqs. (15) and (16), the data points with different $R_{\mathrm{min}}$ still tend to collapse into a single curve (Fig. 2 (b)(c)) due to the angular-momentum-insensitive nature of the short-range phase originating from $R_{\mathrm{min}}\ll r_{\mathrm{vdw}}$ [52]. The marginal disagreement of the large binding energy region $|E|\simeq 1/Mr_{\mathrm{vdw}}^{2}$ in Fig. 2 (c) is ascribed to the breakdown of the low-energy condition used to derive Eqs. (15) and (16).

We note that the universal behaviors of the bosonic system in Fig. 2 (a)(b) were found numerically in Ref. [25]; we have found here that a similar universal behavior holds true for the fermions (Fig. 2 (c)(d)), together with their analytical descriptions in Eqs. (15) (16).

In Fig. 2, we also show the dimer energy of two heavy atoms (green dashed) obtained by solving Eq. (12). Wherever the Efimov trimers lie above the dimer energy, they are in fact resonant trimers embedded in the dimer-atom continuum [25]. In the Born-Oppenheimer approximation, however, the trimers appear as bound states because the channel describing the dimer-atom continuum is neglected.

We show in Fig. 3 the binding energies of the Efimov states in the presence of a weak dipole interaction. As the angular momentum is no longer a good quantum number, different partial wave states mix with each other, resulting in avoided crossings. We have confirmed that the shape of the wave function changes across the avoided crossing via a superposition of different partial waves. We have also confirmed that the width of the avoided crossings increases as $a_{dd}$ gets larger.

As shown in Fig. 3, the numerical results generally get shifted toward larger binding energy. This can be understood by the perturbation theory: we have analytically found that the energy shift induced by the dipole interaction is $E-E^{(0)}\propto a_{dd}^{2}$ for the bosons predominantly in the $L=0$ channel, and $E-E^{(0)}\propto a_{dd}$ for the fermions predominantly in the $L=1$ channel, with both negative shifts $E-E^{(0)}<0$ (see Appendix. C). The difference in the scaling originates from the absence (presence) of the diagonal matrix element of the dipole interaction for $L=0$ ($L=1$) states: $\langle L=0|V_{dd}|L=0\rangle=0$, $\langle L=1|V_{dd}|L=1\rangle\neq 0$. This is confirmed in Fig. 4: the change in binding energy, i.e., three-body parameter of the Efimov states, scales quadratically for the bosons, and linearly for the fermions. Notably, in most previous studies on the Efimov effect, the role of the statistics of the particles has been rather marginal, at best to modify the value of $s_{\ell}$, and thereby change the absence/presence of the Efimov effect through the condition $s_{\ell}^{2}<0$ as the mass ratio is varied [2, 35, 49, 50, 3, 4, 5, 6]. In other words, the statistics of the particles has been essentially irrelevant once its effect is incorporated into the value of $s_{\ell}$. Figure 4 is one of the few examples where the quantum statistics of the particles plays an explicit role in qualitatively changing the universal behavior of the Efimov states.

Erbium atoms interact with a strong dipole interaction owing to their large magnetic moments: the dipole length of ¹⁶⁶Er-¹⁶⁶Er is $a_{dd}=75.5a_{0}$, which is comparable to their van der Waals length $r_{\mathrm{vdw}}=65.5a_{0}$ [57, 58] (see Table 1). Figure 4 shows that this is away from the region where the perturbation theory works well $a_{dd}/r_{\mathrm{vdw}}\lesssim 0.5$.

We have therefore performed the numerical calculations for moderate strengths of the dipole interaction in Fig. 5(a)-(f); we first show the binding energy of the bosonic Efimov states plotted against the heavy-heavy $s$-wave scattering length estimated without the dipole interaction $a_{\mathrm{no-dd}}^{(\mathrm{HH})}$ (Eq. (8)). As the dipole strength is increased from a smaller ((a)(b)) to a larger ((e)(f)) value, the avoided crossing between the $\ell=0$ dominant state and the $\ell=2$ dominant state gets broader. Due to the admixture of the $\ell=0$ state, the slope of the $\ell=2$ dominant state appearing around $a^{(\mathrm{HH})}_{\mathrm{no-dd}}/r_{\mathrm{vdw}}\simeq 1$ gets gentler. Furthermore, the binding energies increase as the dipole interaction gets larger, which is consistent with the perturbative analysis in Sec. III.2. The deviation from the analytical curves of the non-dipolar system in Eq. (15) (black curves) is so significant that they cannot well explain the numerical results of the realistic dipole interaction strength (middle row (c)(d)), let alone stronger ones (bottom row (e)(f)). Still, the binding energy data of different short-range parameters $R_{\mathrm{min}}$ collapse onto a single curve, suggesting the universality of the three-body parameter of the bosonic Efimov states even in the presence of non-perturbatively large dipole interactions.

To further elucidate the nature of this universality, we convert the horizontal axis into an authentic $s$-wave scattering length between the heavy atoms $a^{(\mathrm{HH})}_{\mathrm{with-dd}}$ in Fig. 5(g)-(l); in the presence of the dipole interaction, the $s$-wave scattering length between the heavy atoms $a^{(\mathrm{HH})}_{\mathrm{with-dd}}$ gets significantly modified from a purely van der Waals value of $a^{(\mathrm{HH})}_{\mathrm{no-dd}}$ in Eq. (13) [59, 60]. As the experimentally accessible scattering length between the Er atoms is $a^{(\mathrm{HH})}_{\mathrm{with-dd}}$ rather than $a^{(\mathrm{HH})}_{\mathrm{no-dd}}$, it should be more physical to use $a^{(\mathrm{HH})}_{\mathrm{with-dd}}$ to eliminate $R_{\mathrm{min}}$. In Fig. 6, we show $a^{(\mathrm{HH})}_{\mathrm{with-dd}}$ (circles) obtained by numerically solving the coupled-channel Schrödinger equation between the two heavy atoms

with the boundary condition $\psi(|\bm{r}|=R_{\mathrm{min}})=0$, and compare it with $a^{(\mathrm{HH})}_{\mathrm{no-dd}}$. While they are almost indistinguishable for a weak dipole interaction (Fig. 6(a)) except for the proximity of higher-partial wave induced resonances (dotted vertical lines) [59, 61], $a^{(\mathrm{HH})}_{\mathrm{with-dd}}$ is markedly different from $a^{(\mathrm{HH})}_{\mathrm{no-dd}}$ for a realistic dipole strength (Fig. 6(b)), and stronger one (Fig. 6(c)). This is expected because the realistic dipole strength of the Er atoms $a_{dd}/r_{\mathrm{vdw}}=0.86755....$ in Fig. 6(b) (and stronger one in Fig. 6(c)) is away from the perturbative regime which was confirmed to be $a_{dd}/r_{\mathrm{vdw}}\lesssim 0.5$ in Sec. III.2.

In Fig. 5(g)-(l), we show the binding energies of the Efimov states plotted against $a^{(\mathrm{HH})}_{\mathrm{with-dd}}$. Similar to Fig. 5(a)-(f), the data with different values of $R_{\mathrm{min}}$ mostly collapse into a single curve, demonstrating the universality. One notable difference is that the agreement between the numerical data and the purely van der Waals analytical curves (black) is much better in (g)-(l) than those in (a)-(f): except for the broad avoided crossing at $a_{\mathrm{with-dd}}^{(\mathrm{HH})}/r_{\mathrm{vdw}}\simeq 1$, a significant fraction of the numerical results agree well with the van der Waals curves not only for the weak dipole interaction, but also for the realistic dipole strength (Fig. 5(i)(j)) and stronger dipole interaction ((k)(l)). This suggests that the van der Waals QDT of the bosonic Efimov states in Eq. (15) can capture the universal behavior of the Efimov states in the presence of the dipole interaction, once the effects of the dipole interaction are incorporated into the $s$-wave scattering length between the heavy atoms: $a^{(\mathrm{HH})}=a^{(\mathrm{HH})}_{\mathrm{with-dd}}$. In other words, the major role of the dipole interaction is to change, hence renormalize, the value of the $s$-wave scattering length between the heavy atoms, so that the van der Waals universality of the Efimov states in Eq. (15) persists with the renormalized $s$-wave scattering length $a^{(\mathrm{HH})}_{\mathrm{with-dd}}$. This renormalized van der Waals universality for the dipolar bosonic Efimov states we have found here seems non-pertubative; the dipole interaction is so large that it significantly changes both the binding energies (Fig. 5(c)-(f)) and the $s$-wave scattering length (Fig. 6(b)(c))) from the purely van der Waals system. Some horizontal data points in Fig. 5(g)(l) which are not present in (a)-(f) are artifacts of the higher-partial-wave-induced resonances: the two-body problem between the heavy atoms in Eq. (17) not only shows the $s$-wave dominated resonances, but also shows the resonant enhancement of $a^{(\mathrm{HH})}$ induced by higher-partial-wave channels [59, 61], whose approximate positions are indicated by the vertical black dotted lines in Fig. 6. The narrow nature of these higher-partial-wave-induced resonances leads to an abrupt change of $a^{(\mathrm{HH})}_{\mathrm{with-dd}}$ as $R_{\mathrm{min}}$ is changed, resulting in flat energy spectra in Fig. 5(g)-(l). Equation (15) cannot explain these flat data points. We note, however, that these flat energy spectra are absent in Fig. 5(a)-(b), and they should therefore be interpreted as artifacts originating from the two-body effect of $a^{(\mathrm{HH})}$ rather than a genuine three-body phenomenon. Physically, when we attempt to evaluate the binding energies of the Efimov states from an experimentally obtained value of $a^{(\mathrm{HH})}$ (see Sec. IV), we can simply use the van der Waals curves of Eq. (15) by neglecting the stray flat spectra, except for the case where the two Er atoms are coincidentally in the proximity of the $s$-wave resonance induced by a $d$-wave or higher-angular-momentum molecular state. In the avoided crossing region, the variance of the spectrum is larger than in the other regions (see Fig. 5(l)). This is due to the finite values of $R_{\mathrm{min}}$, as it can be checked that reducing it further does lead to a converged universal behavior. For the finite values of $R_{\mathrm{min}}$ used in the figures (which is varied between 0.25 and 0.4 $r_{\mathrm{vdw}}$), it is still justified to call this region’s behavior universal despite the noticeable energy variance because it remains much smaller than the Efimov period, in sheer contrast with the strength of the van der Waals interaction at short distance, which varies by a factor 16. Yet, a $d$-wave dimer of the two heavy atoms appearing around $a_{\mathrm{with-dd}}^{(\mathrm{HH})}/r_{\mathrm{vdw}}\simeq 1$ (see Fig. 11 in Appendix. D where the dimer energies are presented) seems to significantly affect two-body and three-body behaviors in this region, making the simple renormalized van der Waals universality of Eq. (15) inadequate.

While the physical origin of the renormalized van der Waals universality is yet to be clarified, we suspect that the angular-momentum-insensitive nature of the short-range phase [52] is playing a crucial role; we have confirmed in our three-body calculations that a substantial fraction of the wave function is occupying $\ell\geq 2$ states. If the short-range phase shifts in various angular momentum channels are directly related with each other, as in the angular-momentun-insensitive QDT [52], $a^{(\mathrm{HH})}$ characterizing the $s$-wave short-range phase can universally capture the scatterings occuring in all the angular momentum channels. We also suspect that a universal three-body repulsion appearing at moderately large distance may be another possible mechanism [23, 21, 22, 46, 47]; in the three-body system of identical bosons interacting via the van der Waals interaction, the appearance of a non-adiabatic three-body repulsion at a large distance $r\simeq 2r_{\mathrm{vdw}}$ is ascribed to be the physical origin of the universality of the three-body parameter of the Efimov states [23, 21, 22]. Furthermore, for three identical bosons interacting via purely dipolar interactions, a three-body repulsion has also been found to appear universally at $r\simeq a_{dd}$ [46]. The multi-channel effects with an interplay of the van der Waals and dipole interaction in our system may lead to a similar three-body repulsion scenario. We also remark that our universality of the Efimov states involving higher partial-wave channels is analogous to the partial-wave phase-locking mechanism found in the resonant exchange collisions [62, 63], where the phase shifts of different scattering channels between the atoms are related with each other.

Since the $s$-wave scattering length between the heavy atoms characterizes the universal behavior of the bosonic Efimov states so precisely, it is tempting to use the $p$-wave scattering volume $v_{p}^{(\mathrm{HH})}$ to analyze the fermionic Efimov states. However, $v_{p}^{(\mathrm{HH})}$ cannot be defined in the presence of the dipole interaction, because the dipole interaction is present even at large distance in the $\ell=1$ channel and the asymptotic wave function is no longer a free plane wave, in stark contrast to the $\ell=0$ channel [60]. We therefore consider below an alternative scattering parameter $K_{c,F}^{(3D)}$ which characterizes the large-distance asymptotic wave function of the two heavy fermionic dipoles: we consider the Schrödinger equation between the two heavy atoms

We have only kept the van der Waals interaction and the diagonal $\ell=1$ element of the dipole interaction, which are reasonably assumed to play major roles in determining the short- and long-range behaviors. At large distance, the van der Waals interaction is negligible in Eq. (18), so that the wave function behaves at vanishing energy as [64]

where $J_{3}$ and $Y_{3}$ are the Bessel and Neumann functions corresponding to two independent solutions under the central $1/r^{3}$ potential. $K_{c,F}^{(3D)}$ is introduced in an analogous manner as the QDT parameter $K^{c}$ in the van der Waals case (see Eq. (5)), and hence is expected to be a universal parameter characterizing the short-range two-body phase shift. In Fig. 7, we show $K_{c,F}^{(3D)}$ numerically calculated from Eqs. (18) and (19) with the hard-wall boundary condition $u_{1}(r=R_{\mathrm{min}})=0$. We note that the van der Waals term is essential in determining the short-range phase: $K_{c,F}^{(3D)}$ is markedly different from that obtained by neglecting the van der Waals term (black solid), suggesting both the van der Waals and dipole interactions are relevant for characterizing the scattering between the two heavy fermionic atoms.

In Fig. 8, we show the binding energies of the fermionic Efimov states plotted against $K_{c,F}^{(3D)}$ as obtained above. The numerical data of different short-range $R_{\mathrm{min}}$ values collapse into a single curve, demonstrating the universality of the fermionic Efimov states. The universality holds particularly well for a realistic $a_{dd}$ value of ¹⁶⁷Er-¹⁶⁷Er (Fig. 8(a)(b)), where only a few points of higher angular momentum character deviate from the universal curve. In contrast to the bosonic case, we cannot find an analytical expression which could explain the universal energy spectra owing to the impossibility of using $v_{p}^{(\mathrm{HH})}$ and hence Eq. (16). The avoided crossing at $K_{c,F}^{(3D)}\simeq-1$ originates from the $\ell=3$ channel which is present without the dipole interaction at around $r_{\mathrm{vdw}}^{3}/v_{p}^{(\mathrm{HH})}\simeq-0.3$ (see Fig. 2(c) and Fig. 12 in Appendix D). It shifts toward a larger negative $K_{c,F}^{(3D)}$ side with a decreased slope as the dipole strength is increased from a realistic value of ¹⁶⁷Er-¹⁶⁷Er (Fig. 8(a)(b)) toward larger values (Fig. 8(c)-(f)). The foothill of the avoided crossing shows a larger variance for stronger dipole interactions (Fig. 8(c)-(f)). This is likely be a precursor of the breakdown of our assumption that the single-channel short-range phase in $\ell=1$ universally characterizes the three-body system (see Sec. III.4), especially for a very strong dipole interaction and large $K_{c,F}^{(3D)}$ (Fig. 8(c)(e)(f)).

Although the $K_{c,F}^{(3D)}$ parameter has the advantage of capturing the universal behavior of the fermionic Efimov states, it is marginally useful in predicting the three-body parameter of experimental systems of interest, e.g., ¹⁶⁷Er-¹⁶⁷Er-⁶Li. This is primarily because it is not easy, though not impossible, to determine the value of $K_{c,F}^{(3D)}$; it sensitively depends on the position (i.e. $R_{\mathrm{min}}$) and shape of the short-range repulsive potential, and it is therefore challenging to predict its value by theoretical calculations [57, 65, 66]. Furthermore, as the $K_{c,F}^{(3D)}$ dependence appears as the next-next-next leading order if we perform the asymptotic expansion of the $J_{3}$ and $Y_{3}$ at large distance in Eq. (19), it contributes to the low-energy scattering cross section as a small correction term modifying the dominant quasi-universal contribution [60].

While the universal descriptions via partial-wave scattering parameters of $a^{(\mathrm{HH})}$ (boson) and $K_{c,F}^{(3D)}$ (fermion) are shown in Sec. III.3 to be valid for moderate dipole strengths including the realistic values for Er atoms, they should break down in the limit of strong dipole interactions. This is because the single partial-wave description in defining the short-range two-body phase between the two heavy atoms becomes invalid since all the even (odd) angular momentum channels almost equally contribute for the bosons (fermions).

We therefore need to introduce an alternative scattering parameter which should characterize the short-range scattering phase of two strong dipoles. As the two dipoles prefer to be aligned in a head-to-tail, hence quasi-one-dimensional, configuration in the strong dipole limit, we consider a one-dimensional equation

which corresponds to Eq. (17) constrained to the one-dimensional configuration on the $z$-axis. We introduce the scattering parameter $K_{c}^{(1D)}$ between the two dipoles inspired from the asymptotic behavior of the zero-energy wave function of Eq. (20) at $|z|\rightarrow\infty$ [67]

where $J_{1}$ and $Y_{1}$ terms are two independent solutions of Eq. (20) without the van der Waals interaction. We impose a hard-wall boundary condition $u(z=\pm R_{\mathrm{min}})=0$ corresponding to the three-dimensional counterpart. Since the positive and negative $z$ regions are separated, the bosonic and fermionic systems become equivalent, with the same value of $K_{c}^{(1D)}$; their only difference appears in the overall parity, i.e. $P=0$ for bosons and $P=1$ for fermions in Eq. (21).

In Fig. 9, the binding energies of the Efimov states obtained from Eq. (17) are plotted as a function of $K_{c}^{(1D)}$ in Eqs. (20) and (21). The bosonic (blue circles) and fermionic (orange squares) systems, hypothetically assumed to be the same mass ratio, show rather similar behaviors, in stark contrast to the weak and moderate dipole systems of Secs. III.1-III.3. This is because the parity of the wave function plays a marginal role when the two heavy atoms are mostly aligned in a quasi-one-dimensional configuration along the $z$-axis, and the $z>0$ and $z<0$ regions are separated by the hard-wall boundary condition. This is particularly true when the size of the trimer is much smaller than $a_{dd}$, beyond which the centrifugal repulsion in Eq. (17) representing the three-dimensional nature becomes more dominant than the dipole term. For the tightly bound states in Fig. 9, the bosonic and fermionic systems become almost indistinguishable, supporting this scenario. In other words, the statistics of the particles is irrelevant for the strongly dipolar Efimov states as long as their binding energy is larger than the dipole scale $1/Ma_{dd}^{2}$.

The data points with a small slope in Fig. 9, including the weakly bound states, tend to collapse into a single curve for both the bosonic and fermionic systems. Although the spreading of the data points is larger than those in Figs. 5 and 8, it remains much smaller than the Efimov period, even with an almost 16 times change in the short-range part of the van der Waals interaction. We have also confirmed that the spreading gets smaller as $R_{\mathrm{min}}$ is decreased, in a similar manner to Fig. 5. We can thus conclude that $K_{c}^{(1D)}$ universally characterizes the short-range phase, hence the Efimov states with the strong dipole interactions. On the other hand, the data points with steeper slopes do not seem to show the universality with respect to $K_{c}^{(1D)}$. Since those states have significant contributions from large angular momentum states $\ell\gg 1$ and the angular insensitive nature of their short-range phase can break down, it is not unexpected to find that their features cannot be universally captured by the one-dimensional scattering model in Eqs. (20) and (21) where the angular momentum is absent (see also dimers in Fig. 13 in Appendix D with similar behaviors). Nevertheless, the physical mechanism behind the universal and non-universal behaviors in Sec. III.3 and Fig. 9 and how they change between them have yet to be clarified.

In Fig. 10, we show the amplitude of the wave function. For a bosonic Efimov state with a small binding energy (Fig. 10 (a)), it essentially comprises two regions; the inner region $R\lesssim a_{dd}$ where the dipole interaction is so strong that the two heavy particles are preferentially aligned on the $z$ axis, and the outer region $R\gtrsim a_{dd}$ where the dipole interaction is smaller than the kinetic term, thus forming a three-dimensional halo state of a predominantly $s$-wave nature. They are neither a perfect one-dimensional geometry nor a spherical $s$-wave halo owing to the strong admixture of angular momentum states, but their universal nature is essentially determined by the phase in the short-range region, where the one-dimensional configuration is dominant and hence $K_{c}^{(1D)}$ is a universal parameter. For fermionic Efimov states with a similarly small value of the binding energy (Fig. 10 (b)), the short-range and long-range natures are almost the same, with the only notable difference being the presence of a node in the wave function at $z=0$, representing the anti-symmetrization. In contrast, the wave function in the non-universal regime (Fig. 10 (c)) shows a distinct feature: the wave function in the inner region $R\lesssim a_{dd}$ no longer shows a simple quasi-one-dimensional character, but rather a strongly anisotropic shape, which is likely induced by an admixture of various angular momentum states. This is consistent with the breakdown of the universality with respect to $K_{c}^{(1D)}$. When the binding energy is so large that the size of the wave function is smaller than $a_{dd}$ (Fig. 10 (d)), the whole wave function is localized in the inner region $R\lesssim a_{dd}$, forming an almost perfect quasi-one-dimensional state. In this limit, the hard-wall boundary condition $R=R_{\mathrm{min}}$ necessarily produces a node in the wave function at $z\approx 0$, so that the bososnic and fermionic wave functions become the same except for the sign change. In this limit, the bosonic and fermionic systems show the same universal behaviors. While this limit of boson-fermion duality, analogous to the Tonks-Girardeau limit [68, 69], is fascinating, it is challenging to achieve because it requires a very large dipole interaction strength, which is impossible for cold atoms with magnetic dipole interactions, but may be possible with polar molecules [70, 71] or Rydberg atoms [72, 73].