11institutetext: School of Computing, Australian National University,
Canberra, ACT 2601, Australia. 11email<EMAIL_ADDRESS>
# Isogeometric Analysis of Bound States of a Quantum Three-Body Problem in 1D

Quanling Deng 11

###### Abstract

In this paper, we initiate the study of isogeometric analysis (IGA) of a
quantum three-body problem that has been well-known to be difficult to solve.
In the IGA setting, we represent the wavefunctions by linear combinations of
B-spline basis functions and solve the problem as a matrix eigenvalue problem.
The eigenvalue gives the eigenstate energy while the eigenvector gives the
coefficients of the B-splines that lead to the eigenstate. The major
difficulty of isogeometric or other finite-element-method-based analyses lies
in the lack of boundary conditions and a large number of degrees of freedom
required for accuracy. For a typical many-body problem with attractive
interaction, there are bound and scattering states where bound states have
negative eigenvalues. We focus on bound states and start with the analysis for
a two-body problem. We demonstrate through various numerical experiments that
IGA provides a promising technique to solve the three-body problem.

###### Keywords:

Isogeometric analysis Three-body problem Bound state

##  1 Introduction

While there are still unsolved questions in the classical three-body problem
[5, 26], the quantum mechanical three-body problem also has unanswered
questions [15, 34, 20]. Accurate and rigorous solutions are highly desirable
both for answering these open questions as well as for studies of three-body
correlations within various many-body systems. The two-body problem is
generally considered as “solved" due to the momentum conservation that leads
to a second-order ordinary differential equation which can be solved
effectively. For three-body problem, the space is six-dimensional in the
center of mass system. The total angular momentum conservation leads to three
coupled second-order nonlinear differential equations in classical mechanics
[27]. In quantum mechanics, we have a more complicated system that admits no
analytic solution in general and existing numerical methods are not
satisfactory in the sense of robustness, efficiency, and reliability.

In quantum chemistry and molecular physics, the Born–Oppenheimer (BO)
approximation has been the most well-known and widely-used mathematical
approximation since the early days of quantum mechanics [3, 28, 8, 32]. The
method is based on the assumption that the nuclei are much heavier than the
electrons which consequently leads to that the wave functions of atomic nuclei
and electrons in a molecule can be treated separately. For instance, BO was
used in [18, 16] to study the Efimov effect in few-body systems. BO was
adopted in [21, 22] recently to establish the universality in a one-
dimensional three-body system. A more efficient numerical method based on BO
was developed using the tensor-product structure [35]. BO has been the
standard method to describe the interaction between electrons and nuclei but
it can fail whenever the assumption fails, for example, in graphene [30]. When
the mass ratios of the interacting bodies are close to one, the assumption
fails and BO is generally inaccurate. Other methods such as the pseudospectral
method based on Fourier analysis [4] and Skorniakov and Ter-Martirosian (STM)
method based on exact integral equations [33] have been developed to obtain
the three-body bound states for arbitrary mass ratios.

In this paper, we develop a general numerical method to solve one-dimensional
quantum two- and three-body problems with arbitrary mass ratios and any
interaction potentials that lead to bound states. With this goal in mind, we
initiate the study of finite element analysis (FEA) based methods to find the
bound states of three-body systems. In particular, we adopt the more advanced
method isogeometric analysis (IGA) for this purpose. IGA, first developed in
[23, 11], has been widely-used as a numerical analysis tool for various
simulations that are governed by partial differential equations (PDEs). IGA
adopts the framework of classic Galerkin FEA and uses B-splines or non-uniform
rational basis splines (NURBS) instead of the Lagrange polynomials as its
basis functions. These basis functions have higher-order continuity
(smoothness) which consequently improves the accuracy of the FEA numerical
approximations. The work [10] applied IGA to study a structural vibration
problem that is modeled as a Laplacian eigenvalue problem. It has been shown
that IGA improved the accuracy of the spectral approximation significantly
compared with FEA [25]. Further advantages of IGA over FEA on spectral
accuracy have been studied in [24, 31]. With the advantages in mind, we adopt
IGA to solve the quantum three-body problem as a second-order differential
eigenvalue problem.

The rest of this paper is organized as follows. Section 2 presents the two-
and three-body problems under consideration. We then unify these two problems
as a single differential eigenvalue problem in one or two dimensions where 1D
refers to the two-body problem and 2D refers to the three-body problem. We
show an example of solutions to a two-body problem for both bound and
scattering states, which serves as a motivation of the proposed method that
solves only the bound states over an approximate finite domain. We then
present the IGA discretization method in Section 3 to solve the unified
problem for the bound states. Section 4 collects and discusses various
numerical tests to demonstrate the performance of the proposed method. We also
perform the numerical study of the impact of domain size on the approximation
accuracy of the bound states. Concluding remarks are presented in Section 5.

