Universality in a one-dimensional three-body system

TL;DR

This paper investigates the universal behavior of three-body bound states in a one-dimensional heavy-heavy-light system, demonstrating that near zero two-body energy, the three-body states are independent of potential shape.

Contribution

It provides a comparative analysis of zero-range and finite-range potentials, showing universality in three-body bound states near zero two-body energy in one dimension.

Findings

Three-body bound states become universal near zero two-body energy.

Finite-range and zero-range potentials yield similar three-body behaviors.

Numerical and analytical methods confirm potential shape independence.

Abstract

We study a heavy-heavy-light three-body system confined to one space dimension. Both binding energies and corresponding wave functions are obtained for (i) the zero-range, and (ii) two finite-range attractive heavy-light interaction potentials. In case of the zero-range potential, we apply the method of Skorniakov and Ter-Martirosian to explore the accuracy of the Born-Oppenheimer approach. For the finite-range potentials, we solve the Schr\"odinger equation numerically using a pseudospectral method. We demonstrate that when the two-body ground state energy approaches zero, the three-body bound states display a universal behavior, independent of the shape of the interaction potential.

Tables

Table 1. Table 1: Universal constants ϵ n ⋆ superscript subscript italic-ϵ 𝑛 ⋆ \epsilon_{n}^{\star} for different mass ratios commonly used in experiments performed with ultracold atoms.

	atomic mixture $\left(M/m\right)$
	${}_{}{}^{87}\mathrm{Rb}$–${}_{}{}^{40}K$ (2.2)	${}_{}{}^{87}\mathrm{Rb}$–${}_{}{}^7\mathrm{Li}$ (12.4)	${}_{}{}^{133}\mathrm{Cs}$–${}_{}{}^6\mathrm{Li}$ (22.2)
${ϵ}_0^{\star }$	-2.1966	-2.5963	-2.7515
${ϵ}_1^{\star }$	-1.0520	-1.4818	-1.6904
${ϵ}_2^{\star }$	-	-1.1970	-1.3604
${ϵ}_3^{\star }$	-	-1.0377	-1.1479
${ϵ}_4^{\star }$	-	-1.0002	-1.0525
${ϵ}_5^{\star }$	-	-	-1.0040

Equations

$$\left[-\frac{1}{2}\frac{{\text{d}}^2}{\text{d}{x}^2}+v\left(x\right)\right]{\psi }^{(2)}={\mathcal{E}}^{(2)}{\psi }^{(2)}$$

$$v\left(x\right)=v_{0}f(x)$$

$$\left[-\frac{{\alpha }_x}{2}\frac{{\partial }^2}{\partial x^2}-\frac{{\alpha }_y}{2}\frac{{\partial }^2}{\partial y^2}+v\left(r_{+}\right)+v\left(r_{-}\right)\right]\psi =\mathcal{E}\psi$$

$${\alpha }_x\equiv \frac{1+2M/m}{2(1+M/m)}$$

$${\alpha }_y\equiv \frac{2}{1+M/m}$$

$${ϵ}_n\equiv \frac{{\mathcal{E}}_n}{\left|{\mathcal{E}}_{\mathrm{g}}^{(2)}\right|}$$

$$X\equiv x/{\ell }_{\mathrm{g}},Y\equiv y/{\ell }_{\mathrm{g}}$$

$${\ell }_{\mathrm{g}}\equiv \frac{1}{\sqrt{2\left|{\mathcal{E}}_{\mathrm{g}}^{(2)}\right|}},$$

$$\hat{H}\tilde{\psi }=ϵ\tilde{\psi }$$

$$\hat{H}\equiv-\alpha_{x}\frac{\partial^{2}{}}{\partial X^{2}}-\alpha_{y}\frac{\partial^{2}{}}{\partial Y^{2}}-2\ell_{\mathrm{g}}^{2}\left|{v_{0}}\right|\bigg{[}f(\ell_{\mathrm{g}}R_{+})+f(\ell_{\mathrm{g}}R_{-})\bigg{]}$$

$$R_{\pm }\equiv X\pm Y/2.$$

$$f_{\delta }(x)\equiv \delta (x)$$

$${\mathcal{E}}_{\mathrm{g}}^{(2)}=-\frac{1}{2}{v}_0^2$$

$${\ell }_{\mathrm{g}}^2\left|v_0\right|f_{\delta }\left({\ell }_{\mathrm{g}}{R}_{\pm }\right)=\delta (X\pm Y/2),$$

$${\tilde{\psi }}^{(\mathrm{BO})}(X,Y)\equiv \phi (X|Y)\varphi (Y).$$

$$\left\{-{\alpha }_x\frac{{\partial }^2}{\partial X^2}-2\left[\delta (R_{+})+\delta (R_{-})\right]\right\}\phi =u(Y)\phi$$

$$u_{\pm }(Y)=-\frac{1}{{\alpha }_x}{\left[\frac{{\alpha }_x}{\left|Y\right|}{W}_0\left(\pm \frac{\left|Y\right|}{{\alpha }_x}{\mathrm{e}}^{-|Y|/{\alpha }_x}\right)+1\right]}^2$$

$${\phi }_{\pm }(X|Y)=N_{\pm }\left[{\mathrm{e}}^{-\sqrt{|u_{\pm }|/{\alpha }_x}|R_{-}|}\pm {\mathrm{e}}^{-\sqrt{|u_{\pm }|/{\alpha }_x}|R_{+}|}\right],$$

$$\left[-{\alpha }_y\frac{{\partial }^2}{\partial Y^2}+u_{+}(Y)\right]{\varphi }_{+}={ϵ}^{(\mathrm{BO})}{\varphi }_{+},$$

$${\alpha }_x{u}_{+}(Y\to \pm \infty )\to -1,$$

$$n_{\max }\cong \frac{1}{\pi \sqrt{{\alpha }_y}}\int \mathrm{d}Y\sqrt{\left|u_{+}(Y)+\frac{1}{{\alpha }_x}\right|}-\frac{1}{2}$$

$$n_{\max }\cong 0.8781\cdot \sqrt{1+2\frac{M}{m}}-\frac{1}{2}$$

$$G_{ϵ}^{(2)}(X,Y)\equiv -\frac{1}{2\pi \sqrt{{\alpha }_x{\alpha }_y}}{K}_0\left(\sqrt{\left|ϵ\right|}\sqrt{\frac{1}{{\alpha }_x}{X}^2+\frac{1}{{\alpha }_y}{Y}^2}\right)$$

$$\tilde{\psi }(X,Y)=-$$

$$\times \tilde{\psi }(X^{\prime },Y^{\prime })\left[f\left({\ell }_{\mathrm{g}}{R}_{+}^{\prime }\right)+f\left({\ell }_{\mathrm{g}}{R}_{-}^{\prime }\right)\right],$$

