Frustrated orbital Feshbach resonances in a Fermi gas

TL;DR

This paper explores how frustration in orbital Feshbach resonances affects interaction tuning in fermionic gases, revealing suppressed binding energies and weaker pairing, with implications for superconductivity and superfluidity.

Contribution

It provides a theoretical analysis of frustration effects in orbital Feshbach resonances in a $^{173}$Yb system, highlighting their impact on two-body and many-body interactions.

Findings

Suppressed ground state binding energy due to closed-channel atom addition.

Weaker fermion pairing when a closed-channel Fermi sea is introduced.

Potential relevance to unconventional superfluidity and superconductivity.

Abstract

The orbital Feshbach resonance (OFR) is a novel scheme for magnetically tuning the interactions in closed-shell fermionic atoms. Remarkably, unlike the Feshbach resonances in alkali atoms, the open and closed channels of the OFR are only very weakly detuned in energy. This leads to a unique effect whereby a medium in the closed channel can Pauli block, or frustrate, the two-body scattering processes. Here, we theoretically investigate the impact of frustration in the few- and many-body limits of the experimentally accessible three-dimensional $^{173}$Yb system. We find that by adding a closed-channel atom to the two-body problem, the binding energy of the ground state is significantly suppressed, and by introducing a closed-channel Fermi sea to the many-body problem, we can drive the system towards weaker fermion pairing. These results are potentially relevant to superconductivity in solid-state multiband materials, as well as to the current and continuing exploration of unconventional Fermi-gas superfluids near the OFR.

Equations

$$\begin{gathered} H_0=\sum _{\mathbf{k}}[{ϵ}_{\mathbf{k}}\left(g_{\uparrow \mathbf{k}}^{†}{g}_{\uparrow \mathbf{k}}+e_{\downarrow \mathbf{k}}^{†}{e}_{\downarrow \mathbf{k}}\right) \\ +({ϵ}_{\mathbf{k}}+\delta /2)\left(g_{\downarrow \mathbf{k}}^{†}{g}_{\downarrow \mathbf{k}}+e_{\uparrow \mathbf{k}}^{†}{e}_{\uparrow \mathbf{k}}\right)], \end{gathered}$$

$$V=\sum _{\mathbf{Q},\mathbf{k},{\mathbf{k}}^{\prime }}\left(U_S{S}_{\mathbf{Q}\mathbf{k}}^{†}{S}_{\mathbf{Q}{\mathbf{k}}^{\prime }}+U_T{T}_{\mathbf{Q}\mathbf{k}}^{†}{T}_{\mathbf{Q}{\mathbf{k}}^{\prime }}\right),$$

$$S_{\mathbf{Q}\mathbf{k}}^{†}$$

$$T_{\mathbf{Q}\mathbf{k}}^{†}$$

$$\frac{m}{4\pi a_{S,T}}=\frac{1}{U_{S,T}}+\sum _{\mathbf{k}}^{\Lambda }\frac{1}{2{ϵ}_{\mathbf{k}}},$$

