# Nuclear reactions in artificial traps

Peng Guo<EMAIL_ADDRESS>College of Physics, Sichuan University, Chengdu,
Sichuan 610065, China  Department of Physics and Engineering, California State
University, Bakersfield, CA 93311, USA  Kavli Institute for Theoretical
Physics, University of California, Santa Barbara, CA 93106, USA  Bingwei Long
<EMAIL_ADDRESS>College of Physics, Sichuan University, Chengdu, Sichuan
610065, China

###### Abstract

Coupled-channel two-particle systems bound by a harmonic trap are discussed in
the present paper. We derive the formula that relates the energy levels of
such trapped systems to phase shifts and inelasticity of coupled-channel
reactions. The formula makes it possible to extract amplitudes of inelastic
nuclear reactions from ab initio calculations of discrete levels of many-
nucleon systems in a harmonic trap.

##  I Introduction

Studies of nuclear structure have entered a stage where bound states of many
nucleons can be solved for with microscopic nuclear forces and with high
precision Navratil et al. (2007); Epelbaum et al. (2009); Ekström et al.
(2013); Barrett et al. (2013); Lynn et al. (2014); Ekström et al. (2015);
Piarulli et al. (2018); Hammer et al. (2020). For some of those many-body
solvers, the computation typically involves expansion of wave functions on a
chosen basis, e.g., the harmonic-oscillator (HO) basis. This poses a challenge
for immediate extension of basis-expansion methods to continuum-state problems
like resonances and reactions. The need for synergy of bound-state and
continuum computational methods has been recognized by the nuclear-theory
community Navrátil et al. (2016); Yang (2016); Johnson et al. (2019); Mazur et
al. (2020); Ma et al. (2020).

One of these ab initio methods for elastic scattering of two clusters,
normally stable nuclei, utilizes the energy spectrum of all the constituent
nucleons, trapped by an artificial field, to generate directly the phase
shifts at certain discrete energies decided by the spectrum of the trapped
Rotureau et al. (2010, 2012); Luu et al. (2010); Zhang et al. (2020). This
method is mostly enabled by a set of so-called quantization conditions (QCs)
that map energy levels to the phase shifts at momenta related to those energy
eigenvalues Busch et al. (1998); Zhang (2020). The harmonic trap is an
attractive choice for the artificial field because the HO wave functions are
analytically known, so such mapping formulas can be presented in closed forms
to a great extent. Existing vast resource of software for the HO basis is also
an important incentive Navratil et al. (2000). After the pioneering work in
Ref. Busch et al. (1998) to connect the spectrum of the trapped two-body
systems to their elastic scattering amplitude, referred to as the BERW formula
in the present paper, more theoretical works to extend the formula and
applications to nucleon-nucleon and nucleon-cluster scattering followed Stetcu
et al. (2007, 2010); Rotureau et al. (2010, 2012); Luu et al. (2010); Yang
(2016); Zhang (2020); Zhang et al. (2020). Here we turn our attention to
coupled-channel problems where at least one of the two participating particles
changes its species after the reaction. This is relevant to nuclear reactions
where the nuclei fuse and then break up to different nucleus than before.

Relating discrete energy spectrum of trapped composite particles to their
scattering amplitudes has been a topic of great interests in other fields of
strong-interaction physics. In these applications, the form of traps may
differ, such as periodic cubic boxes in lattice quantum chromodynamics (QCD)
Lüscher (1991) and spherical hard wall in some of the lattice implementations
of chiral effective field theory Bovermann et al. (2019); Elhatisari et al.
(2016); Rokash et al. (2015). Because the underlying theory to compute the
energy spectrum is model-independent and the computation often consumes much
computing resources, it is highly desired that the said connection is
constructed without uncontrolled modeling of interactions between the
composite particles. This can be achieved when two scales are well separated;
one is the wideness of the trap $b$ and the other the interaction range
between the composite particles $R$, $b\gg R$. The separation makes it
possible to locate a domain inside the trap where the composite particles do
not interact appreciably while the confining force is sufficiently weak so
that the trapped particles’ wave function resembles that of the scattering
state in infinite volume. One then matches the asymptotic form of the
scattering wave function, which has undetermined scattering parameters, say,
phase shifts, to the trapped wave function. Doing so pins down the phase shift
at the particular energy associated with the trapped wave function. This
yields the QC, a compact relation between the discrete energy eigenvalues and
the scattering parameters, which often takes the following schematic form:

$\det\left[\cot\delta(E)-g(E;b)\right]_{E=\mathcal{E}_{k}}=0\,,$ (1)

where $\mathcal{E}_{k}$ is the $k$-th the energy eigenvalue and $\delta(E)$
refers to the diagonal matrix with phase shifts as the elements. The analytic
matrix function $g(E;b)$ depends only on the configuration of the trap, such
as the functional form of the confining force and the wideness $b$, but not on
detail of the inter-particle interactions.

Matching wave functions was the fundamental methodology used by Lüscher to
derive his seminal QC formula for periodic cubic boxes Lüscher (1991). Once we
are convinced that the QC is irrespective of the interactions, manipulation of
wave functions, however, does not have to be the ubiquitous technique of
deriving the QC. Actively researched in the field of lattice QCD, other
methods have been developed to tackle the problem (see, e.g., Refs.
Rummukainen and Gottlieb (1995); Christ et al. (2005); Bernard et al. (2008);
He et al. (2005); Lage et al. (2009); Döring et al. (2011); Guo et al. (2013);
Guo (2013); Kreuzer and Hammer (2009); Polejaeva and Rusetsky (2012); Hansen
and Sharpe (2014); Mai and Döring (2017, 2019); Döring et al. (2018); Guo
(2017); Guo and Gasparian (2017, 2018); Guo and Morris (2019); Guo et al.
(2018); Guo (2020a); Guo and Döring (2020); Guo (2020b); Guo and Long (2020a);
Guo (2020c); Guo and Long (2020b); Guo (2020d); Guo and Gasparian (2021)). We
will follow the variational approach spelled out in Refs. Guo (2020c, d); Guo
and Gasparian (2021). In particular, contact interactions are used in
describing dynamics in both harmonic trap and infinite volume, and the QC is
derived by eliminating the contact coupling.