##  2 The two- and three-body problems

In this section, we first present the heavy-light two-body and heavy-heavy-
light three-body problems that are modeled as the dimensionless stationary
Schrödinger equations recently studied in [21, 22]. We then generalize the
problems for any mass ratios and unify them as a single differential
eigenvalue problem. A numerical example is followed to show the bound and
scattering states of a two-body problem. The shape of bound states gives a
motivation to pose the differential eigenvalue problem on a finite domain with
a size to be specified depending on the differential operator and accuracy
tolerance.

The heavy-light quantum two-body system with an attractive interaction via a
potential of finite range, after eliminating the center-of-mass motion, is
modeled as a dimensionless stationary Schrödinger equation

$\Big{[}-\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}-v(x)\Big{]}\psi^{(2)}=E^{(2)}\psi^{(2)},$ (1)

where $E^{(2)}$ is the binding energy and $\psi^{(2)}$ is the two-body wave
function. The corresponding three-body system is modeled as

$\Big{[}-\frac{\alpha_{x}}{2}\frac{\partial^{2}}{\partial
x^{2}}-\frac{\alpha_{y}}{2}\frac{\partial^{2}}{\partial
y^{2}}-v(x+y/2)-v(x-y/2)\Big{]}\psi=E\psi,$ (2)

where $E$ is the eigenenergy and $\psi=\psi(x,y)$ is the three-body wave
function describing the relative motions. The coefficients

$\alpha_{x}=\frac{1/2+m_{h}/m_{l}}{1+m_{h}/m_{l}},\qquad\alpha_{y}=\frac{2}{1+m_{h}/m_{l}},$
(3)

where $m_{h}$ denotes the mass of two heavy particles and $m_{l}$ denotes the
mass of the light particle. The potential

$v(\xi)=\beta f(\xi)$ (4)

with $\beta>0$ denoting a magnitude and $f$ denoting the shape of the
interaction potential. We assume that the $f$ is symmetric and describes a
short-range interaction, that is, $|\xi|^{2}f(|\xi|)\to 0$ as
$|\xi|\to\infty$.

###  2.1 The unified problem

The two-body problem (1) is posed on an infinite domain $\Omega=\mathbb{R}$
while the three-body problem (2) is posed on $\Omega=\mathbb{R}^{2}$.
Mathematically, problems (1) and (2) are differential eigenvalue problems
where the differential operator is a Hamiltonian. Moreover, we observe that
(1) is of one variable while (2) is of two variables which can be regarded as
1D and 2D spatial variables, respectively. With this in mind, we unify
problems (1) and (2) to obtain a differential eigenvalue problem

$\displaystyle-\nabla\cdot(\kappa\nabla u)-\gamma u$ $\displaystyle=\lambda
u\quad$ $\displaystyle\forall\ {\boldsymbol{x}}\in\Omega,$ (5)

where $\nabla$ is the gradient operator, $\nabla\cdot$ is the divergence
operator. $\lambda=E^{(2)},\gamma=v(x),\kappa=\frac{1}{2}$ in 1D while
$\lambda=E,\gamma=v(x+y/2)+v(x-y/2)$ and
$\kappa=(\frac{\alpha_{x}}{2},0;0,\frac{\alpha_{y}}{2})$ being a diagonal
matrix in 2D. Herein, $u$ denotes an eigenstate.

From now on, we focus on the unified problem (5). There are three major
difficulties in solving this problem using a Galerkin FEA-based discretization
method.

  * •

(a) The attractive interaction may lead to negative eigenvalues. Consequently,
the discretization of the differential operator
$\mathcal{L}=-\nabla\cdot(\kappa\nabla)-\gamma$ leads to a stiffness matrix
that is not necessarily positive-definite. This in return brings a potential
issue when solving the resulting linear algebra problem.

  * •

(b) The domain $\Omega$ is infinite. This makes it impossible to discretize
the domain with a finite number of elements with each element being of finite
size.

  * •

(c) There are no boundary conditions provided. A Galerkin FEA-based
discretization method requires setting appropriate boundary conditions for the
resulting linear algebra system to be non-singular.

For (a), an eigenvalue shift will resolve the issue. That is, we rewrite
$-\nabla\cdot(\kappa\nabla u)-\gamma u=\lambda u$ by adding a positive scale
to obtain $-\nabla\cdot(\kappa\nabla
u)-(\gamma-\gamma_{0})u=(\lambda+\gamma_{0})u$ where $\gamma_{0}>0$ is a
constant such that $\gamma-\gamma_{0}<0$ for all $\boldsymbol{x}\in\Omega$.
With a slight abuse of notation, the problem (5) can be rewritten as