$$\displaystyle V=\sum_{{\mathbf{Q}},\,\mathbf{k},\,\mathbf{k}^{\prime}}\Bigl{[}U_{+}\Bigl{(}$$

$$+$$

$$\displaystyle+\,\,U_{-}\Bigl{(}$$

$$+$$

$$T(E)=V+V\frac{1}{E-H_0+i0}T(E),$$

$$\mathbb{T}(E)=\mathbb{V}+\mathbb{VR}\sum _{\mathbf{k}}\left(\begin{array}{cc}\frac{1}{E-2{ϵ}_{\mathbf{k}}+i0}& 0\\ 0& \frac{1}{E-2{ϵ}_{\mathbf{k}}-\delta +i0}\end{array}\right)\mathbb{RT}(E),$$

$$\displaystyle\mathbb{V}=\left(\begin{array}[]{cc}U_{S}&0\\ 0&U_{T}\end{array}\right)\quad\mathrm{and}\quad\mathbb{R}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{cc}1&1\\ 1&-1\end{array}\right)$$

$$\mathbb{T}(E)={\left[{\mathbb{V}}^{-1}-\mathbb{R}\sum _{\mathbf{k}}\left(\begin{array}{cc}\frac{1}{E-2{ϵ}_{\mathbf{k}}+i0}& 0\\ 0& \frac{1}{E-2{ϵ}_{\mathbf{k}}-\delta +i0}\end{array}\right)\mathbb{R}\right]}^{-1}.$$

$$\mathbb{T}(E)=$$

$$\frac{4\pi }{m}{\left[{\left(\begin{array}{cc}{a}_S& 0\\ 0& a_T\end{array}\right)}^{-1}-\mathbb{R}\left(\begin{array}{cc}\sqrt{-mE}& 0\\ 0& \sqrt{m(\delta -E)}\end{array}\right)\mathbb{R}\right]}^{-1},$$

$$\sum _{\mathbf{k}}\left(\frac{1}{E-2{ϵ}_{\mathbf{k}}+i0}+\frac{1}{2{ϵ}_{\mathbf{k}}}\right)=\frac{m}{4\pi }\sqrt{-mE}.$$

$$\displaystyle f(E)=-\frac{m}{8\pi}\left(\begin{array}[]{cc}1&1\end{array}\right)\mathbb{T}(E)\left(\begin{array}[]{c}1\\ 1\end{array}\right)=\frac{a_{S}+a_{T}-2a_{S}a_{T}\sqrt{m(\delta-E)}}{(a_{S}+a_{T})\bigl{[}\sqrt{-mE}+\sqrt{m(\delta-E)}\,\bigr{]}-2a_{S}a_{T}\sqrt{-mE}\sqrt{m(\delta-E)}-2}\hskip 1.49994pt,$$

$$a=-f(0)=\frac{2a_S{a}_T\sqrt{m\delta }-a_S-a_T}{(a_S+a_T)\sqrt{m\delta }-2}.$$

$$\displaystyle\mathrm{det}\bigl{[}\mathbb{T}(-\varepsilon_{b})\bigr{]}^{-1}=0\hskip 1.99997pt.$$

$$|\Psi\rangle=\sum_{{\mathbf{n}_{1}},\,{\mathbf{n}_{2}},\,{\mathbf{n}_{3}}}\bigl{(}\lambda_{{\mathbf{n}_{1}}{\mathbf{n}_{2}}{\mathbf{n}_{3}}}|g_{{\downarrow}{\mathbf{n}_{1}}},e_{{\downarrow}{\mathbf{n}_{2}}},g_{{\uparrow}{\mathbf{n}_{3}}}\rangle\bigr{.}\\ \bigl{.}+\,\,\tfrac{1}{2}\gamma_{{\mathbf{n}_{1}}{\mathbf{n}_{2}}{\mathbf{n}_{3}}}|g_{{\downarrow}{\mathbf{n}_{1}}},e_{{\uparrow}{\mathbf{n}_{2}}},g_{{\downarrow}{\mathbf{n}_{3}}}\rangle\bigr{)}\hskip 1.00006pt,$$

$$\begin{gathered} H_0=\sum _{\mathbf{n}}[{\epsilon }_{\mathbf{n}}\left(g_{\uparrow \mathbf{n}}^{†}{g}_{\uparrow \mathbf{n}}+e_{\downarrow \mathbf{n}}^{†}{e}_{\downarrow \mathbf{n}}\right) \\ +({\epsilon }_{\mathbf{n}}+\delta /2)\left(g_{\downarrow \mathbf{n}}^{†}{g}_{\downarrow \mathbf{n}}+e_{\uparrow \mathbf{n}}^{†}{e}_{\uparrow \mathbf{n}}\right)]. \end{gathered}$$

$$V=\sum _{\mathbf{N},{\mathbf{n}}_r,{\mathbf{n}}_r^{\prime }}\left(U_S{S}_{\mathbf{N}{\mathbf{n}}_r}^{†}{S}_{\mathbf{N}{\mathbf{n}}_r^{\prime }}+U_T{T}_{\mathbf{N}{\mathbf{n}}_r}^{†}{T}_{\mathbf{N}{\mathbf{n}}_r^{\prime }}\right),$$

$$S_{\mathbf{N}{\mathbf{n}}_r}^{†}$$

$$T_{\mathbf{N}{\mathbf{n}}_r}^{†}$$

$$\det \left[\mathbb{R}{\left(\begin{array}{cc}{a}_S& 0\\ 0& a_T\end{array}\right)}^{-1}\mathbb{R}-\left(\begin{array}{cc}{\mathbb{A}}_E& 0\\ 0& {\mathbb{A}}_{E-\delta }-{\mathbb{B}}_{E-\delta }\end{array}\right)\right]=0,$$

$${\mathbb{A}}_E^{\rho ,{\rho }^{\prime }}$$

$$=\sqrt{2}\frac{\Gamma [(\rho -E/2)/\omega ]}{\Gamma [(\rho -E/2-1/2)/\omega ]}{\delta }_{\rho ,{\rho }^{\prime }},$$

$${\mathbb{B}}_E^{\rho ,{\rho }^{\prime }}$$

$${\mathcal{H}}_0=$$

$$\displaystyle+\Bigl{(}\xi_{\mathbf{k}}^{e\uparrow}+\delta/2\Bigr{)}e_{\uparrow\mathbf{k}}^{\dagger}e_{\uparrow\mathbf{k}}\Bigr{]}\hskip 0.50003pt,$$

$${\mu }_{g\uparrow }$$

$${\mu }_{e\downarrow }$$

$${\tilde{\Delta }}_o$$

$${\tilde{\Delta }}_c$$

$${\mathcal{H}}^{\prime }=\frac{U_{+}\left({\Delta }_o^2+{\Delta }_c^2\right)-2U_{-}{\Delta }_o{\Delta }_c}{U_{-}^2-U_{+}^2}+\sum _{\mathbf{k}}\left(2{\xi }_{\mathbf{k}}+\delta /2+h\right)$$

$$\displaystyle+\sum_{\mathbf{k}}\psi_{\mathbf{k}}^{\dagger}\left(\begin{array}[]{cccc}\xi_{\mathbf{k}}&\Delta_{o}&0&0\\ \Delta_{o}&-\xi_{\mathbf{k}}&0&0\\ 0&0&\xi_{\mathbf{k}}+\delta/2-h&\Delta_{c}\\ 0&0&\Delta_{c}&-\xi_{\mathbf{k}}-\delta/2-h\\ \end{array}\right)\psi_{\mathbf{k}}\hskip 1.49994pt,$$

Full text

I Introduction

Ultracold atomic gases are a testbed for new quantum physics owing to their remarkable tunability. Factors such as dimensionality, density, and the atoms’ internal states and interactions, can all be independently controlled and measured. This flexibility has been enabled by concurrent experimental advances including the manipulation of optical potentials, radiofrequency spectroscopic techniques and Feshbach resonances Lewenstein et al. (2007); Bloch et al. (2008); Chin et al. (2010). In Bose vapours, the achievement of Bose-Einstein condensation Anderson et al. (1995); Bradley et al. (1995); Davis et al. (1995) has led notably to observations of the Mott insulator-superfluid phase transition Jaksch et al. (1998); Greiner et al. (2002) and bosonic Luttinger liquids Paredes et al. (2004); Kinoshita et al. (2004). On the other hand, the realization of degeneracy in Fermi gases DeMarco and Jin (1999); Truscott et al. (2001); Schreck et al. (2001) has allowed for demonstrations of the crossover from a Bardeen-Cooper-Schrieffer (BCS) superfluid of weakly bound Cooper pairs to a molecular Bose-Einstein condensate (BEC) Regal et al. (2004); Bartenstein et al. (2004); Zwierlein et al. (2004); Bourdel et al. (2004). Studies of the latter have incorporated systems with unequal populations for the spin-up and spin-down components Zwierlein et al. (2006); Partridge et al. (2006), which feature much richer phase diagrams compared to the equal-spin case Parish et al. (2007); Sheehy and Radzihovsky (2007).

Feshbach resonances are a tool for controlling the interparticle interactions in a quantum fluid Chin et al. (2010). The basic premise relies on changing the relative energy between two coupled and pairwise potentials that vary differently with an applied magnetic field due to their unequal magnetic moments. When the threshold of an energetically accessible potential, or open channel, is tuned to match a bound state in an energetically closed channel, the scattering length becomes strongly dependent on the field. This allows the system to be brought into the unitarity regime of strongest possible interactions. For alkali gases, the two scattering channels correspond to whether the total spin from the unpaired electrons is a singlet or triplet. By contrast, and owing to their completely filled outer shells, the total electronic spin of alkaline-earth(-like) atoms in the ground state is zero. Hence, these species were long thought to be inert against such tuning techniques.

However, in 2015, Zhang et al. Zhang et al. (2015) proposed a different kind of Feshbach resonance for closed-shell fermions Note (1) by considering both electronic-orbital and nuclear-spin angular-momentum degrees of freedom. Their proposal was experimentally feasible since, in addition to the ground state ${}^{1}S_{0}$, the alkaline-earth system features a long-lived first excited state ${}^{3}P_{0}$, where one electron is excited to a $p$ orbital and the two valence electrons form a spin triplet Note (2). In this case, the resonance channels are coupled by so-called interorbital nuclear-spin exchange interactions Scazza et al. (2014); Cappellini et al. (2014); Zhang et al. (2014) which depend on whether the orbital degree of freedom is symmetric or antisymmetric. Thus, the mechanism was named the orbital Feshbach resonance (OFR) Zhang et al. (2015).

Shortly after its prediction, the OFR was experimentally observed in a gas of 173Yb atoms by two separate groups Höfer et al. (2015); Pagano et al. (2015). This has permitted the experimental study of strongly interacting multiorbital Fermi systems, including polarons Darkwah Oppong et al. (2019) and ultracold molecules Cappellini et al. (2019). Concurrently, the OFR has attracted much attention as a promising candidate for realizing a Fermi superfluid in nonalkali atomic gases Zhang et al. (2015); Xu et al. (2016); He et al. (2016); Mondal et al. (2018a, b); Wang et al. (2017); Cheng et al. (2016); Iskin (2017); Cheng et al. (2017); Zhang et al. (2017); Deng et al. (2017); Zou et al. (2018). The orbital degree of freedom has furthermore led to progress towards the implementation of the Kondo model in cold atoms Riegger et al. (2018).

For the familiar case of alkali resonances, the detuning between the open and closed channels far exceeds all other energy scales in the problem Chin et al. (2010). Consequently, the closed-channel state can be treated as a structureless boson and is only virtually involved in the collisional process. However, because the nuclear Zeeman effect is much smaller than its electronic counterpart, the channels of the OFR are only very weakly detuned Zhang et al. (2015); Boyd et al. (2007). As a result, an inert medium in the closed channel can effectively Pauli block — or frustrate — atoms scattering from the open channel. This possibility of frustrating the interactions is truly unique to the OFR but has so far only been considered in the polaron scenario Höfer (2017); Darkwah Oppong et al. (2019); Chen et al. (2018).

In this work, we theoretically study the impact of frustration on both the few- and the many-body physics near the OFR in three-dimensional 173Yb. We first introduce an extra closed-channel atom to the two-body problem and we solve for the new energy spectrum. Importantly, we observe a sizeable suppression of the ground-state binding energy for the frustrated system, which would not occur if the closed channel was strongly detuned. Second, we consider the effect of a closed-channel Fermi sea on the mean-field BCS-BEC crossover physics at zero temperature. Complementing our three-body result, we see that as the Fermi sea is enlarged, both the average chemical potential and the pairing order parameters for the resonance channels move increasingly towards the BCS side of weaker pairing. We discuss the significance and applications of this new closed-channel physics to current experiments.

This paper is organized as follows: In Sec. II, we introduce our model for studying the orbital Feshbach resonance in 173Yb. We conduct our few- and many-body studies on frustration in Secs. III and IV, respectively, and we conclude in Sec. V.

II Orbital Feshbach Resonance

We describe the orbital Feshbach resonance between a pair of alkaline-earth(-like) fermionic atoms in free space by following Zhang et al. Zhang et al. (2015) (also see Refs. Xu et al. (2016) and Cheng et al. (2016)). First, we take as our degrees of freedom two electronic-orbital states, “ground” ${}^{1}S_{0}\equiv|g\rangle$ and “excited” ${}^{3}P_{0}\equiv|e\rangle$, and two nuclear-spin states, $m_{F,1}\equiv|{\uparrow}\rangle$ and $m_{F,2}\equiv|{\downarrow}\rangle$ Note (6). Then the open channel is designated by one atom in the state $|g{\uparrow}\rangle$ with the other in $|e{\downarrow}\rangle$, while the closed channel has one atom in $|g{\downarrow}\rangle$ and the other in $|e{\uparrow}\rangle$ Note (5). We write them as $|g{\uparrow},e{\downarrow}\rangle$ and $|g{\downarrow},e{\uparrow}\rangle$, respectively (and we do not need to antisymmetrize these states since the fermions are distinguishable). Both channels have the same threshold energies in the absence of a magnetic field $B$; however, under a finite field they separate by an amount $\delta\equiv\delta(B)$ as shown in Fig. 1. This effect is attributed to the small difference in nuclear Landé $g$ factor between the $|g\rangle$ and $|e\rangle$ orbital levels, which originates from a weak hyperfine coupling of the ${}^{3}P_{0}$ to the ${}^{3}P_{1}$ state Boyd et al. (2007). When the interchannel separation $\delta$ is magnetically tuned to match the binding energy of a molecule in the closed channel, the open-channel scattering length diverges, giving rise to the OFR.

Now the noninteracting Hamiltonian, which governs the system when the atoms are far apart, must be diagonal in the open- and closed-channel basis and is given by

[TABLE]

where the single-particle dispersion $\epsilon_{\mathbf{k}}=|\mathbf{k}|^{2}/(2m)\equiv k^{2}/(2m)$ and $m$ is the atomic mass Note (8). Here, $g_{\sigma\mathbf{k}}^{\dagger}$ and $e_{\sigma\mathbf{k}}^{\dagger}$ ($g_{\sigma\mathbf{k}}$ and $e_{\sigma\mathbf{k}}$) are the operators that create (annihilate) a fermion in the ground ($g$) or excited ($e$) electronic orbital state, with three-dimensional (3D) momentum vector $\mathbf{k}$ and nuclear pseudospin $\sigma\in\{\uparrow,\,\downarrow\}$. Note that, in this work, we employ units where both the system volume and the reduced Planck constant $\hbar$ are set to unity.

When the atoms collide, fermionic statistics imply that the relative wave function is antisymmetric under their exchange. We know that the spatial component of the wave function is symmetric, because at ultralow temperatures we only have spherical ($s$-wave) scattering. Correspondingly the internal state, which encompasses both orbital and (nuclear) spin degrees of freedom, must be overall antisymmetric. This allows for two types of scattering: one orbital triplet (nuclear-spin singlet) channel and another orbital singlet (nuclear-spin triplet) channel Note (7).

One can therefore diagonalize the interaction potential in orbital singlet ($S$) and triplet ($T$) states as follows:

[TABLE]

with the operators

[TABLE]

and the corresponding interaction strengths, $U_{S}$ and $U_{T}$. This describes the scattering of two particles with center-of-mass momentum ${\mathbf{Q}}$, and relative momenta ${\mathbf{k}^{\prime}}$ and $\mathbf{k}$ before and after the collision, respectively. The full Hamiltonian is then $H=H_{0}+V$.

By the process of renormalization, the two parameters of the model $\{U_{S},\,U_{T}\}$ can be related to the physical parameters of low-energy collisions, the $s$-wave scattering lengths $\{a_{S},\,a_{T}\}$. This yields the relations

[TABLE]

where we set the cutoff $\Lambda$ to infinity at the end of the calculation, while the singlet and triplet scattering lengths Note (4) are fixed by the properties of the atoms.

When rotated into the same basis as $H_{0}$, the interaction potential takes the form

[TABLE]

with $U_{\pm}=(U_{S}\pm U_{T})/2$. The $U_{+}$ term denotes an intrachannel potential for both open and closed interactions, while $U_{-}$ describes an interchannel coupling as illustrated in Fig. 1.

From the Lippmann-Schwinger equation Lippmann and Schwinger (1950), we now construct the $T$-matrix operator:

[TABLE]

where $i0$ denotes an infinitesimal positive imaginary quantity and $E$ the energy. Importantly, in the center-of-mass frame, the matrix representation of Eq. (6) is independent of incoming and outgoing relative momentum. We can write this conveniently in the singlet and triplet basis,

[TABLE]

where

[TABLE]

are, respectively, the interaction matrix and the involutory matrix that transforms between the singlet-triplet and open-closed channel bases.

Equation (7) may be readily solved to give

[TABLE]

The interactions are then renormalized according to the prescription in Eq. (4), yielding

[TABLE]

where we have used the fact that

[TABLE]

To obtain the scattering amplitude $f$ relevant to experiment, we project the $T$ matrix onto the open channel:

[TABLE]

and this immediately affords the open-channel $s$-wave (zero-energy) scattering length,

[TABLE]

Here, we see that $a$ diverges at \delta=4/\bigl{[}m(a_{S}+a_{T})^{2}\bigr{]}\equiv\delta_{\mathrm{res}}, which denotes the magnetic-field position of the OFR.

The pole of the $T$ matrix, Eq. (II), provides the condition for a two-body bound state with binding energy $\varepsilon_{b}$:

[TABLE]

Note, this expression is equivalent to requiring that the open-channel scattering amplitude, Eq. (12), diverges. We thus find two molecular states associated with the two fixed scattering lengths. At zero detuning these have energies \varepsilon_{b}=1/\bigl{(}ma_{S}^{2}\bigr{)} and \varepsilon_{b}=1/\bigl{(}ma_{T}^{2}\bigr{)}.

For most alkaline-earth atoms, the difference in nuclear magnetic moment between the OFR channels, which allows them to be tuned, is very small compared to the alkali case. Fortunately, in 173Yb, $a_{T}$ is unusually large and provides a shallow bound state that can be accessed with magnetic field strengths typical to experiments Cappellini et al. (2014); Scazza et al. (2014). Conversely, 173Yb also features an extremely small value of $a_{S}>0$ that gives rise to a very deeply bound state within our model. It is worth emphasizing that the latter is inaccessible by experiments and not responsible for the OFR. Moreover, because we are considering low-energy physics, our model cannot properly treat the physics associated with this deeper state. Thus, we ignore it throughout our paper.

III Three-Body Problem

As a first study of frustrated orbital Feshbach resonances, we consider the scenario where two fermionic 173Yb atoms start out in either the open channel ($g{\uparrow},\,e{\downarrow}$) or the closed channel ($g{\downarrow},\,e{\uparrow}$) and we add a third atom to the closed channel (say $g{\downarrow}$). We allow the particles to be tightly confined by a three-dimensional and isotropic harmonic-oscillator potential, with frequency $\omega$ Note (3), which introduces an effective interparticle spacing set by the harmonic-oscillator length, $\sqrt{1/m\omega}$. Such a model can mimic, in effect, a many-body system where the physics is governed by the scale of the interparticle separation. Now when two of the atoms, with 3D harmonic-oscillator indices ${\mathbf{n}_{1}}$ and ${\mathbf{n}_{2}}$, scatter from the open channel into the closed channel, the extra atom restricts the number of harmonic-oscillator levels into which they can scatter (see Fig. 2). Here, we investigate this Pauli-blocking effect on two-body collisions for the full range of interaction strengths measured in Ref. Höfer et al. (2015).

We start from the most general three-body state,

[TABLE]

with basis states $|g_{{\downarrow}{\mathbf{n}_{1}}},e_{{\downarrow}{\mathbf{n}_{2}}},g_{{\uparrow}{\mathbf{n}_{3}}}\rangle\equiv g_{{\downarrow}{\mathbf{n}_{1}}}^{\dagger}e_{{\downarrow}{\mathbf{n}_{2}}}^{\dagger}g_{{\uparrow}{\mathbf{n}_{3}}}^{\dagger}|0\rangle$ and $|g_{{\downarrow}{\mathbf{n}_{1}}},e_{{\uparrow}{\mathbf{n}_{2}}},g_{{\downarrow}{\mathbf{n}_{3}}}\rangle\equiv g_{{\downarrow}{\mathbf{n}_{1}}}^{\dagger}e_{{\uparrow}{\mathbf{n}_{2}}}^{\dagger}g_{{\downarrow}{\mathbf{n}_{3}}}^{\dagger}|0\rangle$, where now the subscripts on the creation operators indicate the harmonic-oscillator indices, and $|0\rangle$ is the vacuum. Explicitly, ${\mathbf{n}_{i}}\equiv\{\rho_{i},\,l_{i},\,\varsigma_{i}\}$ comprises the radial $(\rho=0,\,1,\,...,\,\infty)$, angular-momentum $(l=0,\,1,\,...,\,\rho-1)$, and projection $(\varsigma=-l,\,-l+1,\,...,\,l-1,\,l)$ harmonic-oscillator quantum numbers for the atom at position $\mathbf{x}_{i}$. The objects, $\lambda_{{\mathbf{n}_{1}}{\mathbf{n}_{2}}{\mathbf{n}_{3}}}=\langle g_{{\downarrow}{\mathbf{n}_{1}}},e_{{\downarrow}{\mathbf{n}_{2}}},g_{{\uparrow}{\mathbf{n}_{3}}}|\Psi\rangle$ and $\gamma_{{\mathbf{n}_{1}}{\mathbf{n}_{2}}{\mathbf{n}_{3}}}=\langle g_{{\downarrow}{\mathbf{n}_{1}}},e_{{\uparrow}{\mathbf{n}_{2}}},g_{{\downarrow}{\mathbf{n}_{3}}}|\Psi\rangle$, are the wave-function amplitudes and Fermi statistics require $\gamma_{{\mathbf{n}_{1}}{\mathbf{n}_{2}}{\mathbf{n}_{3}}}=-\,\gamma_{{\mathbf{n}_{3}}{\mathbf{n}_{2}}{\mathbf{n}_{1}}}$.

In the presence of the confining potential, the non-interacting Hamiltonian in Eq. (1) is modified to

[TABLE]

Here, $\varepsilon_{{\mathbf{n}}}=(2\rho+l)\omega$ is the single-particle energy eigenvalue for the harmonic oscillator (and we ignore the energy of $\frac{3}{2}\omega$ corresponding to zero-point motion).

The short-range pairwise interactions of the system are diagonal in the basis where the orbital degree of freedom forms a singlet or triplet Zhang et al. (2015). Since they depend only on the relative motion, we describe them by transforming from the coordinates of two colliding atoms, ${\mathbf{x}_{i}}$ and ${\mathbf{x}_{j}}$, to those of their center-of-mass and relative positions, ${\mathbf{X}}=({\mathbf{x}_{i}}+{\mathbf{x}_{j}})/\sqrt{2}$ and ${\mathbf{x}_{r}}=({\mathbf{x}_{j}}-{\mathbf{x}_{i}})/\sqrt{2}$, respectively Note (9). These have the corresponding harmonic-oscillator indices ${\mathbf{N}}$ and ${\mathbf{n}_{r}}$. We then recast the interacting part, Eq. (2), of the Hamiltonian as

[TABLE]

with the singlet ($S$) and triplet ($T$) operators

[TABLE]

At low energy, we only have scattering in the $s$-wave ($l_{r}=\varsigma_{r}=0$) because the centrifugal barrier for higher partial waves prevents particles from approaching each other at short distance. Consequently, $\phi_{\rho_{r}}$ above denotes the relative $s$-wave harmonic-oscillator eigenfunction at zero separation (see Appendix A). The mentioned coordinate transformation is performed via the two-body Clebsch-Gordan coefficient $\langle{\mathbf{n}_{1}},{\mathbf{n}_{2}}|{\mathbf{N}},{\mathbf{n}_{r}}\rangle$, which has the selection rule ${\mathbf{n}_{1}}+{\mathbf{n}_{2}}={\mathbf{N}}+{\mathbf{n}_{r}}$ Smirnov (1962).

Now that we have set up the problem, we relegate the details of its solution to Appendix A and only present the final equations here. By working in the center-of-mass frame and specializing to the sector of zero total angular momentum, we find that the energy $E$ of the system satisfies

[TABLE]

where $\{\mathbb{A}_{E},\,\mathbb{B}_{E}\}$ are square matrices given by

[TABLE]

and $\mathbb{R}$ is the rotation matrix of Eq. (8). The summation indices $\{\kappa,\,\kappa^{\prime}\}$ refer to the radial quantum numbers for relative atom-atom motion between the interacting particles, while the matrix indices $\{\rho,\,\rho^{\prime}\}$ are the radial quantum numbers for relative atom-pair motion between the interacting particles’ center of mass and the noninteracting atom. Equation (19) describes both open-channel scattering (through the two-body term $\mathbb{A}_{E}$) and closed-channel scattering ($\mathbb{A}_{E-\delta}-\mathbb{B}_{E-\delta}$). When all three atoms are in the closed channel, the diagonal matrix $\mathbb{A}_{E-\delta}$ is concerned with one of the two interacting pairs. The nondiagonal kernel $\mathbb{B}_{E-\delta}$ contains the three-body Clebsch-Gordan coefficient $\langle\kappa,\rho|\kappa^{\prime},\rho^{\prime}\rangle$ (with $\kappa+\rho=\kappa^{\prime}+\rho^{\prime}$), which then allows us to rotate our relative coordinate axes from one interacting pair to the other. Pauli blocking from the closed channel is manifested by the minus sign in $\mathbb{A}_{E-\delta}-\mathbb{B}_{E-\delta}$, which follows from the symmetry property of the amplitudes, $\gamma_{{\mathbf{n}_{1}}{\mathbf{n}_{2}}{\mathbf{n}_{3}}}=-\,\gamma_{{\mathbf{n}_{3}}{\mathbf{n}_{2}}{\mathbf{n}_{1}}}$, in Eq. (9). The matrix elements of $\mathbb{A}_{E}$ in Eq. (III) were first treated by Ref. Busch et al. (1998), while a method for the efficient evaluation of $\mathbb{B}_{E}$ is provided in Appendix A.

Equations (19)-(III) can be solved numerically efficiently and the resulting energy spectrum is shown by the solid blue lines ( ) in Fig. 3. To provide a comparison, the dashed green lines ( ) give the energies for two atoms interacting in a 3D harmonic trap, which were first calculated by Ref. Höfer et al. (2015). The deep equidistant bound states in the three-body problem arise from the small and positive singlet scattering length of 173Yb, and are not relevant to the physics of interest, as explained in Sec. II. Instead, we concentrate on the ground state in the region of low-energy physics (orange line, , and highlighted by the box). Significantly, we see that the binding energy of the ground state is strongly suppressed due to closed-channel Pauli blocking — by $\sim 0.3\,\omega$ at the OFR (gray dashed line, ). This physics could, in principle, be probed with clock spectroscopy in a deep optical lattice, similar to the experiment that realized the OFR Höfer et al. (2015). The effects of frustration are stronger and may therefore be easier to observe in this scenario than in free space, where an observation was previously attempted Darkwah Oppong et al. (2019). However, it remains to be seen whether one would encounter difficulties stemming from three-body losses Hofrichter et al. (2016).

IV Many-Body Problem

IV.1 Theory

Now that we have demonstrated how the few-body correlations can be tailored by adding an atom to the closed channel, we consider how frustration manifests at the many-body level. To this end, we modify the standard theory Giorgini et al. (2008) for describing the zero-temperature BCS-BEC crossover to accommodate a single Fermi sea (say $g{\downarrow}$) in the closed channel near an OFR. Intuitively, this means that atoms can no longer scatter into the closed channel at an arbitrary momentum, but can only do so within a restricted range of momenta which is set by the Fermi sea. We stress that the case we consider is different from the more standard scenario of a spin imbalance in the open channel Chen et al. (2016); Deng et al. (2018); Xu and Qi (2018).

The Hamiltonian, $\mathcal{H}=\mathcal{H}_{0}+V$, for this system consists of the noninteracting component,

[TABLE]

and the interaction potential, $V$ in Eq. (II). Above, $\xi_{\mathbf{k}}^{\alpha}=\epsilon_{\mathbf{k}}-\mu_{\alpha}$ is the single-particle kinetic energy measured relative to the chemical potential $\mu_{\alpha}$, while the fermionic operators are explained below Eq. (1). Since we have an interconversion process (see Fig. 1) between the open and closed channels, they must be in chemical equilibrium: $\mu_{g\uparrow}+\mu_{e\downarrow}=\mu_{g\downarrow}+\mu_{e\uparrow}$. However, this does not place any restriction on the differences, $\mu_{g\uparrow}-\mu_{e\downarrow}$ and $\mu_{g\downarrow}-\mu_{e\uparrow}$. To allow for the possibility of a closed-channel Fermi sea, we therefore define the average chemical potential $\mu$ and the ‘Zeeman’ field $h$ such that

[TABLE]

Similar to the approach in Ref. Zhang et al. (2015), we now define the real quantities,

[TABLE]

and we perform a mean-field decoupling of the interaction term, $V$. This leads to the following reduced Hamiltonian:

[TABLE]

with \psi_{\mathbf{k}}^{\dagger}=\bigl{(}g_{\uparrow\mathbf{k}},\,e_{\downarrow-\mathbf{k}}^{\dagger},\,g_{\downarrow\mathbf{k}},\,e_{\uparrow-\mathbf{k}}^{\dagger}\bigr{)} and $\xi_{\mathbf{k}}=\epsilon_{\mathbf{k}}-\mu$. Here, we have defined the pairing order parameters for the open ($o$) and closed ($c$) channels, $\Delta_{o}=U_{+}\widetilde{\Delta}_{o}+U_{-}\widetilde{\Delta}_{c}$ and $\Delta_{c}=U_{-}\widetilde{\Delta}_{o}+U_{+}\widetilde{\Delta}_{c}$. Notice that although $\Delta_{o}$ and $\Delta_{c}$ of Eq. (IV.1) appear to be completely decoupled, they are in fact coupled by the relations in Eq. (IV.1).

We now diagonalize $\mathcal{H}^{\prime}$ by using the Bogoliubov transformation:

[TABLE]

where $\{\beta_{\sigma\mathbf{k}},\,\zeta_{\sigma\mathbf{k}}\,;\,\sigma\in\{\uparrow,\,\downarrow\}\}$ are the operators for (fer-mionic) quasiparticle excitations in both channels, while $\{u_{\mathbf{k}},\,v_{\mathbf{k}}\}$ are real coefficients satisfying $u^{2}_{\mathbf{k}}+v^{2}_{\mathbf{k}}=1$ (and similarly for the primed variables). Eventually, this yields the new Hamiltonian,

[TABLE]

The expression above contains the open- and closed-channel excitation energies,

[TABLE]

In the absence of spin imbalance, the zero-temperature grand potential is

[TABLE]

In this case, the eigenenergies for excitations are simply $\{E_{\mathbf{k}}^{o},\,E_{\mathbf{k}}^{c}\}$ and there are no unpaired quasiparticles in the ground state. However, in the presence of the Zeeman field, we begin populating quasiparticle modes in the ground state once $|h|>\mathrm{min}_{\mathbf{k}}E_{\mathbf{k}}^{c}$. This condition defines a region of momenta, between $k_{\mathrm{min}}=\left\{2m\left[\mu-\delta/2-\left(h^{2}-\Delta_{c}^{2}\right)^{1/2}\right]\right\}^{1/2}$ and $k_{\mathrm{max}}=\left\{2m\left[\mu-\delta/2+\left(h^{2}-\Delta_{c}^{2}\right)^{1/2}\right]\right\}^{1/2}$, where the pairing gap is removed and we build up an unpaired Fermi sea in the closed channel (see Fig. 4). Under these considerations, we write the grand potential for the spin-imbalanced system as

[TABLE]

where the unit step function, $\Theta(\hskip 1.4pt\cdots)$, selects out the range of states with negative quasiparticle energies. (Since this formulation is symmetric with respect to a change in sign of the Zeeman field, we subsequently take $h>0$.)

To study the crossover from BEC to BCS pairing, we minimize the grand potential in Eq. (IV.1) with respect to both order parameters $\{\Delta_{o},\,\Delta_{c}\}$ for given $\{\mu,\,h,\,\delta\}$. By taking $\partial\Omega/\partial\Delta_{o}=0$ and $\partial\Omega/\partial\Delta_{c}=0$ to determine the stationary points, we obtain

[TABLE]

Continuing, we renormalize the interactions in accordance with Eq. (4), and finally, we arrive at the two coupled gap equations:

[TABLE]

with $\chi_{\pm}=a_{S}^{-1}\pm a_{T}^{-1}$ Note (4).

We iteratively solve Eq. (IV.1) and find two solutions: one corresponds to a very deep global minimum of $\Omega$ where $\Delta_{o}$ and $\Delta_{c}$ are in phase, and the other is an excited saddlepoint with the gaps out of phase by $\pi$. The former describes the case when all the particles are in a deep bound state created by the small $a_{S}$ in 173Yb (see Sec. II). Since the chemical potential for this state is large and negative at all values of the detuning He et al. (2016), it is not related to the BCS-BEC crossover near the OFR. Instead, the crossover physics is attributed to the shallow out-of-phase result which is well described at the low-energy scale of our model. We can choose this particular solution by initially setting $U_{S}=0$ in Eq. (IV.1) so that we have $\Delta_{o}=-\Delta_{c}$; then we force the system to adiabatically maintain this phase difference for small $a_{S}$.

Along with the gap equations, we derive two separate number equations from Eq. (IV.1). The average chemical potential $\mu$ is associated with the total number density, $n=n_{o}+n_{c}=-(\partial\Omega/\partial\mu)$, from which we obtain the densities in the open and closed channels,

[TABLE]

Similarly, the extra field $h$ corresponds to the density difference in the closed channel,

[TABLE]

IV.2 Results

We now investigate how the BCS-BEC crossover physics close to the OFR is modified as a result of the closed-channel component in our model. To properly isolate the effect of Pauli blocking, we fix the number of atoms that are involved in pairing,

[TABLE]

and we consider different fixed values of $n_{\mathrm{unpaired}}/n_{\mathrm{paired}}$, where $n_{\mathrm{unpaired}}=\delta n_{c}$ is the number of unpaired atoms in the closed channel. We also define the Fermi energy, $E_{F}=k_{F}^{2}/(2m)$, and the Fermi momentum, k_{F}=\bigl{(}3\pi^{2}n_{\mathrm{paired}}\bigr{)}^{1/3}, with respect to $n_{\mathrm{paired}}$ because this allows a comparison with the spin-balanced system.

In Fig. 5, we plot the crossover for $n_{\mathrm{paired}}=1\times 10^{13}$ atoms per $\mathrm{cm}^{3}$ and with $n_{\mathrm{unpaired}}/n_{\mathrm{paired}}=0$ (red, ), 1 (blue, ), 2 (green, ). These densities are comparable to those achieved in current experiments Höfer et al. (2015); Darkwah Oppong et al. (2019). The upper panel (a) depicts the average chemical potential, $\mu/E_{F}$, against the interchannel detuning, $\delta/\delta_{\mathrm{res}}$, while the lower panel (b) shows the corresponding behaviors of the open- and closed-channel gaps, $|\Delta_{o}|/E_{F}$ and $|\Delta_{c}|/E_{F}$. The red lines ($n_{g{\downarrow}}=n_{e{\uparrow}}$) resemble those in Fig. 3 of Ref. Zhang et al. (2015), although we have used updated values for the scattering lengths Note (4) and have also chosen a different paired-atom density. Notably, here, in the BEC limit of zero detuning, we find that the chemical potential rises by $\sim 0.2\,E_{F}$ when the unpaired-atom density increases by $n_{\mathrm{paired}}$. Furthermore, at most magnetic fields, we see that $\mu/E_{F}$ is visibly shifted upwards, and $\{|\Delta_{o}|,\,|\Delta_{c}|\}/E_{F}$ downwards, for higher $n_{\mathrm{unpaired}}$. This demonstrates how, by introducing a closed-channel medium that frustrates the interactions, we can effectively tune the crossover to suppress pairing.

For additional clarity, in Fig. 6, we show how the atom populations evolve as we sweep the magnetic field through the range in Fig. 5, with $n_{\mathrm{unpaired}}/n_{\mathrm{paired}}=0$ (upper panel), 1 (middle panel), 2 (lower panel). Since we have a fixed total density $n$ in each panel, we observe that increasing the detuning drives atoms from the closed channel into the open channel. Consequently, for large $\delta/\delta_{\mathrm{res}}$, the paired density moves fully into the open channel, $n_{o}\to n_{\mathrm{paired}}$, while the closed channel becomes either empty or completely polarized, $n_{c}\to n_{\mathrm{unpaired}}$. Both channels in the upper plot have equal spin populations: $n_{g{\uparrow}}=n_{e{\downarrow}}$ and $n_{g{\downarrow}}=n_{e{\uparrow}}$. For the two remaining plots, we tune the Zeeman parameter $h/E_{F}$ while increasing $\delta/\delta_{\mathrm{res}}$, so that $n_{\mathrm{unpaired}}$ in the closed channel is fixed to a particular multiple of $n_{\mathrm{paired}}$.

IV.3 BEC Limit

For strong pairing, we can understand the shifts of the chemical potential in Fig. 5 if we treat them as due to an effective interaction between open-channel pairs (or “bosons”) and excess closed-channel atoms (fermions). In the BEC regime where $\delta/\delta_{\mathrm{res}}\to 0$ and $\mu/E_{F}\to-\infty$, we assume that the magnitude of the chemical potential in panel (a) is much greater than the gaps (b). Therefore, we linearize Eq. (IV.1) in the variables $\{\Delta_{o},\,\Delta_{c}\}$, which affords the following implicit function for $\mu$:

[TABLE]

where $\chi_{\pm}$ are defined below Eq. (IV.1). For the spin-balanced system ($\delta n_{c}=0$), the expression above has two solutions: \mu=-1/\bigl{(}2ma_{S}^{2}\bigr{)} and \mu=-1/\bigl{(}2ma_{T}^{2}\bigr{)}. Both of these correspond to binding energies $\varepsilon_{b}$ where $\mu=-\varepsilon_{b}/2$, as required in the deep BEC limit, but only the latter is associated with the crossover (see Sec. II).

Since the density of bosons is fixed across all three cases in Fig. 5, we expand the chemical potential around $-\varepsilon_{b}/2$ in powers of the excess density of fermions, $\delta n_{c}$. To this end, we rewrite the right-hand side of Eq. (35) as a Taylor series in $h-|\mu|$, and then we use Eq. (33) to replace that small parameter by the excess density. To first order this yields

[TABLE]

and solving for the chemical potential we then obtain $\mu\simeq-\varepsilon_{b}/2+g_{\mathrm{BF}}\delta n_{c}$. Here, the binding energy is \varepsilon_{b}=1/\bigl{(}ma_{T}^{2}\bigr{)} and $g_{\mathrm{BF}}=2\pi a_{T}/m$ is the effective coupling constant for the boson-fermion interaction in the Born approximation. These first-order shifts are marked by the orange stars ( ) in panel (a) and well capture the chemical potential behavior for BEC pairing.

V Conclusions and Outlook

We have theoretically studied the orbital Feshbach resonance in the experimentally accessible 173Yb system. Here, the closed channel is available as a set of well resolved intermediate states in the two-atom scattering process. This is in contrast to more familiar alkali resonances, where the closed channel is detuned by an energy much larger than all other scales in the problem. Following from this distinction, we investigated the effect of Pauli-blocking (frustration) caused by a closed-channel medium at both the few- and many-body levels. By harmonically confining three particles to a single lattice site (with one atom always in the closed channel), we found that the lower lying eigenstates are strongly shifted relative to the two-body energy spectrum. In the many-particle continuum picture, we saw that the introduction of a Fermi sea to the closed channel drives the system towards the BCS limit.

Our three-atom prediction could potentially be realized via a clock-spectroscopy scheme very similar to the recent experiment on two-dimensional polarons Darkwah Oppong et al. (2019). Although this experiment attempted the first observation of frustration, the largest Fermi seas that could be achieved, in practice, only led to shifts that were within the experimental uncertainties. By contrast, the magnitude of our energy shifts suggests that an observation might be more feasible in the few-body scenario (on a lattice) than in the free-space system considered by that work. For this reason, we furthermore expect that our many-body result would be enhanced in the context of a state-dependent optical lattice Riegger et al. (2018). The findings of the current study may shed light on BCS pairing in more complicated solid-state multiband materials. Additionally, our model of a Fermi sea in the presence of OFR-induced tunable pairing interactions may assist in the search for exotic superfluids, such as the sought after breached-pair phase Liu and Wilczek (2003).

Acknowledgements.

We are grateful to N. Darkwah Oppong for helpful discussions on the experimental implementation of our three-body results. We acknowledge support from the Australian Research Council (ARC) Centre of Excellence in Future Low-Energy Electronics Technologies (CE170100039). Z.Y.S., M.M.P. and J.L. acknowledge financial support from the ARC via Discovery Project No. DP160102739. J.L. is additionally supported through the ARC Future Fellowship FT160100244.

Appendix A    Derivation of the Three-Body Equations

Here, we work through the method of solving the three-body problem discussed in Sec. III. To parameterize the motion of three equal-mass particles we define two relative coordinates, ${\mathbf{x}_{ij}}=\sqrt{1/2}\,({\mathbf{x}}_{j}-{\mathbf{x}}_{i})$ and ${\mathbf{x}_{k}^{ij}}=\sqrt{2/3}\,[({\mathbf{x}}_{i}+{\mathbf{x}}_{j})/2-{\mathbf{x}}_{k}]$, respectively, between atoms $i$ and $j$, and between their center of mass and atom $k$. The square-root factors ensure the effective masses for the relative motions are both equal to the atomic mass $m$, and throughout this appendix we set $m=\omega=1$. The coordinates, ${\mathbf{x}_{ij}}$ and ${\mathbf{x}_{k}^{ij}}$, are assigned sets of 3D harmonic-oscillator indices, ${\mathbf{n}_{ij}}\equiv\{\rho_{ij},\,l_{ij},\,\varsigma_{ij}\}$ and ${\mathbf{n}_{k}^{ij}}\equiv\{\rho_{k}^{ij},\,l_{k}^{ij},\,\varsigma_{k}^{ij}\}$, as done in the main text. Now, the relative-motion equation decouples from the motion of the three-body center of mass with coordinate ${\mathbf{X}}=({\mathbf{x}}_{i}+{\mathbf{x}}_{j}+{\mathbf{x}}_{k})/3$. Without loss of generality, therefore, we set the corresponding center-of-mass quantum numbers to zero, ${\mathbf{N}}=\{0,\,0,\,0\}$. Since interactions are only in the atom-atom $s$-wave channel, we furthermore have $l_{ij}=\varsigma_{ij}=0$. The atom-pair angular-momentum channels thus decouple, and we consider only the atom-pair $s$-wave states ($l_{k}^{ij}=\varsigma_{k}^{ij}=0$) as these are accessible in experiments. This leads to a vanishing total angular momentum ($L+l_{k}^{ij}+l_{ij}=0$) and, as explained below, we consequently realize a pair of coupled equations depending only on a single index, $\rho_{k}^{ij}$. [Note that $\{\rho_{ij},\,\rho_{k}^{ij}\}$ appear respectively as $\{\kappa,\,\rho\}$ in Eq. (III).]

To begin, we make use of the three-body wave function, $|\Psi\rangle$ in Eq. (9), in addition to the noninteracting and interacting parts of the Hamiltonian, $H_{0}$ and $V$ in Eqs. (10) and (17)-(III). It is convenient to write $V$ in terms of the Huang-Yang pseudopotential Cheng et al. (2016), $4\pi a_{S,\,T}\delta_{\mathrm{reg}}^{(3)}(\mathbf{x}_{i}-\mathbf{x}_{j})$, where the singlet ($S$) and triplet ($T$) scattering lengths are $\{a_{S},\,a_{T}\}$ Note (4) and $\delta_{\mathrm{reg}}^{(3)}(\mathbf{x})\equiv\delta^{(3)}({\mathbf{x}})\partial/\partial x(x\,\cdot\,)$. Then we project the Schrödinger equation, $(E-H_{0})|\Psi\rangle=V|\Psi\rangle$, separately onto both basis states of $|\Psi\rangle$. Following projection onto $|g_{{\downarrow}{\mathbf{n}_{1}}},e_{{\downarrow}{\mathbf{n}_{2}}},g_{{\uparrow}{\mathbf{n}_{3}}}\rangle$, we obtain

[TABLE]

while after projection onto $|g_{{\downarrow}{\mathbf{n}_{1}}},e_{{\uparrow}{\mathbf{n}_{2}}},g_{{\downarrow}{\mathbf{n}_{3}}}\rangle$, we have

[TABLE]

Above, the state $|g_{{\downarrow}{\mathbf{n}_{1}}},e_{{\downarrow}{\mathbf{n}_{2}}},g_{{\uparrow}{\mathbf{n}_{3}}}\rangle$ is defined as the zero of detuning and the single-particle energies of Eq. (10) are replaced by the notation, $E_{{\mathbf{n}}}\equiv\sum_{i\,=\,1}^{3}\varepsilon_{{\mathbf{n}_{i}}}$.

To reduce the number of summation indices we follow the general methods introduced in our earlier works Laird et al. (2017, 2018). The idea is first to transform into the relative-motion frame and then to take the $\delta$-function boundary condition on the wave function that the colliding particles, $i$ and $j$, are superimposed on one another. To this end, we insert a complete set of transformed states, $|{\mathbf{n}_{1}},{\mathbf{n}_{2}},{\mathbf{n}_{3}}\rangle\to|{\mathbf{n}}_{ij},{\mathbf{n}}_{k}^{ij},{\mathbf{N}}\rangle$, on either side of the pseudopotential in Eqs. (A1)-(A). Since the interaction only changes the quantum number for relative atom-atom motion ${\mathbf{n}}_{ij}$, this yields

[TABLE]

where $\widetilde{{\mathbf{n}}}\equiv\{{\mathbf{n}_{1}},\,{\mathbf{n}_{2}},\,{\mathbf{n}_{3}}\}$ (and similarly for the primed variables). Here, $\phi_{\rho_{ij}}=(1/\sqrt{4\pi}\,)\sqrt{2\rho_{ij}!/(\rho_{ij}+1/2)!}\,L_{\rho_{ij}}^{(1/2)}(0)$ are the $s$-wave eigenfunctions for the harmonic oscillator at zero relative separation, and $L$ is the associated Laguerre polynomial.

Next, we reformulate these expressions in terms of reduced atom-pair wave functions that have been subjected to the boundary condition. We define two independent functions,

[TABLE]

by taking atoms 2 and 3 on top of each other. By permuting the particle labels, we also have

[TABLE]

which correspond to when the boundary condition is applied to atoms 1 and 2. After a direct substitution, Eqs. (29)-(29) become

[TABLE]

Continuing, we divide Eq. (A7) by $E-E_{\mathbf{n}}$ and Eq. (30) by $E-E_{\mathbf{n}}-\delta$, and then we act on the left (respectively) with the operators

[TABLE]

We thus arrive at

[TABLE]

After a series of straightforward manipulations Laird et al. (2017, 2018), and by recalling that the angular-momentum ($l,\,m$) quantum numbers are zero, we obtain a matrix equation in terms of the reduced wave functions:

[TABLE]

Above, $\{\eta_{\lambda},\,\eta_{\gamma}\}\equiv\{\eta_{\lambda}^{\rho},\,\eta_{\gamma}^{\rho}\}$ are vectors and $\{\mathbb{A}_{E},\,\mathbb{B}_{E}\}\equiv\{\mathbb{A}_{E}^{\rho,\,\rho^{\prime}},\,\mathbb{B}_{E}^{\rho,\,\rho^{\prime}}\}$ are square matrices indexed by the radial atom-pair quantum number, $\rho\equiv\rho_{k}^{ij}$. Equation (A12) may be solved for the energy $E$ as an eigenvalue problem which we rearrange into Eq. (19) of the main text, with $\{\mathbb{A}_{E},\,\mathbb{B}_{E}\}$ defined by Eq. (III). It only remains to be shown how the nondiagonal kernel $\mathbb{B}_{E}$ can be treated numerically, and the ensuing approach is inspired by Ref. Liu et al. (2010).

By using the notation of Eq. (III), and now replacing $\{\mathbf{x}_{ij},\,\mathbf{x}_{k}^{ij}\}$ with $\{\mathbf{x},\,\mathbf{y}\}$ for clarity, we write the kernel explicitly:

[TABLE]

where the zero-angular-momentum eigenfunctions for the harmonic oscillator are given by

[TABLE]

with $R$ the radial component Smirnov (1962).

For each coordinate in Eq. (A13), there is a Green’s function $G$ which solves the noninteracting Schrödinger equation with a Dirac-delta-function source. As an example, we have (in first quantization)

[TABLE]

for the $\mathbf{x}$ coordinate, where $\nabla_{\mathbf{x}}^{2}$ is the 3D Laplacian operator and the energy $\nu$ is defined relative to the zero-point energy of $\frac{3}{2}$. By expanding both functions in the complete set of harmonic-oscillator eigenfunctions, we obtain

[TABLE]

These relations allow us to convert Eq. (A13) into a one-dimensional integral, as shown:

[TABLE]

and this agrees with the corresponding result in Ref. Liu et al. (2010).

Now, to evaluate the Green’s function, we simply replace $R_{\kappa}(0)$ in Eq. (A16) with its definition from Eq. (A14). Then the identities,

[TABLE]

lead us to write down

[TABLE]

Therefore, we finally arrive at the following equivalent form of the nondiagonal kernel in Eq. (A13):

[TABLE]

where

[TABLE]

Reference Liu et al. (2010) found that for a cutoff $\rho_{\mathrm{max}}>\{\rho,\,\rho^{\prime}\}$, the integral $N_{\kappa,\,\rho,\,\rho^{\prime}}$ vanishes for large enough $\kappa>\kappa_{\mathrm{max}}\sim 4\rho_{\mathrm{max}}$. Numerically, this means that we can tabulate the values of $N$ first and then populate the matrix $\mathbb{B}_{E}$ by a simple and fast summation over $\kappa$. We can subsequently see that the kernel is symmetric, $\mathbb{B}_{E}^{\rho,\,\rho^{\prime}}=\mathbb{B}_{E}^{\rho^{\prime},\,\rho}$.