The paper is organized as follows. To establish our notation, derivation of
the BERW formula is discussed in Sec. II. The extension to coupled-channel
reactions is presented in Sec. III, followed by discussions and summary in
Sec. IV.

##  II BERW formula

Using the approach in Refs. Guo (2020c, d); Guo and Gasparian (2021), we re-
derive the BERW formula for elastic scattering of two spinless particles in
this section. We consider cases where the kinematics is such that the reaction
does not takes places near threshold. So the long-range Coulomb repulsion, if
any, can be neglected and one is left with only the strong interaction, which
is finite-ranged.

###  II.1 Scattering dynamics

In infinite volume, the Lippmann-Schwinger (LS) equation for two-particle
scattering has the following operator form,

$\hat{T}_{\infty}(E)=-\hat{V}+\hat{V}\hat{G}_{\infty}(E)\hat{T}_{\infty}(E)\,,$
(2)

where $\hat{T}_{\infty}$ is the $T$-matrix, $E$ the energy, and $\hat{V}$ the
short-range interaction between the clusters. The infinite-volume Green’s
function is given in the operator form by

$\hat{G}_{\infty}(E)=\frac{1}{E-\hat{H}_{0}}\,,$ (3)

with $\hat{H}_{0}$ defined in terms of the reduced mass $\mu$ and the
Laplacian:

$\hat{H}_{0}=-\frac{\nabla^{2}}{2\mu}\,.$ (4)

The solution to the LS equation can be formally written as

$\hat{T}_{\infty}(E)=-\left[\hat{V}^{-1}-\hat{G}_{\infty}(E)\right]^{-1}\,.$
(5)

The on-shell elements of partial-wave projected $T$-matrix are usually
parameterized in the center-of-mass (CM) frame by the phase shift $\delta_{L}$
as

$T^{(\infty)}_{L}(q,q)=\frac{(4\pi)^{2}}{2\mu
q}\frac{1}{\cot\delta_{L}(q)-i}\,,$ (6)

where $q=\sqrt{2\mu E}$ is relative momentum of particles.

###  II.2 Two particles trapped by a HO potential

One may still define the $T$-matrix for the trapped two-particle system for
any values of $E$, even though there are no scattering states:

$\hat{T}_{\omega}(E)=-\hat{V}+\hat{V}\hat{G}_{\omega}(E)\hat{T}_{\omega}(E)\,,$
(7)

where the Green’s function of the HO potential with angular frequency $\omega$
is given by

$\hat{G}_{\omega}(E)=\frac{1}{E-\hat{H}_{\omega}}=\sum_{n}\frac{|n\rangle\langle
n|}{E-\omega(n+\frac{3}{2})}\,.$ (8)

Here $\hat{H}_{\omega}$ is the HO Hamiltonian,

$\hat{H}_{\omega}=-\frac{\nabla^{2}}{2\mu}+\frac{1}{2}\mu\omega^{2}\mathbf{r}^{2}\,,$
(9)

and $|n\rangle$ stands for eigenstates of $\hat{H}_{\omega}$ operator:

$\hat{H}_{\omega}|n\rangle=\omega(n+\frac{3}{2})|n\rangle,\ \ \ \
n=0,1,2\cdots.$ (10)

The wideness of the harmonic trap can be characterized by the so-called
oscillator length:

$b\equiv\left(\frac{1}{2}\mu\omega\right)^{-\frac{1}{2}}\,.$ (11)

The bound state sustained by the harmonic trap can be determined by poles of
$T_{\omega}(E)$, and

$\hat{T}_{\omega}(E)=-\left[\hat{V}^{-1}-\hat{G}_{\omega}(E)\right]^{-1}\,.$
(12)

Therefore, any eigenvalue of the full Hamiltonian $(\hat{H}_{\omega}+\hat{V})$
$E=\mathcal{E}_{k}$ satisfies

$\det\left[\hat{V}^{-1}-\hat{G}_{\omega}(E)\right]_{E=\mathcal{E}_{k}}=0\,.$
(13)

###  II.3 Contact potentials and the BERW formula

Combining Eqs. (5) and (13) by eliminating $\hat{V}$, one finds symbolically

$\det\left[\hat{T}_{\infty}^{-1}(E)+\hat{G}_{\omega}(E)-\hat{G}_{\infty}(E)\right]_{E=\mathcal{E}_{k}}=0.$
(14)

Although the above equation is formally correct, it is difficult to use in
practical calculations because inverting $T_{\infty}(E)$ directly is highly
non-trivial.

As long as $R/b$ is small, as we argued, a model-independent QC must exist. So
one can use whatever form of two-body interaction to deduce the QC. We choose
a sum of contact interactions, with one term responsible for each partial
wave, to describe the cluster-cluster interaction:

$V(\mathbf{q}^{\prime},\mathbf{q}\,)=\sum_{L}\frac{2L+1}{4\pi}(q^{\prime}q)^{L}V_{L}P_{L}(\hat{\mathbf{q}}^{\prime}\cdot\hat{\mathbf{q}})\,,$
(15)

where $P_{L}(x)$ is the Legendre polynomial and $\mathbf{q}^{\prime}$
($\mathbf{q}$) is the outgoing (incoming) momentum.

In the CM frame for infinite volume, the separable form of the potential in
Eq.(15) allows for a closed-form partial-wave solution to Eq.(5) (see also
Refs. Guo (2020d); Guo and Gasparian (2021)):

$\displaystyle-\frac{2\mu
q^{2L+1}}{(4\pi)^{2}}\left[\cot\delta_{L}(q)-i\right]$
$\displaystyle=\frac{1}{V_{L}}-\frac{2^{2L+1}\Gamma^{2}(L+\frac{3}{2})}{(2\pi)^{3}}\frac{G^{(\infty)}_{L}(r,r^{\prime};q)}{(rr^{\prime})^{L}}|_{r,r^{\prime}\rightarrow
0}\,,$ (16)

where