$\displaystyle-\nabla\cdot(\kappa\nabla u)+\gamma u$ $\displaystyle=\lambda
u\quad$ $\displaystyle\forall\ {\boldsymbol{x}}\in\Omega.$ (6)

To overcome the difficulties (b) and (c), we first present an example of a
solution to the two-body problem in the next subsection.

###  2.2 A solution example of the two-body problem

For attractive interaction, the eigenenergies for certain eigenstates can be
negative. For a potential vanishing at $\pm\infty$, a negative eigenvalue
implies a bound state while a positive eigenvalue implies a scattering state
[19].

Figure 1: Bound and scattering states of two-body problems.

Figure 1 shows an example of state solutions to the two-body problem (1) with
$\kappa=1/2$, potential $f(\xi)=e^{-\xi^{2}}$ and $\beta=1$ for the left plot
while $\beta=2$ for the right plot. The bound states eigenenergies are marked
in the figure. Herein, we apply IGA with 5000 elements and $C^{6}$ septic
B-spline basis functions. We present the details of the IGA method in the next
section. We observe that there is one bound state for $\beta=1$ and two bound
states for $\beta=2$. All other states are scattering states. When
$x\to\pm\infty$, the wavefunctions go to zeros exponentially fast for bound
states while they do not go to zeros for scattering states. A theoretical
explanation can be found in [1].

In this paper, our goal is to find the eigenenergies and eigenstates for bound
states. In the case of $\beta=1$, the bound state solution decays to zero
approximately at $x=\pm 10$ with an error of $5.9\times 10^{-5}$. We observe
similar behaviour for the two bound states of the case $\beta=2$. This
decaying behaviour provides an insight to overcome the difficulties (b) and
(c) listed in section 2.1. The idea is that for a given tolerance
$\epsilon>0$, we propose to solve the problem (6) on a finite domain
$\Omega_{\epsilon}=[-\boldsymbol{x}_{\epsilon},\boldsymbol{x}_{\epsilon}]^{d},d=1,2$
with homogeneous boundary condition

$u=0,\qquad\forall\ {\boldsymbol{x}}\in\partial\Omega_{\epsilon}.$ (7)

###### Remark 1

For smaller error tolerance, one expects to apply a larger finite domain. A
detailed study is presented in section 4.2. We also point out that a simple
transformation such as $\hat{x}=\tanh(x)$ that transfers the infinite domain
$\Omega$ to a finite domain $\hat{\Omega}=(-1,1)^{d},d=1,2$ can not avoid the
difficulties of (b) and (c) listed in section 2.1.

##  3 Isogeometric analysis

In this section, we present the IGA method for the unified problem (6) on
$\Omega_{\epsilon}$ supplied with the boundary condition (7). We also give an
a priori error estimate for the bound states and their eigenenergies.

###  3.1 Continuous level

Let
$\Omega_{\epsilon}=[-\boldsymbol{x}_{\epsilon},\boldsymbol{x}_{\epsilon}]^{d}\subset\mathbb{R}^{d},d=1,2$
be a bounded domain with Lipschitz boundary $\partial\Omega_{\epsilon}$. We
adopt the standard notation for the Hilbert and Sobolev spaces. For a
measurable subset $S\subseteq\Omega_{\epsilon}$, we denote by
$(\cdot,\cdot)_{S}$ and $\|\cdot\|_{S}$ the $L^{2}$-inner product and its
norm, respectively. We omit the subscripts when clear. For an integer $m\geq
1$, we denote the $H^{m}$-norm and $H^{m}$-seminorm as $\|\cdot\|_{H^{m}(S)}$
and $|\cdot|_{H^{m}(S)}$, respectively. In particular, we denote by
$H^{1}_{0}(\Omega_{\epsilon})$ the Sobolev space with functions in
$H^{1}(\Omega_{\epsilon})$ that are vanishing on the boundaries.

The variational formulation of (6) is to find eigenvalue
$\lambda\in\mathbb{R}^{+}$ and eigenfuction $u\in
H^{1}_{0}(\Omega_{\epsilon})$ with $\|u\|_{\Omega_{\epsilon}}=1$ such that

$a(w,u)=\lambda b(w,u),\quad\forall\ w\in H^{1}_{0}(\Omega_{\epsilon}),$ (8)

where the bilinear forms are defined as for $v,w\in
H^{1}_{0}(\Omega_{\epsilon})$

$a(v,w):=(\kappa\nabla v,\nabla w)_{\Omega_{\epsilon}}+(\gamma
v,w)_{\Omega_{\epsilon}},\qquad b(v,w):=(v,w)_{\Omega_{\epsilon}}.$ (9)