$$\tilde{\psi }(X$$

$$\times \tilde{\psi }(X^{\prime },Y^{\prime })\left[\delta \left(X^{\prime }+Y^{\prime }/2\right)+\delta \left(X^{\prime }-Y^{\prime }/2\right)\right].$$

$$\tilde{\psi }(X,Y)=$$

$$+$$

$$\tilde{\psi }(X,-Y)=\pm \tilde{\psi }(X,Y)$$

$$\tilde{\psi }(X,2X)=\int \mathrm{d}{X}^{\prime }{\mathcal{K}}_{ϵ}^{(\pm )}(X,X^{\prime })\tilde{\psi }(X^{\prime },2X^{\prime })$$

$${\mathcal{K}}_{ϵ}^{(\pm )}(X,X^{\prime })\equiv -4$$

Full text

I Introduction

The few-body problem has been of central interest in the physics community since the very beginning of quantum mechanics Born and Oppenheimer (1927); Heitler and London (1927); Bethe and Salpeter (1957); Richter et al. (1993). Continuous efforts have led to theoretical breakthroughs like the Efimov effect Efimov (1970); *Efimov1971; *Efimov1973, that is the appearance of an infinite sequence of universal bound states in the three-dimensional system of three bodies, provided the two-body interactions have a single $s$-wave resonance Nishida (2012). The effect is universal Nielsen et al. (2001); Jensen et al. (2004); Braaten and Hammer (2006); Greene et al. (2017); Naidon and Endo (2017) in the sense that it is independent of the shape of the two-body interaction potential, as long as the latter is tuned to be on $s$-wave resonance.

In the present article we study another class of universal bound states in a three-body system of two identical, heavy particles and a third, light particle, all confined to one spatial dimension (1D) when the heavy-light ground state energy approaches zero. This nearly resonant state is not a virtual state but always weakly bound in the case of an attractive heavy-light interaction. We assume no interaction between the two heavy particles and obtain the binding energies as well as the corresponding wave functions for the zero- and two different finite-range heavy-light interaction potentials. In addition, we prove the universality of these states.

I.1 Dimension of space and symmetry of resonance

The appearance of the Efimov effect crucially depends on the number of spatial dimensions, and on the symmetry of the underlying two-body resonance. Indeed, changing in three dimensions the symmetry of the two-body resonance from an $s$- to a $p$-wave Efremov et al. (2013); Zhu and Tan (2013) results in the reduction of the infinite number of bound states to a finite one.

Moreover, in the case of a two- Bruch and Tjon (1979); Lim and Maurone (1980); Vugal’ter and Zhishin (1983); Levinsen et al. (2014) or one-dimensional Kartavtsev et al. (2009); Mehta (2014) space, a two-body $s$-wave resonance does not lead to the Efimov effect. Again the spectrum of the three-body bound states is finite and determined by the mass ratios between the particles Pricoupenko and Pedri (2010); Bellotti et al. (2013); Ngampruetikorn et al. (2013). However, the two-dimensional system of three particles with a $p$-wave inter-particle resonance can again support an infinite number of universal bound states, the so-called “super Efimov” effect Nishida et al. (2013); Moroz and Nishida (2014); Gridnev (2014); Volosniev et al. (2014); Gao et al. (2015).

Experimentally the changes in the number of space dimensions and the interaction can be implemented. Indeed, the reduction of the dimensionality is achieved by using off-resonant light to confine ultra-cold gases in quasi-1D or quasi-2D geometries Bloch et al. (2008). In addition, the interactions between ultracold atoms can be tuned easily via Feshbach-resonances Chin et al. (2010).

I.2 Methods

We solve the exact integral equations Skorniakov and Ter-Martirosian (1957) of Skorniakov and Ter-Martirosian (STM) for the zero-range heavy-light interaction potential and obtain the three-body bound states for arbitrary mass ratios. Based on these exact results, we investigate the accuracy of the Born-Oppenheimer (BO) approximation Born and Oppenheimer (1927) for the three-body problem depending on the mass ratio.

By considering finite-range potentials of Gaussian and cubic Lorentzian shape, we explore the universal regime. For these finite-range potentials we obtain the bound states of the three-particle system numerically using a pseudospectral method Boyd (2000); Trefethen (2000); Baye (2015) based on the roots of the rational Chebyshev functions.

I.3 Overview

Our article is organized as follows. In Section II we briefly summarize the essential ingredients of the two- and three-body system. We then focus in Section III on the case of the zero-range heavy-light interaction and utilize the BO approximation and the STM method. Next we dedicate Section IV to a study of the universal behavior for two different finite-range potentials. In Section V we then demonstrate the universality of the three-body bound states for any heavy-light interaction. We conclude in Section VI by summarizing our results and by presenting an outlook.

In order to keep our article self-contained but focused on the central ideas, we present more detailed calculations in two appendices. Appendix A is focused on the derivation of the diagonal correction to the BO approximation. In Appendix B, we introduce a grid based on the roots of the rational Chebyshev functions and recall briefly the pseudospectral method applied in Section IV.

II The three-body system

In this section we first briefly discuss the validity of 1D models to describe quasi-1D systems. We then introduce the quantities determining an interacting mass-imbalanced two-body system in 1D. Next, we extend this system to the case of three interacting particles using dimensionless Jacobi coordinates. Finally, we discuss the corresponding Schrödinger equation and its symmetries, which serves as the basis for the studies presented in the subsequent sections.

II.1 1D and quasi-1D models

Many theoretical studies Mehta et al. (2007); Kartavtsev et al. (2009); Mehta (2014); Nishida (2018); Guijarro et al. (2018) of three-body systems confined along two directions are performed using 1D models. This reduction offers the advantage of a simple and intuitive description revealing the underlying three-body properties. However, it is important to emphasize that experiments on these confined systems are always performed in quasi-1D.

In the case of a zero-range interaction the effective interaction potential of two particles in a tight cylindrical symmetric trap (quasi-1D setup) is given by the zero-range potential with the 1D scattering length determined by the 3D scattering length and the harmonic potential width, as shown in Ref. Olshanii (1998). Moreover, the dependence of universal three-body bound states on the dimensionality has been investigated in Refs. Levinsen et al. (2014); Sandoval et al. (2015); Yamashita et al. (2018); Pricoupenko (2018). In particular, when reducing the dimensionality from 3D to quasi-2D, the conditions to reproduce the results obtained by a 2D model are presented.

These results justify the relevance of 1D models for quasi-1D experiments and we thus analyze in the present article the three-body system using a 1D model.

II.2 Two interacting particles

We consider a two-body system consisting of a heavy particle of mass $M$ and a light one of mass $m$, both constrained to 1D and interacting via a potential of range $\xi_{0}$.

After eliminating the heavy-light center-of-mass coordinate, the system is governed by the stationary Schrödinger equation

[TABLE]

for the two-body wave function $\psi^{(2)}=\psi^{(2)}(x)$ of the relative motion presented in dimensionless units. Indeed, $x$ denotes the relative coordinate of the light particle with respect to the heavy one in units of the characteristic length $\xi_{0}$.

The two-body binding energy $\mathcal{E}^{(2)}$ and the potential

[TABLE]

are both given in units of $\hbar^{2}/\mu\xi_{0}^{2}$ with the Planck’s constant $\hbar$ and the reduced mass $\mu\equiv Mm/(m+M)$ of the heavy-light system. Here $v_{0}$ denotes the magnitude and $f=f(x)$ the shape of the interaction potential.

We assume an attractive interaction, $v_{0}<0$, as well as a symmetric shape $f$, that is $f(x)=f(\left|{x}\right|)$. Moreover, we choose $v$ such that (i) it describes a short-range interaction, i.e. $\left|{x}\right|^{2}f(|x|)\to 0$ as $\left|{x}\right|\to\infty$, and (ii) the potential $v$ supports only a single bound state with energy $\mathcal{E}^{(2)}_{\mathrm{g}}$ and even wave function $\psi^{(2)}_{\mathrm{g}}(x)=\psi^{(2)}_{\mathrm{g}}(-x)$.

II.3 Three interacting particles

We now add a third particle of mass $M$, also constrained to 1D and identical to the heavy particle in the heavy-light system considered above. Accordingly, we assume the same interaction potential $v$ between the additional heavy and the light particle, but no interaction between the two heavy ones.

Next, we introduce dimensionless Jacobi coordinates Greene et al. (2017) $x$ and $y$ as displayed in Fig. 1, where $y$ is the relative coordinate between the two heavy particles, and $x$ denotes the coordinate of the light particle with respect to the center-of-mass $C$ of the two heavy ones, both in units of $\xi_{0}$.

Eliminating again the center-of-mass motion of this heavy-heavy-light system, we arrive at the dimensionless stationary Schrödinger equation

[TABLE]

for the three-body wave function $\psi=\psi(x,y)$ describing only the relative motions with $r_{\pm}\equiv x\pm y/2$.

The coefficients

[TABLE]

and

[TABLE]

depend only on the mass ratio $M/m$, and $\mathcal{E}$ denotes the dimensionless three-body energy in units of $\hbar^{2}/\mu\xi_{0}^{2}$.

We notice that Eq. (3) is invariant under the transformation $y\rightarrow-y$, that is an exchange of the two heavy particles. Hence, we distinguish even solutions $\psi(x,-y)=\psi(x,y)$ corresponding to two heavy bosonic particles, and odd solutions $\psi(x,-y)=-\psi(x,y)$ corresponding to two heavy fermionic particles. Moreover, also the transformation $x\rightarrow-x$ leaves Eq. (3) invariant, and leads to the additional symmetry $\psi(-x,y)=\pm\psi(x,y)$.

II.4 Formulation of the problem

Now our aim is to solve Eq. (3) for the three-body bound state with the wave function $\psi_{n}$ and the corresponding energy $\mathcal{E}_{n}$ for $n=0,1,\ldots$. In particular, we are interested in the case when the two-body interaction described by the potential $v$ is close to a resonance, that is the energy $\mathcal{E}_{\mathrm{g}}^{(2)}$ of the two-body ground state approaches zero. Under the assumptions on $v$ presented in Sec. II.2, we then expect a universal behavior of the spectrum $\mathcal{E}_{n}$, namely that in the limit $\mathcal{E}_{\mathrm{g}}^{(2)}\to 0$ the ratio

[TABLE]

is independent of the shape $f$ of the interaction potential.

In order to obtain a Hamiltonian with the eigenenergies $\epsilon_{n}$ given by Eq. (6), we introduce the rescaled variables

[TABLE]

with

[TABLE]

and rewrite Eq. (3) as the equation

[TABLE]

for the wave function $\tilde{\psi}=\tilde{\psi}(X,Y)$ with

[TABLE]

and

[TABLE]

We emphasize that the eigenvalues $\epsilon$ correspond to the ratio of the dimensional three-body and two-body energies and are hence accessible in an experiment.

III Contact interaction

We start our analysis by considering a contact interaction

[TABLE]

between the light particle and each heavy one, with $\delta(x)$ being the Dirac delta function. For the two-body problem, this interaction potential $v_{0}f_{\delta}(x)$ has only one bound state with the energy

[TABLE]

determined by the magnitude $v_{0}$ of the potential.

Using this relation and the scaling property of the Dirac delta function, $\delta(\alpha x)=\delta(x)/\left|{\alpha}\right|$, we obtain

[TABLE]

and the three-body Schrödinger equation, Eq. (9), becomes independent of the interaction strength $v_{0}$. Hence $\epsilon$, the three-body binding energy in units of the two-body ground state energy, does not depend on $v_{0}$.

We now solve Eq. (9) with $f=f_{\delta}$ using two different methods: the BO approximation Born and Oppenheimer (1927) and an approach based on the exact STM integral equation Skorniakov and Ter-Martirosian (1957). We then compare the results of the two techniques to quantify the error of the BO approximation.

III.1 Born-Oppenheimer approximation

The Born-Oppenheimer (BO) approach relies on approximating Efremov et al. (2009); Fonseca et al. (1979) the total three-body wave function $\tilde{\psi}$ in Eq. (9) by the product

[TABLE]

Here, the wave function $\varphi(X|Y)$ describes the dynamics of the light particle in the potential of the two heavy ones, which are assumed to stay at a fixed distance $Y$.

The physical motivation of the ansatz Eq. (15) is that for a large mass ratio, $M/m\gg 1$, the change of distance between the heavy particles is negligible on the relevant timescales of the light-particle dynamics. Hence, $Y$ does not change and enters in $\varphi$ only as a parameter, indicated by the vertical bar, giving rise to the Schrödinger equation

[TABLE]

for the wave function $\varphi$ of the light particle, determining the so-called BO potential $u=u(Y)$.

In Appendix A.1 we solve Eq. (16) analytically and obtain

[TABLE]

expressed in terms of the Lambert function $W_{0}$ Abramowitz and Stegun (1972), and the corresponding wave functions

[TABLE]

where $N_{\pm}$ is a normalization factor.

The two potentials $u_{\pm}=u_{\pm}(Y)$ are displayed in Fig. 2. Only the lower curve $u_{+}$, corresponding to the symmetric light-particle state $\varphi_{+}$, provides an attractive potential for the two heavy particles and therefore supports bound states of the total three-body system.

The wave function $\phi_{+}=\phi_{+}(Y)$ of the heavy particles then obeys the Schrödinger equation

[TABLE]

where $u_{+}$ indeed plays the role of a potential, and $\epsilon^{(\mathrm{BO})}$ is the scaled three-body energy within the BO approach.

Using the attractive potential $u_{+}$ given by Eq. (17), we calculate the values of $\epsilon_{n}^{(\mathrm{BO})}$ numerically to a precision of $10^{-6}$ applying a pseudospectral method based on the Chebyshev grid introduced in Appendix B. It is important to mention that the scaled three-body bound state energies $\epsilon_{n}^{(\mathrm{BO})}$ satisfy the inequality $\epsilon^{(\mathrm{BO})}<-1/\alpha_{x}$, where the upper bound is the value of the BO potential at infinity,

[TABLE]

as shown in Fig. 2.

The number $n_{\mathrm{max}}$ of bound states supported by $u_{+}$ depends on the mass ratio $M/m$ and is depicted in Fig. 3 as a blue line together with the semiclassical Schleich (2001) estimation

[TABLE]

or

[TABLE]

depicted by an orange line.

With increasing mass ratio $M/m$, additional bound states appear. A detailed comparison of the critical mass ratios required for the formation of a new bound state within the BO approximation and an hyperspherical approach can be found in Refs. Kartavtsev et al. (2009); Mehta (2014).

III.2 Integral equation of Skorniakov and Ter-Martirosian

In this section we apply the method Skorniakov and Ter-Martirosian (1957) of Skorniakov and Ter-Martirosian (STM) to the three-body problem described by Eqs. (9) and (10) with a contact interaction. In contrast to the BO approach, this method does not involve any approximation, and in principle provides an exact solution for any mass ratio $M/m$.

We introduce the Green function

[TABLE]

for the two-dimensional free-particle Schrödinger equation with $\epsilon<0$, where $K_{0}$ denotes the modified Bessel function of the second kind Abramowitz and Stegun (1972). We can then cast Eqs. (9) and (10) in integral form

[TABLE]

where $R_{\pm}^{\prime}\equiv X^{\prime}\pm Y^{\prime}/2$.

With the help of Eq. (14) this expression simplifies in the special case of the contact interaction $f=f_{\delta}$ to

[TABLE]

The delta functions then allow us to immediately perform the integration over $Y^{\prime}$ and to obtain the one-dimensional integral equation

[TABLE]

Since the Hamiltonian $\hat{H}$, defined by Eq. (10), is invariant under the transformation $Y\to-Y$, that is under the exchange of the two heavy particles, the solution $\tilde{\psi}$ has to be either even or odd with the symmetry relations

[TABLE]

corresponding to the case of bosonic (plus sign) or fermionic (minus sign) heavy particles.

Evaluating both sides of Eq. (III.2) at $Y=2X$, we arrive at the integral equation

[TABLE]

with the kernel

[TABLE]

where we have used the symmetry relation given by Eq. (27) in writing $\tilde{\psi}(X^{\prime},-2X^{\prime})=\pm\tilde{\psi}(X^{\prime},2X^{\prime})$.

By rescaling the coordinates $X$ and $X^{\prime}$ by $\sqrt{\left|{\epsilon}\right|}$, Eq. (28) can be cast into an eigenvalue problem for the eigenfunction $\tilde{\psi}(X/\sqrt{\left|{\epsilon}\right|},2X/\sqrt{\left|{\epsilon}\right|})$ with eigenvalue $\sqrt{\left|{\epsilon}\right|}$, where the condition $\epsilon<-1$ determines the three-body bound states. Hence, the desired spectrum $\epsilon_{n}^{\star}$ of three-body bound states in units of the two-body ground state energy $\mathcal{E}_{\mathrm{g}}^{(2)}$ can be efficiently computed.

The three-body wave function $\tilde{\psi}(X,Y)$ is more difficult to obtain and requires an additional step. Together with the spectrum $\epsilon_{n}^{\star}$, we first obtain $\tilde{\psi}_{n}^{\star}(X,2X)$ from Eq. (28), that is $\tilde{\psi}_{n}^{\star}$ along the lines of interaction $Y=\pm 2X$. Then, we insert both $\epsilon_{n}^{\star}$ and $\tilde{\psi}_{n}^{\star}(X,2X)$ into the right-hand side of Eq. (III.2). Taking into account the symmetry property $\tilde{\psi}_{n}^{\star}(X,-2X)=(-1)^{n}\tilde{\psi}_{n}^{\star}(X,2X)$ (even $n$ correspond to bosonic heavy particles, whereas odd $n$ represent the fermionic case) and performing the integration over $X^{\prime}$, we finally obtain the entire three-body wave function $\tilde{\psi}_{n}^{\star}=\tilde{\psi}_{n}^{\star}(X,Y)$.

In Fig. 4 we depict the four ($n=0,1,2,3$) lowest three-body bound states obtained via the STM method for $M/m=20$. We emphasize again the scaling property of the delta function yielding the wave function

[TABLE]

in the unscaled variables $x$ and $y$.

III.3 BO approximation vs. STM approach

In the preceding subsections we have applied the BO and STM methods to solve the 1D three-body problem with contact interaction. Now we compare the dependence of the resulting spectra and wave functions on the mass ratio $M/m$. Common experimental mass ratios range from $M/m=1$ for identical particles via $M/m\cong 2.2$ and $M/m\cong 12.4$ for ${}^{87}\mathrm{Rb}$–${}^{40}\mathrm{K}$ and ${}^{87}\mathrm{Rb}$–${}^{7}\mathrm{Li}$ mixtures respectively, to more extreme values of $M/m\cong 22.2$ in case of ${}^{133}\mathrm{Cs}$–${}^{6}\mathrm{Li}$ mixtures Naidon and Endo (2017). Therefore we choose this range for our analysis.

III.3.1 Energy spectrum

In Fig. 5 we display the three-body ground state energy $\epsilon_{0}$ obtained by the BO (blue dots) and the STM (yellow diamonds) method as a function of the mass ratio $M/m$. In addition, the relative error

[TABLE]

is depicted as black squares for $n=0$.

For both methods the energies are computed numerically and are accurate up to $10^{-6}$. We find that our results are in excellent agreement with the values in the literature. Indeed, for $M/m=1$ we obtain $\epsilon_{0}^{(\mathrm{BO})}=-2.42267$ being within accuracy of the value $-2.4227$ in Ref. Mehta (2014), as well as $\epsilon_{0}^{\star}=-2.087719$ which matches the value $-2.087719$ and is very close to $-2.08754$ found in Refs. Kartavtsev et al. (2009) and Gaudin and Derrida (1975), respectively.

Likewise, we obtain the previously reported Mehta (2014) error of about $16\,\%$ for the BO ground state energy $\epsilon_{0}^{(\mathrm{BO})}$ for $M/m=1$, a regime in which the BO approximation is not expected to provide reasonable results. Moreover, the relative error $\delta\epsilon_{0}$ decreases monotonically with increasing mass ratio $M/m$ and drops below $2\,\%$ for $M/m\cong 22.2$, as shown in Fig. 5.

We depict in the top row of Fig. 6 the relative errors $\delta\epsilon_{n}$ as a function of the mass ratio $M/m$ for the lowest four bound states $(n=0,1,2,3)$. Excited states ($n\geq 1$) appear with increasing mass ratio, as shown in Fig. 3. The higher excited a state is, the larger the corresponding error gets. This behavior can be understood from the fact that the BO approximation involves neglecting derivatives, effecting more strongly higher excited states as they are more oscillatory.

III.3.2 Wave functions

After comparing the energy spectra calculated within the BO and the STM method, we now turn to the wave functions obtained by both methods and study the fidelity Nielsen and Chuang (2000) which for pure states simplifies to the spatial overlap

[TABLE]

of the wave functions $\tilde{\psi}_{n}^{(\mathrm{BO})}$ and $\tilde{\psi}_{n}^{\star}$. The BO wave function $\tilde{\psi}_{n}^{(\mathrm{BO})}$ is a product of $\varphi_{+}$, Eq. (18), and $\phi_{+n}$, obtained from Eq. (19), with $u_{+}$ given by Eq. (17), whereas the STM wave function $\tilde{\psi}_{n}^{\star}$, Eq. (30), is calculated directly from Eq. (III.2). Both functions are evaluated numerically on the same Chebyshev grid introduced in Appendix B, and the computed fidelities are accurate up to $10^{-4}$.

In the bottom row of Fig. 6 we present $\mathcal{F}_{n}$ for the four lowest $(n=0,1,2,3)$ bound states as a function of $M/m$. As expected, the fidelity increases monotonically for all bound states with increasing mass ratio $M/m$, that is the BO approximation becomes more accurate. However, it is remarkable that the fidelity $\mathcal{F}_{n}$ starts already from 0.988 at $M/m=1$ and reaches values up to 0.999 for $M/m=25$.

Moreover, $\mathcal{F}_{n}$ does not show a clear dependence on $n$: higher excited states do not always have a lower fidelity, in contrast to the expectation that the BO approximation should be worse for higher excited states. It is nevertheless possible that this behavior arises for much larger values of the ratio $M/m$.

III.3.3 Diagonal energy correction

We emphasize that the comparisons of the energy spectra and of the wave functions are based on different measures. The fidelity indicates how well a state can mimic another one in a measurement. As the fidelity is almost unity, we expect only a small deviation in the spectrum with respect to $\epsilon_{n}^{\star}$, if the exact state $\tilde{\psi}_{n}^{\star}$ is replaced by $\tilde{\psi}_{n}^{(\mathrm{BO})}$ leading to the expression

[TABLE]

for the mean value of the energy. Here $\hat{H}$ is the full three-body Hamiltonian defined in Eq. (10) with $f=f_{\delta}$.

As shown in Appendix A.2, $\bar{\epsilon}_{n}^{(\mathrm{BO})}$ coincides with the BO energies including diagonal correction terms, and is depicted in Fig. 7 for the two lowest states ($n=0,1$). Compared to the zero-order BO approximation $\epsilon_{n}^{(\mathrm{BO})}$, the deviation with respect to $\epsilon_{n}^{\star}$ is reduced by up to an order of magnitude. Hence, this feature suggests that the major contribution to the deviation between $\epsilon_{n}^{(\mathrm{BO})}$ and $\epsilon_{n}^{\star}$ stems from the Hamiltonian itself, and not from the wave functions.

III.3.4 Summary

In summary, for the contact interaction, the BO approximation works surprisingly well in estimating the bound state energies, and even better for the corresponding wave functions. Moreover, for the contact interaction, the accuracy of the BO approximation is determined solely by the mass ratio between heavy and light particles and provides reasonable results even in the case of equal masses.

IV General interaction potentials

So far we have only studied the case of a contact interaction between heavy and light particles. In this section we focus on different short-range interaction potentials and apply a pseudospectral method based on the roots of rational Chebyshev functions.

In particular, we consider a two-body system close to a resonance and analyze the emergence of the universality in the mass-imbalanced three-body system. In this regime we retrieve for both the energy spectrum and the corresponding wave functions the results obtained for the case of a contact interaction.

IV.1 Two-body interaction

In this section we find numerically the relation between the two-body binding energy $\mathcal{E}^{(2)}_{\mathrm{g}}$ of the ground state, and the potential depth $v_{0}$ for different shapes $f$ of the interaction potential, Eq. (2). For this purpose, we apply a pseudospectral method Boyd (2000); Trefethen (2000); Baye (2015) using a grid based on the roots of the rational Chebyshev functions Boyd (1987).

According to Appendix B, we represent the dimensionless Schrödinger equation, Eq. (1), for the two-body system as a generalized eigenvalue problem

[TABLE]

with the generalized eigenvalue $-v_{0}$ and the generalized eigenvector $\vec{\psi}_{\mathrm{g}}^{(2)}$ of size $\mathcal{N}$ containing the values of the function $\psi_{\mathrm{g}}^{(2)}(x)$ evaluated at the grid points. The matrices $\boldsymbol{D_{2}}$ and $\boldsymbol{F}$ of size $\mathcal{N}\times\mathcal{N}$ are introduced in Appendix B.1, and $\boldsymbol{\mathds{1}}$ denotes the identity matrix.

For a given two-body binding energy $\mathcal{E}^{(2)}_{\mathrm{g}}$ we determine the lowest generalized eigenvalue $-v_{0}=\left|v_{0}\right|$ of Eq. (34) which specifies the potential depth, as well as the corresponding generalized eigenvector $\vec{\psi}_{\mathrm{g}}^{(2)}$, yielding an approximation to the wave function $\psi_{\mathrm{g}}^{(2)}(x)$ of the lowest state with the energy $\mathcal{E}^{(2)}_{\mathrm{g}}$.

We perform this calculation for two different interaction potentials $v=v(x)$, namely a potential

[TABLE]

with a Gaussian shape and a potential

[TABLE]

being characterized by the cube of a Lorentzian.

In order to reach sufficient convergence, we use $\mathcal{N}=2500$ and numerically obtain from Eq. (34) the potential depth $\left|{v_{0}}\right|$ as a function of $\left|{\mathcal{E}_{\mathrm{g}}^{(2)}}\right|$, displayed in Fig. 8 by empty blue and filled red circles corresponding to $f_{\mathrm{G}}$ and $f_{\mathrm{L}}$.

In the limit $\left|v_{0}\right|\to 0$, the binding energy of the ground state is approximated by the expression Simon (1976); Berger et al. (2008)

[TABLE]

From Eq. (37) we obtain for our two test potentials the approximation

[TABLE]

for $f=f_{\textrm{G}}$, and

[TABLE]

for $f=f_{\textrm{L}}$, depicted by a dashed blue and red line in Fig. 8, accordingly.

In the case of a contact interaction with $f=f_{\delta}$, Eq. (12), the relation given by Eq. (37) not only provides an approximation, but is exact as presented in Eq. (13). The pure quadratic dependence of the two-body binding energy on the potential depth $v_{0}$, and the fact that the contact interaction gives rise to only a single bound state for any value $v_{0}<0$ confirms that the corresponding two-body system is exactly on resonance. In the next sections we show that this unique feature has important consequences for the respective three-body system.

IV.2 Universal limit

We now consider the three-body problem in 1D with the heavy-light interaction potentials having the shape $f$ given by Eqs. (35) and (36), and compare the results to those obtained for the contact interaction, Eq. (12). In particular, we study the universal limit, that is $\mathcal{E}_{\mathrm{g}}^{(2)}\to 0$.

Following Appendix B, we apply a pseudospectral method based on the roots of the rational Chebyshev functions and represent the Schrödinger equation, Eq. (3), as the eigenvalue problem

[TABLE]

for the eigenvalue $\mathcal{E}$ and the eigenvector $\vec{\psi}$ of size $\mathcal{M}$ approximating the three-body wave function $\psi(x,y)$. Here, the matrices $\boldsymbol{D_{xx}}$ and $\boldsymbol{D_{yy}}$ of size $\mathcal{M}\times\mathcal{M}$ correspond to the discretized second-order derivative with respect to $x$ and $y$, respectively. Moreover, the diagonal matrices $v_{0}\boldsymbol{F_{+}}$ and $v_{0}\boldsymbol{F_{-}}$ describe accordingly the interaction potentials along the lines $x+y/2$ and $x-y/2$, as shown in Appendix B.2.

We solve the finite-dimensional eigenvalue problem, Eq. (40), on the Data Vortex system DV206 Dat by employing a parallelized version of the ARPACK software ARP , including an implementation of the Implicitly Restarted Arnoldi Method Lehoucq and Sorensen (1996). In order to obtain sufficient convergence of all energies and the corresponding wave functions, we use a grid of size $\mathcal{M}=\mathcal{N}_{x}\cdot\mathcal{N}_{y}$ with $\mathcal{N}_{x}=512$ and $\mathcal{N}_{y}=256$.

IV.2.1 Energy spectrum

For the mass ratio $M/m=20$ we present in Fig. 9 (a) and (b) the scaled energies $\epsilon_{n}$, Eq. (6), of the first four three-body bound states ($n=0,1,2,3$), as a function of the two-body binding energy $\left|\mathcal{E}_{\mathrm{g}}^{(2)}\right|$. In particular, we display by empty blue and filled red symbols the results for a Gaussian-shaped potential $f_{\mathrm{G}}$, Eq. (35), and for a cubic Lorentzian-shaped potential $f_{\mathrm{L}}$, Eq. (36). The values of $\epsilon_{n}$ corresponding to the contact interaction $f_{\delta}$, Eq. (12), are independent of $\left|\mathcal{E}_{\mathrm{g}}^{(2)}\right|$ and are shown by gray lines, which reflects the feature of this interaction potential to support an exact two-body resonance.

In Fig. 9 we separate the cases of bosonic heavy particles (a) associated with $n=0,2$ and fermionic heavy particles (b) represented by $n=1,3$. For the interaction potentials with the shapes $f_{\mathrm{G}}$ and $f_{\mathrm{L}}$, we have used the numerically obtained relation between the two-body binding energy $\mathcal{E}_{\mathrm{g}}^{(2)}$ and the potential depth $v_{0}$ presented in Fig. 8.

In the case that the two-body heavy-light subsystem is close to a resonance, that is in the limit $\mathcal{E}_{\mathrm{g}}^{(2)}\to 0$, we observe a universal behavior of the scaled energies $\epsilon_{n}$ for all presented three-body bound states, that is for $n=0,1,2,3$. Moreover, we point out that different states approach the universal regime in different ways. Indeed, Fig. 9 (a) and (b) show clearly that the difference of the energy $\epsilon_{n}$ and the corresponding universal limit for a fixed value of the two-body binding energy $\mathcal{E}_{\mathrm{g}}^{(2)}$ is usually smaller in the case of higher excited states.

IV.2.2 Wave functions

Now we are in the position to compare not only the energies of the three-body bound states for different interaction potentials, but also the corresponding wave functions. For this purpose, we use again the fidelity

[TABLE]

as a measure of the spatial overlap between the wave function $\psi_{n}^{\star}$, Eq. (30), obtained for the case of a contact interaction $f_{\delta}$, Eq. (12), and the three-body wave function $\psi_{n}$ obtained as solution of Eq. (40) for the interaction potential with the shape $f_{\mathrm{G}}$ and $f_{\mathrm{L}}$, respectively.

Using the relation between the potential depth $v_{0}$ and the two-body binding energy $\mathcal{E}_{\mathrm{g}}^{(2)}$ presented in Fig. 8 we show in Fig. 9 the fidelities $\mathfrak{F}_{n}$ for (c) bosonic and (d) fermionic heavy particles as a function of the two-body binding energy $\mathcal{E}_{\mathrm{g}}^{(2)}$ by empty blue for $f_{\mathrm{G}}$, and filled red symbols for $f_{\mathrm{L}}$. The value $\mathfrak{F}_{n}=1$ obtained for a contact interaction for any value $\left|{\mathcal{E}_{\mathrm{g}}^{(2)}}\right|$ is displayed by gray lines. Independent of the shape $f$ of the interaction potential, the fidelity $\mathfrak{F}_{n}$ for $n=0,1,2,3$ approaches unity as $\mathcal{E}_{\mathrm{g}}^{(2)}\to 0$, and thus describes a perfect overlap of $\psi_{n}$ and $\psi_{n}^{\star}$.

As shown by Eq. (30) the coordinates of the wave function $\psi_{n}^{\star}$ are simply rescaled as $\mathcal{E}_{\mathrm{g}}^{(2)}$ is varied due to the scaling property, Eq. (14), for $f=f_{\delta}$. Thus, close to the two-body resonance, $\mathcal{E}_{\mathrm{g}}^{(2)}\to 0$, the same behavior is true for the wave function $\psi_{n}$, revealing again the universal limit.

IV.2.3 Summary

In summary, the universal behavior of the three-body system is shown to appear both in the scaled energies $\epsilon_{n}$, Eq. (6), as well as in the scaling of the corresponding wave functions $\psi_{n}$. Moreover, we emphasize that the universal behavior manifests itself for all presented three-body bound states.

V Proof of universality

In the preceding section we have explored the universal behavior of the three-body bound state energies and the corresponding wave functions for interaction potentials of Gaussian and cubic Lorentzian shape, with the two-body ground state energy $\mathcal{E}_{\mathrm{g}}^{(2)}$ approaching zero. In this limit, we now prove the universality of the 1D three-body system for an arbitrary mass ratio $M/m$ and any short-range interaction potential.

For this purpose we consider a heavy-light interaction of shape $f$ and recall the three-body Schrödinger equation in integral form, given by Eq. (III.2). Next, we perform the substitutions $Z\equiv\ell_{\mathrm{g}}R_{+}^{\prime}=\ell_{\mathrm{g}}(X^{\prime}+Y^{\prime}/2)$ on the first summand and $Z\equiv\ell_{\mathrm{g}}R_{-}^{\prime}=\ell_{\mathrm{g}}(X^{\prime}-Y^{\prime}/2)$ on the second one, to arrive at

[TABLE]

Here, we have used $\ell_{\mathrm{g}}$ defined by Eq. (8). With the approximate expression, Eq. (37), for the two-body binding energy, valid in the regime $\mathcal{E}_{\mathrm{g}}^{(2)}\ll 1$, we obtain

[TABLE]

In the limit $\mathcal{E}_{\mathrm{g}}^{(2)}\to 0$, that is $\ell_{\mathrm{g}}\to\infty$, $G_{\epsilon}^{(2)}$ and $\tilde{\psi}$ become independent of $Z$, and as a result any dependence on the potential shape $f$ cancels in Eq. (V). Indeed, we retrieve Eq. (III.2) valid for the contact interaction, with the solutions $\epsilon_{n}=\epsilon_{n}^{\star}$ and $\tilde{\psi}_{n}(X,Y)=\tilde{\psi}_{n}^{\star}(X,Y)$, as considered in Section III.2. We emphasize that this is a consequence of the fact that Eq. (37) is exact for this particular interaction potential.

As a result, these universal constants $\epsilon_{n}^{\star}=\epsilon_{n}^{\star}(M/m)$ depend only on the mass ratio and can be used to formulate the relation

[TABLE]

for the three-body binding energies as a function of the two-body interaction, valid for $\left|{v_{0}}\right|\to 0$.

Hence, we have shown explicitly that all scaled energies $\epsilon_{n}$, as well as the wave functions $\tilde{\psi}_{n}$ coincide with the results for the contact interaction, for any short-range heavy-light interaction potential of shape $f$, provided we approach the two-body resonance defined by $\mathcal{E}_{\mathrm{g}}^{(2)}\to 0$.

VI Conclusion and outlook

In this article, we have presented a quantum mechanical treatment of a heavy-heavy-light system confined to 1D. For a zero-range heavy-light interaction we have studied the three-body energy spectrum and the corresponding wave functions using two different methods: (i) the Born-Oppenheimer approximation, and (ii) the exact integral equations of Skorniakov and Ter-Martirosian. In addition, for finite-range interactions, we have investigated the universal limit of the three-body energies and the corresponding wave functions when the ground state energy of the heavy-light subsystem approaches zero.

In particular, for the case of a contact interaction we have explored the accuracy of the BO approximation in a regime of experimentally feasible mass ratios and found that the error in the energy spectrum drops rapidly from around 20% in case of equal masses to below 2% for rather extreme mass ratios $M/m\cong 22.2$ like in ${}^{133}\mathrm{Cs}$–${}^{6}\mathrm{Li}$ mixtures Naidon and Endo (2017). In addition, the ground state energy presented in Ref. Mehta (2014) for $M/m=1$ agrees with our result.

The approximate BO wave functions are very close to the exact ones, since for $M/m=25$ the fidelity reaches values up to $0.999$. As a result, the use of the approximate BO wave functions to calculate the mean value of the total Hamiltonian has significantly improved the accuracy of the three-body binding energies.

Moreover, by applying a pseudospectral method based on the roots of rational Chebyshev functions we have obtained the three-body energies and wave functions for the short-range interaction potentials of Gaussian and cubic Lorentzian shape. When the ground state energy of the heavy-light potential approaches zero, the universal behavior is apparent for both potentials, that is each three-body binding energy converges to the limit value determined by the zero-range contact interaction. We have also compared the associated wave functions to the ones provided by the contact interaction and we found that they follow a universal scaling law when the two-body resonance is approached.

Finally, we have demonstrated the universality of all three-body bound states for any short-range interaction potential when the heavy-light ground state energy approaches zero. Here, we recover the results for the contact interaction obtained within the STM approach. Hence, the three-body bound states for an arbitrary short-range interaction on resonance can be obtained by using a zero-range potential and applying the BO approximation, provided the mass ratio is sufficiently large. For experimentally relevant mass ratios, we present in Tab. 1 the universal constants determining the three-body binding energies in case of weak interactions.

We conclude by raising a few interesting generalizations of our approach. A nearly resonant excited state in the two-body system might lead to different features compared to the ones induced by the two-body ground state. According to Ref. Barlette et al. (2000) the two-body scattering in 1D only depends on the symmetry of the state. Based on this argument one might conclude that universal behavior in the three-body system only depends on the symmetry of the underlying two-body resonance. However, a rigorous study of this case is necessary to arrive at a definite statement. Moreover, we emphasize that further features might appear within a 3D consideration for the quasi-1D three-body system. Needless to say these questions go beyond the scope of the present article but will be addressed in a future publication.

Acknowledgements.

We thank N.L. Harshman for fruitful discussions. Moreover, we gratefully acknowledge the support of Data Vortex Technologies which provided the Data Vortex system used for the numerical calculations. L.H. and M.A.E. thank the Center for Integrated Quantum Science and Technology (IQST) for financial support. This work is funded in part by the German-Israeli Project Cooperation DIP (Project No. AR 924/1-1, DU 1086/2-1). W.P.S. is most grateful to Texas A$\&$M University for a Faculty Fellowship at the Hagler Institute for Advanced Study at the Texas A$\&$M University as well as to Texas A$\&$M AgriLife Research. The research of the IQST is financially supported by the Ministry of Science, Research and Arts Baden-Württemberg.

Appendix A Born-Oppenheimer approximation for the three-body problem

In this appendix we recall the main ideas of applying Efremov et al. (2009); Fonseca et al. (1979) the BO approximation to the three-body problem in the presence of a contact interaction. In particular, we derive the relevant formulas in zero order as well as the diagonal correction to the energy spectrum.

In the BO approach we represent the total wave function as the product

[TABLE]

where we assign $\varphi$ and $\phi$ to the dynamics of the light and heavy particles, respectively.

When we apply the complete three-body Hamiltonian

[TABLE]

with

[TABLE]

onto $\tilde{\psi}$, we obtain

[TABLE]

So far our calculation is exact.

A.1 Zero-order consideration

In the zero-order approximation we neglect all derivatives of $\varphi$ with respect to the relative coordinate $Y$ of the heavy particles. This fact is emphasized by the vertical bar in the notation $\varphi(X,Y)\to\varphi(X|Y)$ and suggests to choose $\left\{\varphi\right\}$ as eigenbasis of $\hat{H}_{0}$, summarized by the light-particle Schrödinger equation

[TABLE]

We obtain the BO potential $u=u(Y)$ by rewriting Eq. (49) with $\hat{H}_{0}$ given by Eq. (47) in integral form

[TABLE]

where

[TABLE]

is the Green function of the one-dimensional free-particle Schrödinger equation for $u<0$.

Due to the delta functions, the integration over $X^{\prime}$ can be performed immediately. By evaluating both sides of Eq. (A.1) at the points $X=\pm Y/2$, we arrive at the transcendental equations

[TABLE]

for the BO potentials $u_{\pm}=u_{\pm}(Y)$ with the solutions

[TABLE]

in terms of the Lambert function $W_{0}$ Abramowitz and Stegun (1972), and the corresponding wave functions

[TABLE]

where

[TABLE]

is a normalization factor.

The wave function $\phi_{\pm}=\phi_{\pm}(Y)$ is then a solution of the Schrödinger equation

[TABLE]

where the potential $u_{\pm}$ is given by Eq. (53), and $\epsilon^{(\mathrm{BO})}$ denotes the three-body binding energy in the zero-order BO approximation.

A.2 Diagonal correction to the energy spectrum

In the zero-order BO approximation, we neglect the last two terms in Eq. (A). However, in order to find corrections to these zero-order expressions, we have to consider now all terms in Eq. (A). In this section, we derive the diagonal correction to the BO binding energies and find the connection to the mean value $\bar{\epsilon}^{(\mathrm{BO})}$ defined in Eq. (33)

To distinguish between the so-called diagonal and non-diagonal contributions, we consider Eq. (A) with the zero-order solutions $\varphi=\varphi_{i}$, $\phi=\phi_{in}$, where the subscript $i$ labels the light-particle channels, and $n$ numbers the state in each channel. We then multiply Eq. (A) by $\varphi^{\ast}_{k}=\varphi_{k}$ from the left-hand side, perform the integration over $X$, and use the orthonormality of the light-particle states $\varphi_{i}$ ($i,k=\pm$) to write

[TABLE]

where $\delta_{ik}$ is the Kronecker delta and

[TABLE]

Here we have used Eqs. (49) and (56) in order to identify the zero-order contribution $\epsilon^{(\mathrm{BO})}$. For $k=i$ ($k\neq i$) we speak of the diagonal (non-diagonal) part.

In the diagonal case, the expression for $\delta H_{ii}$ simplifies, as the second term in Eq. (58) vanishes

[TABLE]

due to normalization.

In order to calculate the mean value

[TABLE]

defined by Eq. (33), we multiply Eq. (A) by $\tilde{\psi}_{n}^{(\mathrm{BO})}\equiv\varphi_{+}\phi_{+n}$ from the left, integrate over $X$, and use Eq. (57) for $k=i$ and Eq. (59), to arrive at

[TABLE]

For the mean value, Eq. (60), only the diagonal term ($k=i$) of Eq. (58) contributes to the correction. Hence, the expression $\bar{\epsilon}_{n}^{(\mathrm{BO})}$ given by Eq. (61) equals the corrected BO binding energies, if couplings between different states are neglected, that is if $\delta\hat{H}_{ik}=0$ for $k\neq i$.

Appendix B Pseudospectral methods

Pseudospectral methods Boyd (2000); Trefethen (2000); Baye (2015) are an efficient tool to obtain a numerical solution of an ordinary or partial differential equation. In the following, we focus only on linear equations, where we represent the differential operators by matrices, and the unknown eigenfunctions by vectors. The corresponding eigenvalue problem of finite size can then be solved numerically.

Indeed, a key advantage of pseudospectral methods is the exponential convergence of the approximate solution to the exact one as the matrix size increases. For problems on a finite domain, the convergence rate is usually geometric, whereas convergence for problems on an infinite domain Boyd (2000) is usually subgeometric. However, for a given matrix size the accuracy of the approximate solution is crucially determined by the deployed set of basis functions. Throughout this article we follow the suggestion of Boyd Boyd (2000, 1987) and choose the rational Chebyshev functions as a basis.

In Appendix B.1 we present the matrices used for a finite dimensional representation of a linear ordinary differential equation of second order. A generalization of these matrices is obtained in Appendix B.2 for the case of a linear partial differential equation depending on two variables. Finally, we consider in Appendix B.3 the discretization of the eigenvalue problems analyzed in this article.

B.1 Matrix representation of 1D-problems

We begin by reviewing matrix representations of differential operators defined on the finite domain $(-1,1)$, where Chebyshev polynomials Boyd (2000) are used as basis functions. Next, we apply an algebraic map Boyd (1987) and obtain a finite dimensional representation of these operators on the complete real domain.

B.1.1 Finite domain

First, we consider a grid based on the roots

[TABLE]

of the Chebyshev polynomial $T_{\mathcal{N}}=T_{\mathcal{N}}(\eta)$ of the first kind with degree $\mathcal{N}$ and $i=0,\ldots,\mathcal{N}-1$.

This polynomial is defined by the recurrence relation

[TABLE]

for $\mathcal{N}>1$ with $T_{0}(\eta)\equiv 1$ and $T_{1}(\eta)\equiv\eta$, where the argument $\eta$ is restricted to the finite interval $(-1,1)$.

For this grid the first-order derivative is represented Boyd (2000) by the matrix

[TABLE]

with $i,j=0,1,\ldots,\mathcal{N}-1$.

Similarly, the matrix representation of the second-order derivative reads

[TABLE]

B.1.2 Infinite domain

Next, we consider the variable $x\in(-\infty,\infty)$ and extend the previous grid to an infinite domain. For this purpose we introduce the new grid points

[TABLE]

obtained from the old ones $\eta_{i}$ given by Eq. (62) by applying an algebraic map Boyd (1987, 2000). The mapping parameter $L$ determines the effective size of the grid.

The grid points $x_{i}$ are the roots of the rational Chebyshev functions

[TABLE]

given in terms of the Chebyshev polynomials $T_{\mathcal{N}}$ of degree $\mathcal{N}$ as defined by Eq. (63).

As a result, the discrete representation of the differential operator $\frac{\textrm{d}}{\textrm{d}x}$ is given by the matrix

[TABLE]

of size $\mathcal{N}\times\mathcal{N}$ where the elements of the diagonal matrix $\boldsymbol{A}$ read

[TABLE]

Similarly, the discrete representation

[TABLE]

of the second-order differential operator $\frac{\textrm{d}^{2}}{\textrm{d}x^{2}}$ is also determined by the diagonal matrix $\boldsymbol{B}$ with elements

[TABLE]

Additionally, a matrix representation of any function $f=f(x)$ is given by the diagonal matrix $\boldsymbol{F}$ of size $\mathcal{N}\times\mathcal{N}$ with the elements

[TABLE]

obtained by evaluating $f$ at the grid points $x_{i}$.

B.2 Matrix representation of 2D-problems

Now, we generalize our grid to accommodate a partial differential equation depending on the two independent variables $x$ and $y$.

We introduce the grid points $(x_{i},y_{j})$ with

[TABLE]

and

[TABLE]

where $\eta_{x,i}\equiv\eta_{i}$ for $i=0,1,\ldots,\mathcal{N}_{x}-1$ and $\eta_{y,j}\equiv\eta_{j}$ for $j=0,1,\ldots,\mathcal{N}_{y}-1$. The grid points $\eta_{i}$ and $\eta_{j}$ are defined by Eq. (62) and the integers $\mathcal{N}_{x}$ and $\mathcal{N}_{y}$ denote the number of grid points used for the variables $x$ and $y$, respectively, with the corresponding mapping parameters $L_{x}$ and $L_{y}$.

The discrete representation of the partial second-order derivative $\frac{\partial^{2}}{\partial x^{2}}$ reads

[TABLE]

Here, the matrix $\boldsymbol{D_{2,x}}\equiv\boldsymbol{D_{2}}$ of the size $\mathcal{N}_{x}\times\mathcal{N}_{x}$ is given by Eq. (70), whereas $\boldsymbol{\mathds{1}_{y}}$ denotes the identity matrix of size $\mathcal{N}_{y}\times\mathcal{N}_{y}$. Thus, the matrix $\boldsymbol{D_{xx}}$ has the size $\mathcal{M}\times\mathcal{M}$ with $\mathcal{M}\equiv\mathcal{N}_{x}\cdot\mathcal{N}_{y}$.

In a similar way, the partial second-order derivative $\frac{\partial^{2}}{\partial y^{2}}$ is represented by the $\mathcal{M}\times\mathcal{M}$ matrix

[TABLE]

where $\boldsymbol{\mathds{1}_{x}}$ is the identity matrix of size $\mathcal{N}_{x}\times\mathcal{N}_{x}$, and the matrix $\boldsymbol{D_{2,y}}\equiv\boldsymbol{D_{2}}$ of the size $\mathcal{N}_{y}\times\mathcal{N}_{y}$ is given by Eq. (70).

Moreover, similar to Eq. (72), a function depending on the variables $x$ and $y$ is represented by a diagonal matrix. In particular, the function $f(x\pm y/2)$ is given by the diagonal matrix $\boldsymbol{F_{\pm}}$ of size $\mathcal{M}\times\mathcal{M}$ with the elements

[TABLE]

B.3 Eigenvalue problem

Finally, we possess all ingredients to represent the stationary Schrödinger equation for the two-body, and the three-body system as an eigenvalue problem in terms of matrices provided by a pseudospectral method being determined by the roots of rational Chebyshev functions.

We start by discussing the two-body system described the 1D Schrödinger equation

[TABLE]

given by Eq. (1) where we have used the definition, Eq. (2), of the interaction potential.

Using the matrices $\boldsymbol{D_{2}}$ and $\boldsymbol{F}$ of size $\mathcal{N}\times\mathcal{N}$ defined by Eqs. (70) and (72), we arrive at the eigenvalue problem

[TABLE]

for the eigenvector

[TABLE]

where the grid points $x_{0},\ldots,x_{\mathcal{N}-1}$ are given by Eqs. (62) and (66).

Next, we consider the three-body system governed by the Schrödinger equation

[TABLE]

given by Eq. (3) with $r_{\pm}\equiv x\pm y/2$ where $\alpha_{x}$ and $\alpha_{y}$ are defined by Eqs. (4) and (5), respectively.

Using Eqs. (75), (76) and (77) for the matrices $\boldsymbol{D_{xx}}$, $\boldsymbol{D_{yy}}$ and $\boldsymbol{F_{\pm}}$ of size $\mathcal{M}\times\mathcal{M}$, we obtain the eigenvalue problem

[TABLE]

for the eigenvector

[TABLE]

determining the values of $\psi(x,y)$ at the grid points $(x_{i},y_{j})$ given by Eqs. (73) and (74).

We emphasize that the size of the matrices used in Eq. (82) reduces by a factor of four when we take advantage of the symmetries of Eq. (81) with respect to the transformations $x\rightarrow-x$ and $y\rightarrow-y$. In our calculations we have made use of these symmetries and modified a method Fornberg (1995); Trefethen (2000) originally suggested to improve pseudospectral grids for polar and spherical geometries. In this way we could reduce the size of our matrices while keeping the same accuracy.