$G^{(\infty)}_{L}(r,r^{\prime};q)=-i2\mu qj_{L}(qr_{<})h_{L}^{(+)}(qr_{>})$

is the free-particle Green’s function in the $L$-th partial wave.

In the harmonic trap, the separable form also facilitates a closed-form
solution to Eq.(13):

$\frac{1}{V_{L}}=\frac{2^{2L+1}\Gamma^{2}(L+\frac{3}{2})}{(2\pi)^{3}}\frac{G^{(\omega)}_{L}(r,r^{\prime};\epsilon^{(\omega)})}{(rr^{\prime})^{L}}|_{r,r^{\prime}\rightarrow
0}\,,$ (17)

where $G^{(\omega)}_{L}$ is the partial-wave HO Green’s function and is given
in Ref. Blinder (1984) by

$\displaystyle
G^{(\omega)}_{L}(r,r^{\prime};\epsilon^{(\omega)})=-\frac{1}{\omega(rr^{\prime})^{\frac{3}{2}}}\frac{\Gamma(\frac{L}{2}+\frac{3}{4}-\frac{\epsilon^{(\omega)}}{2\omega})}{\Gamma(L+\frac{3}{2})}$
$\displaystyle\times\mathcal{M}_{\frac{\epsilon^{(\omega)}}{2\omega},\frac{L}{2}+\frac{1}{4}}(\mu\omega
r^{2}_{<})\mathcal{W}_{\frac{\epsilon^{(\omega)}}{2\omega},\frac{L}{2}+\frac{1}{4}}(\mu\omega
r^{2}_{>})\,.$

$\mathcal{M}$ and $\mathcal{W}$ are the Whittaker functions as defined in Ref.
Whittaker and Watson (1996). $\epsilon^{(\omega)}$ represents the energy shift
relative to the HO levels and is related to relative momentum $q$ by

$\frac{q^{2}}{2\mu}=\epsilon^{(\omega)}+\omega(n+\frac{3}{2})=E\,.$

where $\omega(n+\frac{3}{2})$ is the CM energy in the harmonic trap, and the
CM frame in infinite volume is assumed.

By eliminating $V_{L}$ from Eqs. (16) and (17), and using asymptotic forms of
partial-wave Green’s functions $G^{(\omega)}_{L}$ and $G^{(\infty)}_{L}$, one
obtains the BERW formula:

$\cot\delta_{L}(q)=g_{L}(q),$ (18)

where

$g_{L}(q)=(-1)^{L+1}\left(\frac{4\mu\omega}{q^{2}}\right)^{L+\frac{1}{2}}\frac{\Gamma(\frac{3}{4}+\frac{L}{2}-\frac{\epsilon^{(\omega)}}{2\omega})}{\Gamma(\frac{1}{4}-\frac{L}{2}-\frac{\epsilon^{(\omega)}}{2\omega})}.$
(19)

##  III Two-channel reactions

For nuclear reactions where two nuclei fuse, break up, and rearrange, the BERW
formula can not be immediately applied. Extension of the derivation in Sec. II
is needed to extract inelastic scattering amplitudes. We focus on reactions
that involve two channels and assume that each channel is made up of two
distinguishable particles $A+B$ and $C+D$. Labelling two channels by $1$ and
$2$, we can write, for instance, reaction $1\to 2$ as

$A+B\to C+D\,.$ (20)

The kinetic energies in both channels are related by

$\frac{\mathbf{q}_{1}^{2}}{2\mu_{1}}+\frac{\mathbf{P}^{2}}{2M_{1}}=\frac{\mathbf{q}_{2}^{2}}{2\mu_{2}}+\frac{\mathbf{P}^{2}}{2M_{2}}+\Delta=E\,,$
(21)

where $\mathbf{P}$ is the total momentum, $\mathbf{q}_{1,2}$ relative momenta,
$M_{1,2}$ the total mass of each channel,

$\Delta\equiv M_{2}-M_{1},$

and $\mu_{1,2}$ the reduced mass. For simplicity, we will assume the CM frame
in what follows: $\mathbf{P}=0$.

###  III.1 Coupled-channel energy levels in a harmonic trap

We start with the coupled Schrödinger equations for two interacting spinless
particles bound by the HO potential:

$\left[\hat{H}_{\omega}+V(\mathbf{r})\right]\begin{bmatrix}\psi^{(1,\omega)}(\mathbf{r})\\\
\psi^{(2,\omega)}(\mathbf{r})\end{bmatrix}=\begin{bmatrix}\epsilon^{(\omega)}_{1}&0\\\
0&\epsilon^{(\omega)}_{2}\end{bmatrix}\begin{bmatrix}\psi^{(1,\omega)}(\mathbf{r})\\\
\psi^{(2,\omega)}(\mathbf{r})\end{bmatrix},$ (22)

where $\hat{H}_{\omega}$ is given by

$\hat{H}_{\omega}=\begin{bmatrix}-\frac{\nabla_{\mathbf{r}}^{2}}{2\mu_{1}}+\frac{1}{2}\mu_{1}\omega^{2}\mathbf{r}^{2}&0\\\
0&-\frac{\nabla_{\mathbf{r}}^{2}}{2\mu_{2}}+\frac{1}{2}\mu_{2}\omega^{2}\mathbf{r}^{2}\end{bmatrix}.$
(23)

Here $\mathbf{r}$ is the relative coordinate and the potential matrix
$V(\mathbf{r})$ is short-ranged and describes inter-cluster interactions and
transition between the channels:

$V(\mathbf{r})=\begin{bmatrix}V^{(1,1)}(\mathbf{r})&V^{(1,2)}(\mathbf{r})\\\
V^{(2,1)}(\mathbf{r})&V^{(2,2)}(\mathbf{r})\end{bmatrix}.$ (24)

$\epsilon^{(\omega)}_{1}$ and $\epsilon^{(\omega)}_{2}$ are energy shifts in
their respective channel and are related to the full energy $E$ by

$\epsilon^{(\omega)}_{1}+\omega(n_{1}+\frac{3}{2})=\epsilon^{(\omega)}_{2}+\omega(n_{2}+\frac{3}{2})+\Delta=E\,.$