The eigenvalue problem (8) with
$\gamma=\gamma({\boldsymbol{x}})>0,\forall{\boldsymbol{x}}\in\Omega_{\epsilon},$
has a countable set of positive eigenvalues (see, for example, [6, Sec. 9.8])

$0<\lambda_{1}<\lambda_{2}\leq\lambda_{3}\leq\cdots$

with an associated set of orthonormal eigenfunctions
$\\{u_{j}\\}_{j=1}^{\infty}$. Thus, there holds $(u_{j},u_{k})=\delta_{jk},$
where $\delta_{jk}=1$ is the Kronecker delta. As a consequence, the
eigenfunctions are also orthogonal in the energy inner-product as there holds
$a(u_{j},u_{k})=\lambda_{j}b(u_{j},u_{k})=\lambda_{j}\delta_{jk}.$

###  3.2 IGA discretized level

At the discretized level, we first discretize the domain $\Omega_{\epsilon}$
with a uniform tensor-product mesh. We denote a general element as $\tau$ and
its collection as $\mathcal{T}_{h}$ such that
$\overline{\Omega}_{\epsilon}=\cup_{\tau\in\mathcal{T}_{h}}\tau$. Let
$h=\max_{\tau\in\mathcal{T}_{h}}\text{diameter}(\tau)$. In the IGA setting,
for simplicity, we use the B-splines. The B-spline basis functions in 1D are
given as the Cox-de Boor recursion formula; we refer to [12, 29] for details.
Let $X=\\{x_{0},x_{1},\cdots,x_{m}\\}$ be a knot vector with a nondecreasing
sequential knots $x_{j}$. The $j$-th B-spline basis function of degree $p$,
denoted as $\phi^{j}_{p}(x)$, is defined recursively as

$\displaystyle\phi^{j}_{0}(x)$ $\displaystyle=\begin{cases}1,\quad\text{if}\
x_{j}\leq x<x_{j+1},\\\ 0,\quad\text{otherwise},\\\ \end{cases}$ (10)
$\displaystyle\phi^{j}_{p}(x)$
$\displaystyle=\frac{x-x_{j}}{x_{j+p}-x_{j}}\phi^{j}_{p-1}(x)+\frac{x_{j+p+1}-x}{x_{j+p+1}-x_{j+1}}\phi^{j+1}_{p-1}(x).$

A tensor-product of these 1D B-splines produces the B-spline basis functions
in multiple dimensions. We define the multi-dimensional approximation space as
$V^{h}_{p}\subset H^{1}_{0}(\Omega_{\epsilon})$ with:

$V^{h}_{p}=\text{span}\\{\phi^{j}_{p}\\}_{j=1}^{N_{h}}=\begin{cases}\text{span}\\{\phi^{j_{x}}_{p_{x}}(x)\\}_{j_{x}=1}^{N_{x}},&\text{in
1D},\\\
\text{span}\\{\phi^{j_{x}}_{p_{x}}(x)\phi^{j_{y}}_{p_{y}}(y)\\}_{j_{x},j_{y}=1}^{N_{x},N_{y}},&\text{in
2D},\\\ \end{cases}$

where $p_{x},p_{y}$ specify the approximation order in each dimension.
$N_{x},N_{y}$ is the total number of basis functions in each dimension and
$N_{h}$ is the total number of degrees of freedom. The isogeometric analysis
of (6) in variational formulation seeks $\lambda^{h}\in\mathbb{R}$ and
$u^{h}\in V^{h}_{p}$ with $\|u^{h}\|_{\Omega_{\epsilon}}=1$ such that

$a(w^{h},u^{h})=\lambda^{h}b(w^{h},u^{h}),\quad\forall\ w^{h}\in V^{h}_{p}.$
(11)

###  3.3 Algebraic level

At the algebraic level, we approximate the eigenfunctions as a linear
combination of the B-spline basis functions, i.e.,

$u^{h}=\sum_{j=1}^{N_{h}}\nu_{j}\phi^{j}_{p},$

where $\nu_{j},j=1,\cdots,N_{h}$ are the coefficients. We then substitute all
the B-spline basis functions for $w^{h}$ in (11). This leads to the
generalized matrix eigenvalue problem

$\mathbf{K}\mathbf{U}=\lambda^{h}\mathbf{M}\mathbf{U},$ (12)

where
$\mathbf{K}_{kl}=a(\phi_{p}^{k},\phi_{p}^{l}),\mathbf{M}_{kl}=b(\phi_{p}^{k},\phi_{p}^{l}),$
and $\mathbf{U}$ is the corresponding representation of the eigenvector as the
coefficients of the B-spline basis functions. The homogeneous Dirichlet
boundary condition (7) can be set by removing the rows and columns
corresponding to the degrees of freedom associated with the nodes at the
boundary. This matrix eigenvalue problem is to be solved in a computing
program.