The integral representation of Eq.(22) is the coupled LS equation:

$\begin{bmatrix}\psi^{(1,\omega)}(\mathbf{r})\\\
\psi^{(2,\omega)}(\mathbf{r})\end{bmatrix}=\int
d\mathbf{r}^{\prime}G^{(\omega)}(\mathbf{r},\mathbf{r}^{\prime};E)V(\mathbf{r}^{\prime})\begin{bmatrix}\psi^{(1,\omega)}(\mathbf{r}^{\prime})\\\
\psi^{(2,\omega)}(\mathbf{r}^{\prime})\end{bmatrix},$ (25)

where $G^{(\omega)}$ denotes the coupled-channel version of the HO Green’s
function:

$\displaystyle
G^{(\omega)}(\mathbf{r},\mathbf{r}^{\prime};E)=\begin{bmatrix}G^{(1,\omega)}(\mathbf{r},\mathbf{r}^{\prime};\epsilon^{(\omega)}_{1})&0\\\
0&G^{(2,\omega)}(\mathbf{r},\mathbf{r}^{\prime};\epsilon^{(\omega)}_{2})\end{bmatrix},$
(26)

and $G^{(\alpha,\omega)}$ satisfies

$\left[\epsilon^{(\omega)}_{\alpha}+\frac{\nabla_{\mathbf{r}}^{2}}{2\mu_{\alpha}}-\frac{1}{2}\mu_{\alpha}\omega^{2}\mathbf{r}^{2}\right]G^{(\alpha,\omega)}(\mathbf{r},\mathbf{r}^{\prime};\epsilon^{(\omega)}_{\alpha})=\delta(\mathbf{r}-\mathbf{r}^{\prime}).$
(27)

The partial-wave projection of the HO Green’s function is similar to the
uncoupled-channel case:

$G^{(\alpha,\omega)}(\mathbf{r},\mathbf{r}^{\prime};\epsilon^{(\omega)}_{\alpha})=\sum_{L}\frac{2L+1}{4\pi}G^{(\alpha,\omega)}_{L}(r,r^{\prime};\epsilon^{(\omega)}_{\alpha})P_{L}(\mathbf{\hat{r}}\cdot\mathbf{\hat{r}}^{\prime}),$
(28)

where

$\displaystyle
G^{(\alpha,\omega)}_{L}(r,r^{\prime};\epsilon^{(\omega)}_{\alpha})=-\frac{1}{\omega(rr^{\prime})^{\frac{3}{2}}}\frac{\Gamma(\frac{L}{2}+\frac{3}{4}-\frac{\epsilon^{(\omega)}_{\alpha}}{2\omega})}{\Gamma(L+\frac{3}{2})}$
$\displaystyle\times\mathcal{M}_{\frac{\epsilon^{(\omega)}_{\alpha}}{2\omega},\frac{L}{2}+\frac{1}{4}}(\mu_{\alpha}\omega
r^{2}_{<})\mathcal{W}_{\frac{\epsilon^{(\omega)}_{\alpha}}{2\omega},\frac{L}{2}+\frac{1}{4}}(\mu_{\alpha}\omega
r^{2}_{>}).$ (29)

We are going to take one more step to write the integral equation (25) in
momentum space:

$T^{(\omega)}(\mathbf{q}^{\prime})=\int\frac{d\mathbf{q}}{(2\pi)^{3}}\frac{d\mathbf{k}}{(2\pi)^{3}}\widetilde{V}(\mathbf{q}^{\prime}-\mathbf{k})\widetilde{G}^{(\omega)}(\mathbf{k},\mathbf{q};E)T^{(\omega)}(\mathbf{q}),$
(30)

where $T^{(\omega)}(\mathbf{q}^{\prime})$ matrix is defined by

$T^{(\omega)}(\mathbf{q}^{\prime})=-\int
d\mathbf{r}e^{-i\mathbf{q}^{\prime}\cdot\mathbf{r}}V(\mathbf{r})\begin{bmatrix}\psi^{(1,\omega)}(\mathbf{r})\\\
\psi^{(2,\omega)}(\mathbf{r})\end{bmatrix}.$ (31)

$\widetilde{V}$ and $\widetilde{G}^{(\omega)}$ are the Fourier transform of
$V$ and $G^{(\omega)}$, respectively. The partial wave expansion of Eq.(30)
yields

$T^{(\omega)}_{L}(q^{\prime})=\int_{0}^{\infty}\frac{q^{2}dq}{(2\pi)^{3}}\frac{k^{2}dk}{(2\pi)^{3}}\widetilde{V}_{L}(q^{\prime},k)\widetilde{G}^{(\omega)}_{L}(k,q;E)T^{(\omega)}_{L}(q).$
(32)

With the assumption of separable potential,

$\widetilde{V}_{L}(q^{\prime},q)=(q^{\prime}q)^{L}V_{L}\,,$ (33)

where $V_{L}$ is a constant $2\times 2$ matrix, it follows that
$T^{(\omega)}_{L}(q^{\prime})$ matrix must be of the form

$T^{(\omega)}_{L}(q)=q^{L}t^{(\omega)}_{L}\,,$

where $t^{(\omega)}_{L}$ matrix does not depend on $q$. Therefore, the
quantization condition is given by

$\det\left[V_{L}^{-1}-\int_{0}^{\infty}\frac{q^{2}dq}{(2\pi)^{3}}\frac{k^{2}dk}{(2\pi)^{3}}(kq)^{L}\widetilde{G}^{(\omega)}_{L}(k,q;E)\right]=0.$
(34)

Using the limiting form of the spherical Bessel functions

$q^{L}=\frac{2^{L+1}\Gamma(L+\frac{3}{2})}{\sqrt{\pi}}\frac{j_{L}(qr)}{r^{L}}|_{r\rightarrow
0}\,,$ (35)

one can show that

$\displaystyle\int_{0}^{\infty}\frac{q^{2}dq}{(2\pi)^{3}}\frac{k^{2}dk}{(2\pi)^{3}}(kq)^{L}\widetilde{G}^{(\omega)}_{L}(k,q;E)$
$\displaystyle=\frac{2^{2L+1}\Gamma^{2}(L+\frac{3}{2})}{(2\pi)^{3}}\frac{G^{(\omega)}_{L}(r,r^{\prime};E)}{(rr^{\prime})^{L}}|_{r,r^{\prime}\rightarrow
0}\,.$ (36)

Therefore, we can rewrite Eq.(34) as

$\det\left[V_{L}^{-1}-\frac{2^{2L+1}\Gamma^{2}(L+\frac{3}{2})}{(2\pi)^{3}}\frac{G^{(\omega)}_{L}(r,r^{\prime};E)}{(rr^{\prime})^{L}}|_{r,r^{\prime}\rightarrow
0}\right]=0,$ (37)

which is just the generalization of Eq.(17) in the case of coupled-channel
reactions.

###  III.2  Coupled-channel reaction amplitude

In infinite volume, the coupled-channel reaction is described by the two-
channel LS equation:

$\displaystyle
T^{(\infty)}(\mathbf{q}^{\prime},\mathbf{q})=-\widetilde{V}(\mathbf{q}^{\prime}-\mathbf{q})$
$\displaystyle+\int\frac{d\mathbf{k}}{(2\pi)^{3}}\widetilde{V}(\mathbf{q}^{\prime}-\mathbf{k})\widetilde{G}^{(\infty)}(\mathbf{k};E)T^{(\infty)}(\mathbf{k},\mathbf{q})\,,$
(38)

where $T^{(\infty)}$ is a $2\times 2$ matrix,

$T^{(\infty)}=\begin{bmatrix}T^{((1,1),\infty)}&T^{((1,2),\infty)}\\\
T^{((2,1),\infty)}&T^{((2,2),\infty)}\end{bmatrix}\,,$ (39)

and the free Green’s function is also in matrix form,

$\widetilde{G}^{(\infty)}(\mathbf{k};E)=\begin{bmatrix}\frac{1}{\frac{q_{1}^{2}}{2\mu_{1}}-\frac{\mathbf{k}^{2}}{2\mu_{1}}}&0\\\
0&\frac{1}{\frac{q_{2}^{2}}{2\mu_{2}}-\frac{\mathbf{k}^{2}}{2\mu_{2}}}\end{bmatrix},$
(40)

with

$\frac{q_{1}^{2}}{2\mu_{1}}=\frac{q_{2}^{2}}{2\mu_{2}}+\Delta=E\,.$

The partial-wave expansion of Eq.(38) produces

$\displaystyle
T_{L}^{(\infty)}(q^{\prime},q)=-\widetilde{V}_{L}(q^{\prime},q)$
$\displaystyle+\int_{0}^{\infty}\frac{k^{2}dk}{(2\pi)^{3}}\widetilde{V}_{L}(q^{\prime},k)\widetilde{G}^{(\infty)}(k;E)T_{L}^{(\infty)}(k,q).$
(41)

The separable form of the potential (33) enables semi-closed form solutions to
the above integration equation:

$T_{L}^{(\infty)}(q^{\prime},q)=(q^{\prime}q)^{L}t_{L}(E),$ (42)

where

$t_{L}^{-1}(E)=-V_{L}^{-1}+\int_{0}^{\infty}\frac{k^{2}dk}{(2\pi)^{3}}k^{2L}\widetilde{G}^{(\infty)}(k;E).$
(43)

When it is on-shell, the matrix $t_{L}$ is usually parameterized by phase
shifts, $\delta_{L}^{(1,2)}$, and inelasticity $\eta_{L}$,

$\frac{t_{L}(E)}{(4\pi)^{2}}=\begin{bmatrix}\
\frac{\frac{1}{2i}\left(\eta_{L}e^{2i\delta_{L}^{(1)}}-1\right)}{2\mu_{1}q_{1}^{2L+1}}&\frac{\frac{\sqrt{1-\eta^{2}_{L}}}{2}e^{i\left(\delta_{L}^{(1)}+\delta_{L}^{(2)}\right)}}{2\mu_{1}q_{1}(q_{1}q_{2})^{L}}\\\
\frac{\frac{\sqrt{1-\eta^{2}_{L}}}{2}e^{i\left(\delta_{L}^{(1)}+\delta_{L}^{(2)}\right)}}{2\mu_{2}q_{2}(q_{1}q_{2})^{L}}&\frac{\frac{1}{2i}\left(\eta_{L}e^{2i\delta_{L}^{(2)}}-1\right)}{2\mu_{2}q_{2}^{2L+1}}\end{bmatrix}.$
(44)

Using the following identity,

$\displaystyle\int_{0}^{\infty}\frac{k^{2}dk}{(2\pi)^{3}}k^{2L}\widetilde{G}^{(\infty)}(k;E)$
$\displaystyle=\frac{2^{2L+1}\Gamma^{2}(L+\frac{3}{2})}{(2\pi)^{3}}\frac{G^{(\infty)}_{L}(r,r^{\prime};E)}{(rr^{\prime})^{L}}|_{r,r^{\prime}\rightarrow
0},$ (45)

where

$\displaystyle
G^{(\infty)}_{L}(r,r^{\prime};E)=\begin{bmatrix}G^{(1,\infty)}_{L}(r,r^{\prime};q_{1})&0\\\
0&G^{(2,\infty)}_{L}(r,r^{\prime};q_{2})\end{bmatrix},$ $\displaystyle
G^{(\alpha,\infty)}_{L}(r,r^{\prime};q_{\alpha})=-i2\mu_{\alpha}q_{\alpha}j_{L}(q_{\alpha}r_{<})h_{L}^{(+)}(q_{\alpha}r_{>})\,,$
(46)

we find the coupled-channel version of Eq.(16):

$t_{L}^{-1}(E)=-V_{L}^{-1}+\frac{2^{2L+1}\Gamma^{2}(L+\frac{3}{2})}{(2\pi)^{3}}\frac{G^{(\infty)}_{L}(r,r^{\prime};E)}{(rr^{\prime})^{L}}|_{r,r^{\prime}\rightarrow
0}.$ (47)

###  III.3 BERW formula for coupled-channel reactions