###  3.4 A priori error estimates

IGA is a Galerkin finite element discretization method. On a rectangular
domain with tensor-product grids, the only difference of IGA from the
classical FEA is the basis functions. FEA adopts $C^{0}$ polynomials as basis
function while IGA adopts $C^{k},k\geq 1$ polynomials. We observe that $C^{0}$
is a larger space, i.e., $C^{k}(\Omega_{\epsilon})\subset
C^{0}(\Omega_{\epsilon})$. In general, for an a priori error estimate that is
established in the Galerkin FEA framework, the estimate also holds for IGA.
Thus, we expect optimal convergence rates for the eigenvalues and
eigenfunctions as in FEA [9, 2, 17]. We present the following estimate without
a theoretical proof. Instead, we show numerical validation in Section 4.

Given the mesh configuration and IGA setting described above, let
$(\lambda_{j},u_{j})\in\mathbb{R}^{+}\times H^{1}_{0}(\Omega)$ solve (8) for
bound states and let $(\lambda_{j}^{h},u_{j}^{h})\in\mathbb{R}^{+}\times
V^{h}_{p}$ solve (11) for bound states with the normalizations
$\|u_{j}\|_{\Omega_{\epsilon}}=1$ and $\|u_{j}^{h}\|_{\Omega_{\epsilon}}=1$.
Assuming elliptic regularity on the operator
$\mathcal{L}=-\nabla\cdot(\kappa\nabla)+\gamma$ and high-order smoothness of
the eigenfunctions $u_{j}$ on $\Omega_{\epsilon}$, there holds:

$\big{|}\lambda_{j}^{h}-\lambda_{j}\big{|}\leq
Ch^{2p},\qquad|u_{j}-u_{j}^{h}|_{H^{1}(\Omega)}\leq Ch^{p},$ (13)

where $C$ is a positive constant independent of the mesh-size $h$. We remark
that these estimates only hold for bound states and do not necessarily hold
for scattering states.

##  4 Numerical experiments

In this section, we present various numerical examples to demonstrate the
performance of IGA. We first show the IGA approximation optimal convergence
accuracy with a domain $\Omega_{\epsilon}$ of large size. Then we study the
impact of the domain size on accuracy and give an approximate formula that
determines the size of the domain given a certain accuracy tolerance.

###  4.1 IGA discretization accuracy

We focus on the two- and three-body problems with a potential with polynomial
decay

$f(\xi)=\frac{1}{(1+\xi^{2})^{3}}$ (14)

of the cube of a Lorentzian and one with exponential decay

$f(\xi)=e^{-\xi^{2}}$ (15)

of a Gaussian. For these potentials, finding the exact analytical solutions is
impossible. For the purpose of characterizing the errors, we use, as a
reference solution to the exact one, the solution of IGA with a septic $C^{6}$
B-spline basis functions and a fine mesh. We focus on the eigenvalue error
that is defined as

$e_{j}=|\lambda^{h}_{j}-\hat{\lambda}_{j}|,$ (16)

where $\lambda^{h}_{j}$ is an IGA eigenvalue and $\hat{\lambda}_{j}$ is a
reference eigenvalue that is of high accuracy approximating the exact one
$\lambda_{j}$.

Figure 2 shows the eigenvalue error convergence rates for IGA of the two-body
problem with $\kappa=\frac{1}{2}$ in 1D. We consider $C^{0}$ linear, $C^{1}$
quadratic, and $C^{2}$ cubic IGA elements. We study the problem with both
potentials (14) and (15) and a fixed magnitude $\beta=1$. For both potentials
(14) and (15), there is one bound state. The state reference eigenvalue is
$\hat{\lambda}_{1}=-0.31658012845$ for (14) and
$\hat{\lambda}_{1}=-0.47738997738$ for (15), respectively. We solve the
problem for the bound state using $C^{6}$ septic IGA with 5000 uniform
elements over the domain $\Omega_{\epsilon}=[-20,20]$. We observe optimal
error convergence rates in all the scenarios. This confirms the theoretical
prediction (13) in section 3.4.

Figure 2: Eigenvalue error convergence rates of IGA for the two-body problem
with $\beta=1$ on a domain $\Omega_{\epsilon}=[-20,20].$ The potential has
polynomial decay (14) for the left plot while exponential decay (15) for the
right plot.