Using Eq.(47) in Eq.(37) and , one finds

$\displaystyle\det$
$\displaystyle\bigg{[}\frac{2^{2L+2}\Gamma^{2}(L+\frac{3}{2})}{\pi}\frac{G^{(\omega)}_{L}(r,r^{\prime};E)-G^{(\infty)}_{L}(r,r^{\prime};E)}{(rr^{\prime})^{L}}|_{r,r^{\prime}\rightarrow
0}$ $\displaystyle+(4\pi)^{2}t_{L}^{-1}(E)\bigg{]}=0\,.$ (48)

The cancellation of ultraviolet divergences manifested by the Green’s
functions is made clear if we use the limiting forms of both Green’s
functions:

$\displaystyle\frac{2^{2L+2}\Gamma(L+\frac{3}{2})^{2}}{\pi}\frac{G^{(\alpha,\omega)}_{L}(r,r^{\prime};\epsilon^{(\omega)}_{\alpha})}{(rr^{\prime})^{L}}|_{r,r^{\prime}\rightarrow
0}$
$\displaystyle=-\frac{(\mu_{\alpha}\omega)^{L+\frac{3}{2}}}{\omega}2^{2L+2}(-1)^{L+1}\frac{\Gamma(\frac{L}{2}+\frac{3}{4}-\frac{\epsilon^{(\omega)}_{\alpha}}{2\omega})}{\Gamma(\frac{1}{4}-\frac{L}{2}-\frac{\epsilon^{(\omega)}_{\alpha}}{2\omega})}$
$\displaystyle-\frac{2^{2L+1}\Gamma(L+\frac{1}{2})\Gamma(L+\frac{3}{2})}{\pi}\frac{2\mu_{\alpha}}{r^{2L+1}}|_{r\rightarrow
0},$ (49)

and

$\displaystyle\frac{2^{2L+2}\Gamma(L+\frac{3}{2})^{2}}{\pi}\frac{G_{L}^{(\alpha,\infty)}(r,r^{\prime};q_{\alpha})}{(rr^{\prime})^{L}}|_{r,r^{\prime}\rightarrow
0}$
$\displaystyle=-i2\mu_{\alpha}q_{\alpha}^{2L+1}-\frac{2^{2L+1}\Gamma(L+\frac{1}{2})\Gamma(L+\frac{3}{2})}{\pi}\frac{2\mu_{\alpha}}{r^{2L+1}}|_{r\rightarrow
0}\,.$ (50)

Equation (48) is thus reduced to the more compact form:

$\displaystyle\quad\eta_{L}\left(1+g_{L}^{(1)}g_{L}^{(2)}\right)\cos\left(\delta_{L}^{(1)}-\delta_{L}^{(2)}\right)$
$\displaystyle+\left(1-g_{L}^{(1)}g_{L}^{(2)}\right)\cos\left(\delta_{L}^{(1)}+\delta_{L}^{(2)}\right)$
$\displaystyle-\eta_{L}\left(g_{L}^{(1)}-g_{L}^{(2)}\right)\sin\left(\delta_{L}^{(1)}-\delta_{L}^{(2)}\right)$
$\displaystyle-\left(g_{L}^{(1)}+g_{L}^{(2)}\right)\sin\left(\delta_{L}^{(1)}+\delta_{L}^{(2)}\right)=0,$
(51)

where

$g_{L}^{(\alpha)}=(-1)^{L+1}\left(\frac{4\mu_{\alpha}\omega}{q_{\alpha}^{2}}\right)^{L+\frac{1}{2}}\frac{\Gamma(\frac{L}{2}+\frac{3}{4}-\frac{\epsilon^{(\omega)}_{\alpha}}{2\omega})}{\Gamma(\frac{1}{4}-\frac{L}{2}-\frac{\epsilon^{(\omega)}_{\alpha}}{2\omega})}.$
(52)

The coupled-channel QC in a harmonic trap expressed by Eq.(51) resembles that
in a periodic cube [see Eq.(25) in Guo et al. (2013)]. In the harmonic trap,
$g_{L}^{(\alpha)}$ functions play the same role of Lüscher zeta functions in
the periodic cube.

Given one energy eigenvalue $E$, we can not determine all of
$\delta_{L}^{(1,2)}$ and $\eta_{L}$ simultaneously. A possible way out is to
gather several pairs of $(E,\omega)$ to build more constraints on
$\delta_{L}^{(1,2)}$ and $\eta_{L}$, provided that energy eigenvalues $E$’s
are somewhat close to each other and the ratio $E/\omega$ is distinct from one
another. Though this approach is model- independent, the requirement may only
be met by small portion of the calculated levels. Another useful approach
widely used in the lattice QCD community is to model $\delta_{L}^{(1,2)}$ and
$\eta_{L}$ with a few free parameters, which can be fitted to more energy
levels by using QC as defined in Eq.(51), (see, e.g., discussions in Refs. Guo
et al. (2013); Guo (2013)).

##  IV Discussions and Summary

The harmonic trap may be replaced with artificial fields of different
configuration or boundary condition. For instance, hard-sphere boundary
condition is used in some Monte Carlo simulations Bovermann et al. (2019);
Elhatisari et al. (2016); Rokash et al. (2015). We can repeat all the
discussions in previous sections by substituting the following external
potential in place of the HO potential:

$V_{h.s.}(r)=\begin{cases}0,&r<b\,,\\\ \infty,&r>b\,,\end{cases}$ (53)

where $b$ is the radius of the sphere. The Green’s function must satisfy the
boundary condition:

$\displaystyle G_{L}^{(h.s.)}(r,r^{\prime};q)$ $\displaystyle=-2\mu
qj_{L}(qr_{<})j_{L}(qr_{>})\left[\frac{n_{L}(qb)}{j_{L}(qb)}-\frac{n_{L}(qr_{>})}{j_{L}(qr_{>})}\right].$
(54)

For instance, in the case of the single-channel elastic scattering, Eqs. (16)
and (17) are replaced with