Now we consider a case where there are two bound states in the two-body
problem. Let $\beta=5$ and we apply a potential of polynomial decay (14).
Figure 3 shows the two bound states solutions and their eigenenergies are
$-2.9149185630$ and $-0.25417134380$. Herein, the numerical eigenstates are
computed using $C^{6}$ septic IGA with 5000 uniform elements over the domain
$\Omega_{\epsilon}=[-20,20]$. The plot shows the eigenstate over $[-10,10]$
for better focus while the problem is solved over the larger domain
$\Omega_{\epsilon}=[-20,20]$ for high accuracy.

Figure 3: The bound states of the two-body problem with $\beta=5$ and a
potential of polynomial decay (14).  Figure 4:  Eigenvalue error convergence
rates of IGA for the two-body problem with $\beta=5$ and a polynomially
decaying potential (14) on a domain $\Omega_{\epsilon}=[-20,20].$ The left
plot shows the eigenenergy errors of first bound state while the right plot
shows these of the second bound state.

Figure 4 shows the eigenenergy error convergence rates of the problem
described above (also shown in Figure 3). Again, we observe optimal error
convergence rates that verify the theoretical prediction. Moreover, for IGA
with higher-order elements, the eigenvalue errors reach small errors faster
with coarser meshes. This validates that the reference solutions obtained by
using $C^{6}$ septic IGA with 5000 elements are of high accuracy and can be
used as highly accurate approximations to the exact solutions.

###  4.2 A study on domain size

For bound states of the two-body problems as discussed in section 2.2, the
state values approach zero exponentially fast. The IGA discretization requires
a finite domain with homogeneous boundary condition (7). The accuracy depends
on the domain size $x_{\epsilon}$ for
$\Omega_{\epsilon}=[-x_{\epsilon},x_{\epsilon}]$.

Figure 5: Eigenvalue errors versus domain size $x_{\epsilon}$ for the two-body
problem with polynomially decaying potential (14) (left plot) and
exponentially decaying potential (15) (right plot).

To study the impact of the domain size $x_{\epsilon}$ on the accuracy, we
apply the high-accuracy IGA method with $C^{6}$ septic B-spline elements. We
apply uniform mesh grids with a fixed grid size $h=0.01$. This setting of
using a high-order element with fine grid size is to guarantee that the errors
are dominated by the choice of the domain size. Figure 5 shows how the
eigenvalue errors decrease when the domain size $x_{\epsilon}$ increases. We
set $\kappa=1/2$ in (6) and the potential magnitude $\beta=1$. The left plot
of Figure 5 shows that the eigenvalue error decays exponentially when
$x_{\epsilon}$ increases for the potential (14) while the right plot of Figure
5 shows that of the potential (15). The fitted functions that establish the
relation between the error and the domain size are
$e_{1}=10^{-0.72x_{\epsilon}+0.37}$ and $e_{1}=10^{-0.87x_{\epsilon}+0.61}$
for (14) and (15), respectively. Figure 6 shows the case where there are two
bound states. Therein, potential (15) is used with $\beta=5$. We observe a
similar behaviour. The errors of the ground state reach an order of $10^{-12}$
when $x_{\epsilon}\geq 6.5$. This is due to that the IGA discretization error
dominates the overall error (from discretization and approximation of the
domain). The fitted functions give guidance for choosing the domain
$\Omega_{\epsilon}$ appropriately. For example, for the two-body problem with
a potential (15) and $\beta=5$, the fitted function for the second bound state
(a larger domain is required to compute this mode) is
$e_{2}=10^{-0.68x_{\epsilon}+0.69}$. Thus, to achieve an accuracy of error
$10^{-15}$, we set $10^{-15}=10^{-0.68x_{\epsilon}+0.69}$ and solve for
$x_{\epsilon}$ to get the domain
$\Omega_{\epsilon}=[-x_{\epsilon},x_{\epsilon}]=[-23.1,23.1]$. This means that
we require to solve the problem with a minimal domain size $x_{\epsilon}=23.1$
to get an accuracy of order $10^{-15}$.

Figure 6: Eigenvalue errors versus domain size $x_{\epsilon}$ for the two
bound states of the two-body problem with exponentially decaying potential
(15).

###  4.3 Three-body problem

Now, we consider the three-body problem with a heavy-light body ratio
$m_{h}/m_{l}=20$ that is studied in [21, 22, 35]. With such a mass ratio,
$\kappa$ in the unified problem (6) is a matrix with entries
$\kappa=(41/84,0;0,2/21)$. Our goal is to approximate the eigenenergies and
eigenstates obtained by using the classical BO approximation in these papers.
This preliminary numerical study demonstrates that the proposed IGA method is
a promising alternative to the classical BO approximation method that a strong
assumption is posed on the mass ratio.

To solve the three-body problem, we apply the highly accurate IGA method with
$C^{6}$ septic B-spline elements. We set the domain as
$\Omega_{\epsilon}=[-20,20]$ and apply a non-uniform grid with $80\times 80$
elements. Table 1 shows the eigenvalues of the bound states of the three-body
problem with the exponentially decaying potential (15). The potential
magnitude is $\beta=0.344595351$. The IGA eigenvalues are close to the ones
(scaled) shown in Table 1 of [35]. Figure 7 shows the first four bound state
eigenfunctions. The eigenstate solution shapes match well with the ones
obtained using the BO approximation in Figure 4 of [21]. Moreover, we observe
a similar universality behaviour as in [21] and we will present a detailed
study in future work. In conclusion, the IGA method with a small mesh grid has
the ability to approximate well both eigenenergies and eigenfunctions of the
bound states of the three-body problem.

method | $\beta$ |  $j$ (Bosons) | $\lambda_{j}^{h}$ |  $j$ (Fermions) | $\lambda_{j}^{h}$
---|---|---|---|---|---
|  | 0 | -0.2476034576 | 1 | -0.1825896533
IGA | 0.344595351 | 2 | -0.1412793292 | 3 | -0.1182591543
|  | 4 | -0.1060931444 | 5 | -0.1005294105
|  | 0 | -0.247603458 | 1 | -0.182589653
BO in [35] | 0.34459535 | 2 | -0.141279329 | 3 | -0.118259157
|  | 4 | -0.106093864 | 5 | -0.102845702
Table 1: Eigenenergies of the three-body problem with a mass ratio 20 and an
exponentially decaying potential (15) when using IGA and BO approximation in
[35]. Figure 7:  Eigenstate wavefunctions of the three-body problem with a
mass ratio 20, potential magnitude $\beta=0.344595351$, and the exponentially
decaying potential (15).

##  5 Concluding remarks

In this paper, we initiated IGA of the quantum two- and three-body problems.
IGA is developed based on the classical Galerkin FEA that has the advantages
of a mature theoretical understanding of the error estimates, stabilities, and
robustness. IGA is successfully applied to solve the bound states of the two-
and three-body problems in 1D with arbitrary mass ratios and potential shapes.

As for future work, the first possible direction is a generalization to the
two- and three-body problems in multiple dimensions. Tensor-product structures
may be applied to reduce the computational costs [35]. Another direction of
future work is that one may use the recently developed softFEM [14] and
dispersion-minimized blending quadratures [13, 7] to solve the three-body
problem with higher accuracy and efficiency.