$\displaystyle\frac{2^{2L+2}\Gamma^{2}(L+\frac{3}{2})}{\pi}\frac{G^{(h.s.)}_{L}(r,r^{\prime};q)-G^{(\infty)}_{L}(r,r^{\prime};q)}{(rr^{\prime})^{L}}|_{r,r^{\prime}\rightarrow
0}$ $\displaystyle=-2\mu q^{2L+1}[\cot\delta_{L}(q)-i].$ (55)

Using the asymptotic form of $G^{(\infty)}_{L}$ (50) and

$\displaystyle\frac{1}{2\mu}\frac{2^{2L+2}\Gamma^{2}(L+\frac{3}{2})}{\pi}\frac{G^{(h.s.)}_{L}(r,r^{\prime};q)}{(rr^{\prime})^{L}}|_{r,r^{\prime}\rightarrow
0}$
$\displaystyle=-q^{2L+1}\frac{n_{L}(qb)}{j_{L}(qb)}-\frac{2^{2L+1}\Gamma(L+\frac{1}{2})\Gamma(L+\frac{3}{2})}{\pi}\frac{1}{r^{2L+1}}|_{r\rightarrow
0},$ (56)

the quantization condition of the single-channel case for the spherical hard
wall is thus given by a form analogous to the BERW formula (or Lüscher’s for
that matter):

$\displaystyle\cot\delta_{L}(q)=\frac{n_{L}(qb)}{j_{L}(qb)}.$ (57)

The hard-sphere geometry is encoded in $n_{L}(qb)/j_{L}(qb)$, which plays the
same role as $g_{L}(q)$ to the harmonic trap. Similarly, the coupled-channel
QC can be worked out for the spherical hard wall by replacing
$g_{L}^{(\alpha)}$ in Eq.(52) as follows

$g_{L}^{(\alpha)}\rightarrow\frac{n_{L}(q_{\alpha}b)}{j_{L}(q_{\alpha}b)}\,.$
(58)

In summary, we have derived the quantization condition for coupled-channel
two-particle systems trapped by the harmonic potential. The formula by itself
relates the coupled-channel reaction amplitudes in infinite volume,
parametrized by two phase shifts and one inelasticity, to energy eigenvalues
of the corresponding trapped system. Although only the quantization condition
for the case of distinguishable spinless particles is shown, with the
technical detail spelled out in the paper, the extension to include slight
complications due to spin, identical particles and so on should be
straightforward. Coulomb repulsion is important, however, if reactions take
place near the threshold of charged particles. Quantization conditions may be
derived for that circumstance in the harmonic trap, following the methodology
developed in Ref. Guo and Gasparian (2021).

###### Acknowledgements.

B.L. thanks Sebastian König and Baishan Hu for useful discussions at early
stage of this work. P.G. acknowledges support from the Department of Physics
and Engineering, California State University, Bakersfield, CA. The work (BL)
was supported in part by the National Natural Science Foundation of China
(NSFC) under Grants No.11775148 and No.11735003, and (PG) was also supported
in part by the National Science Foundation (US) under Grant No. NSF
PHY-1748958.

## References

  * Navratil et al. (2007) P. Navratil, V. Gueorguiev, J. Vary, W. Ormand, and A. Nogga, Phys. Rev. Lett. 99, 042501 (2007), eprint nucl-th/0701038.
  * Epelbaum et al. (2009) E. Epelbaum, H.-W. Hammer, and U.-G. Meissner, Rev. Mod. Phys. 81, 1773 (2009), eprint 0811.1338.
  * Ekström et al. (2013) A. Ekström et al., Phys. Rev. Lett. 110, 192502 (2013), eprint 1303.4674.
  * Barrett et al. (2013) B. R. Barrett, P. Navratil, and J. P. Vary, Prog. Part. Nucl. Phys. 69, 131 (2013).
  * Lynn et al. (2014) J. Lynn, J. Carlson, E. Epelbaum, S. Gandolfi, A. Gezerlis, and A. Schwenk, Phys. Rev. Lett. 113, 192501 (2014), eprint 1406.2787.
  * Ekström et al. (2015) A. Ekström, G. Jansen, K. Wendt, G. Hagen, T. Papenbrock, B. Carlsson, C. Forssén, M. Hjorth-Jensen, P. Navrátil, and W. Nazarewicz, Phys. Rev. C 91, 051301 (2015), eprint 1502.04682.
  * Piarulli et al. (2018) M. Piarulli et al., Phys. Rev. Lett. 120, 052503 (2018), eprint 1707.02883.
  * Hammer et al. (2020) H.-W. Hammer, S. König, and U. van Kolck, Rev. Mod. Phys. 92, 025004 (2020), eprint 1906.12122.
  * Navrátil et al. (2016) P. Navrátil, S. Quaglioni, G. Hupin, C. Romero-Redondo, and A. Calci, Phys. Scripta 91, 053002 (2016), eprint 1601.03765.
  * Yang (2016) C.-J. Yang, Phys. Rev. C 94, 064004 (2016), eprint 1610.01350.
  * Johnson et al. (2019) C. W. Johnson et al., in _From Bound States to the Continuum: Connecting bound state calculations with scattering and reaction theory_ (2019), eprint 1912.00451.
  * Mazur et al. (2020) I. Mazur, A. Shirokov, I. Shin, A. Mazur, Y. Kim, P. Maris, and J. Vary (2020), eprint 2001.08898.
  * Ma et al. (2020) Y. Ma, F. Xu, N. Michel, S. Zhang, J. Li, B. Hu, L. Coraggio, N. Itaco, and A. Gargano, Phys. Lett. B 808, 135673 (2020), eprint 2008.01420.
  * Rotureau et al. (2010) J. Rotureau, I. Stetcu, B. Barrett, M. Birse, and U. van Kolck, Phys. Rev. A 82, 032711 (2010), eprint 1006.3820.
  * Rotureau et al. (2012) J. Rotureau, I. Stetcu, B. Barrett, and U. van Kolck, Phys. Rev. C 85, 034003 (2012), eprint 1112.0267.
  * Luu et al. (2010) T. Luu, M. J. Savage, A. Schwenk, and J. P. Vary, Phys. Rev. C 82, 034003 (2010), eprint 1006.0427.
  * Zhang et al. (2020) X. Zhang, S. Stroberg, P. Navrátil, C. Gwak, J. Melendez, R. Furnstahl, and J. Holt, Phys. Rev. Lett. 125, 112503 (2020), eprint 2004.13575.
  * Busch et al. (1998) T. Busch, B.-G. Englert, K. Rzażewski, and M. Wilkens, Found. Phys. 28, 549–559 (1998).
  * Zhang (2020) X. Zhang, Phys. Rev. C 101, 051602 (2020), eprint 1905.05275.
  * Navratil et al. (2000) P. Navratil, J. P. Vary, and B. R. Barrett, Phys. Rev. C 62, 054311 (2000).
  * Stetcu et al. (2007) I. Stetcu, B. Barrett, U. van Kolck, and J. Vary, Phys. Rev. A 76, 063613 (2007), eprint 0705.4335.
  * Stetcu et al. (2010) I. Stetcu, J. Rotureau, B. Barrett, and U. van Kolck, Annals Phys. 325, 1644 (2010), eprint 1001.5071.
  * Lüscher (1991) M. Lüscher, Nucl. Phys. B354, 531 (1991).
  * Bovermann et al. (2019) L. Bovermann, E. Epelbaum, H. Krebs, and D. Lee, Phys. Rev. C 100, 064001 (2019), eprint 1905.02492.
  * Elhatisari et al. (2016) S. Elhatisari, D. Lee, U.-G. Meißner, and G. Rupak, Eur. Phys. J. A 52, 174 (2016), eprint 1603.02333.
  * Rokash et al. (2015) A. Rokash, M. Pine, S. Elhatisari, D. Lee, E. Epelbaum, and H. Krebs, Phys. Rev. C 92, 054612 (2015), eprint 1505.02967.
  * Rummukainen and Gottlieb (1995) K. Rummukainen and S. A. Gottlieb, Nucl. Phys. B450, 397 (1995), eprint hep-lat/9503028.
  * Christ et al. (2005) N. H. Christ, C. Kim, and T. Yamazaki, Phys. Rev. D72, 114506 (2005), eprint hep-lat/0507009.
  * Bernard et al. (2008) V. Bernard, M. Lage, U.-G. Meißner, and A. Rusetsky, JHEP 08, 024 (2008), eprint 0806.4495.
  * He et al. (2005) S. He, X. Feng, and C. Liu, JHEP 07, 011 (2005), eprint hep-lat/0504019.
  * Lage et al. (2009) M. Lage, U.-G. Meißner, and A. Rusetsky, Phys. Lett. B681, 439 (2009), eprint 0905.0069.
  * Döring et al. (2011) M. Döring, U.-G. Meißner, E. Oset, and A. Rusetsky, Eur. Phys. J. A47, 139 (2011), eprint 1107.3988.
  * Guo et al. (2013) P. Guo, J. Dudek, R. Edwards, and A. P. Szczepaniak, Phys. Rev. D88, 014501 (2013), eprint 1211.0929.
  * Guo (2013) P. Guo, Phys. Rev. D88, 014507 (2013), eprint 1304.7812.
  * Kreuzer and Hammer (2009) S. Kreuzer and H. W. Hammer, Phys. Lett. B673, 260 (2009), eprint 0811.0159.
  * Polejaeva and Rusetsky (2012) K. Polejaeva and A. Rusetsky, Eur. Phys. J. A48, 67 (2012), eprint 1203.1241.
  * Hansen and Sharpe (2014) M. T. Hansen and S. R. Sharpe, Phys. Rev. D90, 116003 (2014), eprint 1408.5933.
  * Mai and Döring (2017) M. Mai and M. Döring, Eur. Phys. J. A53, 240 (2017), eprint 1709.08222.
  * Mai and Döring (2019) M. Mai and M. Döring, Phys. Rev. Lett. 122, 062503 (2019), eprint 1807.04746.
  * Döring et al. (2018) M. Döring, H.-W. Hammer, M. Mai, J.-Y. Pang, A. Rusetsky, and J. Wu, Phys. Rev. D 97, 114508 (2018), eprint 1802.03362.
  * Guo (2017) P. Guo, Phys. Rev. D95, 054508 (2017), eprint 1607.03184.
  * Guo and Gasparian (2017) P. Guo and V. Gasparian, Phys. Lett. B774, 441 (2017), eprint 1701.00438.
  * Guo and Gasparian (2018) P. Guo and V. Gasparian, Phys. Rev. D97, 014504 (2018), eprint 1709.08255.
  * Guo and Morris (2019) P. Guo and T. Morris, Phys. Rev. D99, 014501 (2019), eprint 1808.07397.
  * Guo et al. (2018) P. Guo, M. Döring, and A. P. Szczepaniak, Phys. Rev. D98, 094502 (2018), eprint 1810.01261.
  * Guo (2020a) P. Guo, Phys. Lett. B 804, 135370 (2020a), eprint 1908.08081.
  * Guo and Döring (2020) P. Guo and M. Döring, Phys. Rev. D 101, 034501 (2020), eprint 1910.08624.
  * Guo (2020b) P. Guo, Phys. Rev. D101, 054512 (2020b), eprint 2002.04111.
  * Guo and Long (2020a) P. Guo and B. Long, Phys. Rev. D 101, 094510 (2020a), eprint 2002.09266.
  * Guo (2020c) P. Guo (2020c), eprint 2007.04473.
  * Guo and Long (2020b) P. Guo and B. Long, Phys. Rev. D 102, 074508 (2020b), eprint 2007.10895.
  * Guo (2020d) P. Guo, Phys. Rev. D 102, 054514 (2020d), eprint 2007.12790.
  * Guo and Gasparian (2021) P. Guo and V. Gasparian (2021), eprint 2101.01150.
  * Blinder (1984) S. Blinder, J. Math. Phys. 25, 905 (1984).
  * Whittaker and Watson (1996) E. T. Whittaker and G. N. Watson, _A Course of Modern Analysis_ , Cambridge Mathematical Library (Cambridge University Press, 1996), 4th ed.