## References

  * [1] Agmon, S.: Lectures on exponential decay of solutions of second-order elliptic equations. Princeton University Press (2014)
  * [2] Babuška, I., Osborn, J.: Eigenvalue problems. In: Handbook of Numerical Analysis, Vol. II, pp. 641–787. Handb. Numer. Anal., II, North-Holland, Amsterdam (1991)
  * [3] Baer, M.: Beyond Born–Oppenheimer: electronic nonadiabatic coupling terms and conical intersections. John Wiley & Sons (2006)
  * [4] Boyd, J.P.: Chebyshev and Fourier spectral methods. Courier Corporation (2001)
  * [5] Breen, P.G., Foley, C.N., Boekholt, T., Zwart, S.P.: Newton versus the machine: solving the chaotic three-body problem using deep neural networks. Monthly Notices of the Royal Astronomical Society 494(2), 2465–2470 (2020)
  * [6] Brezis, H.: Functional analysis, Sobolev spaces and partial differential equations. Universitext, Springer, New York (2011)
  * [7] Calo, V., Deng, Q., Puzyrev, V.: Dispersion optimized quadratures for isogeometric analysis. J. Comput. Appl. Math. 355, 283–300 (2019)
  * [8] Cederbaum, L.S.: Born–Oppenheimer approximation and beyond for time-dependent electronic processes. The Journal of chemical physics 128(12), 124101 (2008)
  * [9] Ciarlet, P.G.: Finite Element Method for Elliptic Problems. Society for Industrial and Applied Mathematics, Philadelphia, PA, USA (2002)
  * [10] Cottrell, J.A., Reali, A., Bazilevs, Y., Hughes, T.J.R.: Isogeometric analysis of structural vibrations. Comput. Methods Appl. Mech. Engrg. 195(41-43), 5257–5296 (2006)
  * [11] Cottrell, J.A., Hughes, T.J.R., Bazilevs, Y.: Isogeometric analysis: toward integration of CAD and FEA. John Wiley & Sons (2009)
  * [12] De Boor, C.: A practical guide to splines, vol. 27. Springer-Verlag New York (1978)
  * [13] Deng, Q., Calo, V.: Dispersion-minimized mass for isogeometric analysis. Comput. Methods Appl. Mech. Engrg. 341, 71–92 (2018)
  * [14] Deng, Q., Ern, A.: SoftFEM: revisiting the spectral finite element approximation of second-order elliptic operators. Computers & Mathematics with Applications 101, 119–133 (2021)
  * [15] Drut, J.E., McKenney, J.R., Daza, W.S., Lin, C.L., Ordóñez, C.R.: Quantum anomaly and thermodynamics of one-dimensional fermions with three-body interactions. Physical Review Letters 120(24), 243002 (2018)
  * [16] Efremov, M.A., Plimak, L., Berg, B., Ivanov, M.Y., Schleich, W.P.: Efimov states in atom-molecule collisions. Physical Review A 80(2), 022714 (2009)
  * [17] Ern, A., Guermond, J.L.: Finite Elements II: Galerkin approximation, elliptic and mixed PDEs. Springer-Verlag, New York (2020), in press
  * [18] Fonseca, A.C., Redish, E.F., Shanley, P.: Efimov effect in an analytically solvable model. Nuclear Physics A 320(2), 273–288 (1979)
  * [19] Griffiths, D.J., Schroeter, D.F.: Introduction to quantum mechanics. Cambridge university press (2018)
  * [20] Guo, P., Gasparian, V.: Numerical approach for finite volume three-body interaction. Physical Review D 97(1), 014504 (2018)
  * [21] Happ, L., Zimmermann, M., Betelu, S.I., Schleich, W.P., Efremov, M.A.: Universality in a one-dimensional three-body system. Physical Review A 100(1), 012709 (2019)
  * [22] Happ, L., Zimmermann, M., Efremov, M.A.: Universality of excited three-body bound states in one dimension. Journal of Physics B: Atomic, Molecular and Optical Physics 55(1), 015301 (2022)
  * [23] Hughes, T.J.R., Cottrell, J.A., Bazilevs, Y.: Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer methods in applied mechanics and engineering 194(39), 4135–4195 (2005)
  * [24] Hughes, T.J.R., Evans, J.A., Reali, A.: Finite element and NURBS approximations of eigenvalue, boundary-value, and initial-value problems. Computer Methods in Applied Mechanics and Engineering 272, 290–320 (2014)
  * [25] Hughes, T.J.R., Reali, A., Sangalli, G.: Duality and unified analysis of discrete approximations in structural dynamics and wave propagation: comparison of $p$-method finite elements with $k$-method NURBS. Comput. Methods Appl. Mech. Engrg. 197(49-50), 4104–4124 (2008)
  * [26] Letellier, C.: Chaos in nature, vol. 94. World Scientific (2019)
  * [27] Nielsen, E., Fedorov, D.V., Jensen, A.S., Garrido, E.: The three-body problem with short-range interactions. Physics Reports 347(5), 373–459 (2001)
  * [28] Panati, G., Spohn, H., Teufel, S.: The time-dependent Born–Oppenheimer approximation. ESAIM: Mathematical Modelling and Numerical Analysis 41(2), 297–314 (2007)
  * [29] Piegl, L., Tiller, W.: The NURBS book. Springer Science & Business Media (1997)
  * [30] Pisana, S., Lazzeri, M., Casiraghi, C., Novoselov, K.S., Geim, A.K., Ferrari, A.C., Mauri, F.: Breakdown of the adiabatic Born–Oppenheimer approximation in graphene. Nature materials 6(3), 198–201 (2007)
  * [31] Puzyrev, V., Deng, Q., Calo, V.M.: Dispersion-optimized quadrature rules for isogeometric analysis: modified inner products, their dispersion properties, and optimally blended schemes. Computer Methods in Applied Mechanics and Engineering 320, 421–443 (2017)
  * [32] Scherrer, A., Agostini, F., Sebastiani, D., Gross, E., Vuilleumier, R.: On the mass of atoms in molecules: Beyond the Born–Oppenheimer approximation. Physical Review X 7(3), 031035 (2017)
  * [33] Skorniakov, G., Ter-Martirosian, K.: Three body problem for short range forces. i. scattering of low energy neutrons by deuterons. Soviet Phys. JETP 4 (1957)
  * [34] Sukhareva, O., Grigorenko, L., Kostyleva, D., Zhukov, M.: Validity of quasi-classical approaches to true three-body decays. In: International Conference on Few-Body Problems in Physics. pp. 283–286. Springer (2018)
  * [35] Thies, J., Hof, M.T., Zimmermann, M., Efremov, M.: Exploiting tensor structure for computing bound states of the quantum mechanical three-body problem. arXiv preprint arXiv:2111.02534 (2021)