[1]\fnmJia \surWang 1]\orgdivCentre for Quantum Technology Theory, \orgnameSwinburne University of Technology, \orgaddress\cityMelbourne, \postcode3122, \countryAustralia

Hyperspherical Analysis of Dimer-Dimer Scattering in One-Dimensional Systems

[email protected]    \fnmHui \surHu    \fnmXia-Ji \surLiu [
Abstract

We present a comprehensive analysis of four-body scattering in one-dimensional (1D) quantum systems using the adiabatic hyperspherical representation (AHR). Focusing on dimer-dimer collisions between two species of fermions interacting via the sinh-cosh potential, we implement the slow variable discretization (SVD) method to overcome numerical challenges posed by sharp avoided crossings in the potential curves. Our numerical approach is benchmarked against exact analytical results available in integrable regimes, demonstrating excellent agreement. We further explore non-integrable regimes where no analytical solutions exist, revealing novel features such as resonant enhancement of the scattering length associated with tetramer formation. These results highlight the power and flexibility of the AHR+SVD framework for accurate few-body scattering calculations in low-dimensional quantum systems, and establish a foundation for future investigations of universal few-body physics in ultracold gases.

1 Introduction

Since the groundbreaking development of laser cooling and trapping of atoms [1], the field of ultracold quantum gases has rapidly expanded and made a profound impact on the study of fundamental quantum phenomena [2, 3, 4]. At temperatures as low as a billionth of a degree Kelvin (nano-Kelvin, nK) and densities ranging from 1012superscript101210^{12}10 start_POSTSUPERSCRIPT 12 end_POSTSUPERSCRIPT to 1015superscript101510^{15}10 start_POSTSUPERSCRIPT 15 end_POSTSUPERSCRIPT cm3superscriptcm3{\rm cm}^{-3}roman_cm start_POSTSUPERSCRIPT - 3 end_POSTSUPERSCRIPT, these systems reach the quantum degenerate regime, where the atomic de Broglie wavelength exceeds the typical interparticle spacing [5, 6]. A variety of trapping techniques enable the realization of low-dimensional systems, while Feshbach resonances allow for precise and tunable control over interparticle interactions [7, 8]. This extraordinary level of control, combined with recent advances in theoretical and computational techniques, makes ultracold gases an ideal platform for exploring few-body quantum physics and performing quantitative comparisons between theory and experiment.

One of the most striking few-body phenomena predicted in this context is the Efimov effect. In the early 1970s, Vitaly Efimov predicted that three identical resonant bosons can form a series of three-body bound states—now known as Efimov trimers—characterized by a discrete scaling symmetry [9]. Remarkably, the scaling factor is independent of short-range interaction details, making the phenomenon universal. Once considered purely theoretical, Efimov states have since been observed in numerous ultracold gas experiments [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]. Moreover, new forms of universality have emerged in these systems, with Efimov trimers displaying properties that reflect the long-range van der Waals tails of the underlying two-body potentials [21, 22, 23, 24, 25]. Other exotic few-body states, including super-Efimov [26, 27, 28] and semi-super Efimov states [29], have been theoretically proposed and are being actively explored in ultracold gases, particularly in systems with tunable dimensionality.

In this work, we investigate the four-body scattering problem in one-dimensional (1D) systems. A system of four fermions—two in one component and two in another—is the simplest nontrivial few-body configuration that captures the essential physics of ultracold fermionic gases. We numerically obtain essentially exact solutions to this problem using the adiabatic hyperspherical representation (AHR). AHR is a well-established method in atomic, molecular, and optical physics, as well as in nuclear and condensed matter contexts. It has been particularly successful in describing few-body scattering in ultracold gases, notably in the seminal work by Greene and collaborators on ultracold collisions [30, 31, 32, 33, 34, 35, 36] and Efimov physics [37, 38, 23]. While AHR can also be applied to four-body problems in three dimensions, doing so typically requires carefully tailored variational basis sets due to the increased number of degrees of freedom [39, 40, 41].

In contrast, the reduced dimensionality of 1D systems offers computational advantages, allowing us to apply numerically exact diagonalization methods using simple B-spline basis functions. Although some adiabatic hyperspherical potential curves for 1D four-body systems have been presented in previous studies [39], scattering observables such as the dimer-dimer scattering length have not yet been computed. A key challenge in these calculations is the presence of sharp avoided crossings in the adiabatic potential curves, which complicate numerical integration of the coupled channel equations.

To overcome this difficulty, we employ the slow variable discretization (SVD) method, which enables stable and accurate computation of scattering observables [42, 43]. In particular, we focus on a model interaction known as the sinh\sinhroman_sinh-cosh\coshroman_cosh potential, for which analytical solutions are available for certain parameter regimes [44]. We benchmark our numerical results against these exact solutions, and further extend the analysis to parameter ranges where no analytical solutions exist.

2 Theory

2.1 Hyperspherical Coordinate

We consider four particles confined to one dimension (1D), with coordinates rjsubscript𝑟𝑗r_{j}italic_r start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT and masses mjsubscript𝑚𝑗m_{j}italic_m start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT (for j=1,2,3,4𝑗1234j=1,2,3,4italic_j = 1 , 2 , 3 , 4). It is convenient to begin by transforming to “H”-type Jacobi coordinates, where the configuration is illustrated in Fig. 1 (a). We define

y1=μ1μ4bx1=μ1μ4b(r1r2),y2=μ2μ4bx2=μ2μ4b(r3r4),y3=μ3μ4bx3=μ3μ4b(r1+r22r3+r42),y4=x4=r1+r2+r3+r44.formulae-sequencesubscript𝑦1subscript𝜇1subscript𝜇4bsubscript𝑥1subscript𝜇1subscript𝜇4bsubscript𝑟1subscript𝑟2subscript𝑦2subscript𝜇2subscript𝜇4bsubscript𝑥2subscript𝜇2subscript𝜇4bsubscript𝑟3subscript𝑟4subscript𝑦3subscript𝜇3subscript𝜇4bsubscript𝑥3subscript𝜇3subscript𝜇4bsubscript𝑟1subscript𝑟22subscript𝑟3subscript𝑟42subscript𝑦4subscript𝑥4subscript𝑟1subscript𝑟2subscript𝑟3subscript𝑟44\begin{gathered}y_{1}=\sqrt{\frac{\mu_{1}}{\mu_{{\rm 4b}}}}x_{1}=\sqrt{\frac{% \mu_{1}}{\mu_{{\rm 4b}}}}\left(r_{1}-r_{2}\right),\\ y_{2}=\sqrt{\frac{\mu_{2}}{\mu_{{\rm 4b}}}}x_{2}=\sqrt{\frac{\mu_{2}}{\mu_{{% \rm 4b}}}}\left(r_{3}-r_{4}\right),\\ y_{3}=\sqrt{\frac{\mu_{3}}{\mu_{{\rm 4b}}}}x_{3}=\sqrt{\frac{\mu_{3}}{\mu_{{% \rm 4b}}}}\left(\frac{r_{1}+r_{2}}{2}-\frac{r_{3}+r_{4}}{2}\right),\\ y_{4}=x_{4}=\frac{r_{1}+r_{2}+r_{3}+r_{4}}{4}.\end{gathered}start_ROW start_CELL italic_y start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = square-root start_ARG divide start_ARG italic_μ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_ARG start_ARG italic_μ start_POSTSUBSCRIPT 4 roman_b end_POSTSUBSCRIPT end_ARG end_ARG italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = square-root start_ARG divide start_ARG italic_μ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_ARG start_ARG italic_μ start_POSTSUBSCRIPT 4 roman_b end_POSTSUBSCRIPT end_ARG end_ARG ( italic_r start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - italic_r start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ) , end_CELL end_ROW start_ROW start_CELL italic_y start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = square-root start_ARG divide start_ARG italic_μ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT end_ARG start_ARG italic_μ start_POSTSUBSCRIPT 4 roman_b end_POSTSUBSCRIPT end_ARG end_ARG italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = square-root start_ARG divide start_ARG italic_μ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT end_ARG start_ARG italic_μ start_POSTSUBSCRIPT 4 roman_b end_POSTSUBSCRIPT end_ARG end_ARG ( italic_r start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT - italic_r start_POSTSUBSCRIPT 4 end_POSTSUBSCRIPT ) , end_CELL end_ROW start_ROW start_CELL italic_y start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT = square-root start_ARG divide start_ARG italic_μ start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT end_ARG start_ARG italic_μ start_POSTSUBSCRIPT 4 roman_b end_POSTSUBSCRIPT end_ARG end_ARG italic_x start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT = square-root start_ARG divide start_ARG italic_μ start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT end_ARG start_ARG italic_μ start_POSTSUBSCRIPT 4 roman_b end_POSTSUBSCRIPT end_ARG end_ARG ( divide start_ARG italic_r start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT + italic_r start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG - divide start_ARG italic_r start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT + italic_r start_POSTSUBSCRIPT 4 end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG ) , end_CELL end_ROW start_ROW start_CELL italic_y start_POSTSUBSCRIPT 4 end_POSTSUBSCRIPT = italic_x start_POSTSUBSCRIPT 4 end_POSTSUBSCRIPT = divide start_ARG italic_r start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT + italic_r start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT + italic_r start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT + italic_r start_POSTSUBSCRIPT 4 end_POSTSUBSCRIPT end_ARG start_ARG 4 end_ARG . end_CELL end_ROW (1)

Here, xisubscript𝑥𝑖x_{i}italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT are the conventional Jacobi coordinate, while yisubscript𝑦𝑖y_{i}italic_y start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT are their corresponding mass-scaled forms. The centre-of-mass coordinate y4=x4subscript𝑦4subscript𝑥4y_{4}=x_{4}italic_y start_POSTSUBSCRIPT 4 end_POSTSUBSCRIPT = italic_x start_POSTSUBSCRIPT 4 end_POSTSUBSCRIPT can usually be separated from the dynamics. The reduced masses are defined by μ11=m11+m21superscriptsubscript𝜇11superscriptsubscript𝑚11superscriptsubscript𝑚21\mu_{1}^{-1}=m_{1}^{-1}+m_{2}^{-1}italic_μ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT = italic_m start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT + italic_m start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT, μ21=m31+m41superscriptsubscript𝜇21superscriptsubscript𝑚31superscriptsubscript𝑚41\mu_{2}^{-1}=m_{3}^{-1}+m_{4}^{-1}italic_μ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT = italic_m start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT + italic_m start_POSTSUBSCRIPT 4 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT, and μ31=(m1+m2)1+(m3+m4)1superscriptsubscript𝜇31superscriptsubscript𝑚1subscript𝑚21superscriptsubscript𝑚3subscript𝑚41\mu_{3}^{-1}=\left(m_{1}+m_{2}\right)^{-1}+\left(m_{3}+m_{4}\right)^{-1}italic_μ start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT = ( italic_m start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT + italic_m start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT + ( italic_m start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT + italic_m start_POSTSUBSCRIPT 4 end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT. The four-body reduced mass μ4bsubscript𝜇4b\mu_{{\rm 4b}}italic_μ start_POSTSUBSCRIPT 4 roman_b end_POSTSUBSCRIPT can essentially be defined arbitrarily as long as the relative kinetic energy operator are defined consistently. We choose μ4b=Πjmj/(2jmj)3subscript𝜇4b3subscriptΠ𝑗subscript𝑚𝑗2subscript𝑗subscript𝑚𝑗\mu_{{\rm 4b}}=\sqrt[3]{\Pi_{j}m_{j}/\left(2\sum_{j}m_{j}\right)}italic_μ start_POSTSUBSCRIPT 4 roman_b end_POSTSUBSCRIPT = nth-root start_ARG 3 end_ARG start_ARG roman_Π start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT italic_m start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT / ( 2 ∑ start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT italic_m start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT ) end_ARG, and the relative kinetic energy operator reads as

T^=i=1322μi2xi2=22μ4bi=132yi2.^𝑇superscriptsubscript𝑖13superscriptPlanck-constant-over-2-pi22subscript𝜇𝑖superscript2superscriptsubscript𝑥𝑖2superscriptPlanck-constant-over-2-pi22subscript𝜇4bsuperscriptsubscript𝑖13superscript2superscriptsubscript𝑦𝑖2\hat{T}=-\sum_{i=1}^{3}\frac{\hbar^{2}}{2\mu_{i}}\frac{\partial^{2}}{\partial x% _{i}^{2}}=-\frac{\hbar^{2}}{2\mu_{{\rm 4b}}}\sum_{i=1}^{3}\frac{\partial^{2}}{% \partial y_{i}^{2}}.over^ start_ARG italic_T end_ARG = - ∑ start_POSTSUBSCRIPT italic_i = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT divide start_ARG roman_ℏ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 italic_μ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT end_ARG divide start_ARG ∂ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG ∂ italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG = - divide start_ARG roman_ℏ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 italic_μ start_POSTSUBSCRIPT 4 roman_b end_POSTSUBSCRIPT end_ARG ∑ start_POSTSUBSCRIPT italic_i = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT divide start_ARG ∂ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG ∂ italic_y start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG . (2)

Hereafter, we focus on the equal-mass case: of mj=msubscript𝑚𝑗𝑚m_{j}=mitalic_m start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT = italic_m for all j𝑗jitalic_j. Then μ4b=m/2subscript𝜇4b𝑚2\mu_{{\rm 4b}}=m/2italic_μ start_POSTSUBSCRIPT 4 roman_b end_POSTSUBSCRIPT = italic_m / 2, μ1=m/2subscript𝜇1𝑚2\mu_{1}=m/2italic_μ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = italic_m / 2, μ2=m/2subscript𝜇2𝑚2\mu_{2}=m/2italic_μ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = italic_m / 2, and μ3=msubscript𝜇3𝑚\mu_{3}=mitalic_μ start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT = italic_m. We assume particles 1 and 2 are identical fermions, as are particles 3 and 4, but the two pairs are distinguishable from each other. The pairwise distances rijsubscript𝑟𝑖𝑗r_{ij}italic_r start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT between the particles i𝑖iitalic_i and j𝑗jitalic_j become:

r12=y1,r13=y1/2y2/2+y3/2,r14=y1/2+y2/2+y32,r23=y1/2y2/2+y3/2,r24=y1/2+y2/2+y3/2,r34=y2.formulae-sequencesubscript𝑟12subscript𝑦1formulae-sequencesubscript𝑟13subscript𝑦12subscript𝑦22subscript𝑦32formulae-sequencesubscript𝑟14subscript𝑦12subscript𝑦22subscript𝑦32formulae-sequencesubscript𝑟23subscript𝑦12subscript𝑦22subscript𝑦32formulae-sequencesubscript𝑟24subscript𝑦12subscript𝑦22subscript𝑦32subscript𝑟34subscript𝑦2\begin{gathered}r_{12}=y_{1},\\ r_{13}=y_{1}/2-y_{2}/2+y_{3}/\sqrt{2},\\ r_{14}=y_{1}/2+y_{2}/2+y_{3}\sqrt{2},\\ r_{23}=-y_{1}/2-y_{2}/2+y_{3}/\sqrt{2},\\ r_{24}=-y_{1}/2+y_{2}/2+y_{3}/\sqrt{2},\\ r_{34}=y_{2}.\end{gathered}start_ROW start_CELL italic_r start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = italic_y start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , end_CELL end_ROW start_ROW start_CELL italic_r start_POSTSUBSCRIPT 13 end_POSTSUBSCRIPT = italic_y start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT / 2 - italic_y start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT / 2 + italic_y start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT / square-root start_ARG 2 end_ARG , end_CELL end_ROW start_ROW start_CELL italic_r start_POSTSUBSCRIPT 14 end_POSTSUBSCRIPT = italic_y start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT / 2 + italic_y start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT / 2 + italic_y start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT square-root start_ARG 2 end_ARG , end_CELL end_ROW start_ROW start_CELL italic_r start_POSTSUBSCRIPT 23 end_POSTSUBSCRIPT = - italic_y start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT / 2 - italic_y start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT / 2 + italic_y start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT / square-root start_ARG 2 end_ARG , end_CELL end_ROW start_ROW start_CELL italic_r start_POSTSUBSCRIPT 24 end_POSTSUBSCRIPT = - italic_y start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT / 2 + italic_y start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT / 2 + italic_y start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT / square-root start_ARG 2 end_ARG , end_CELL end_ROW start_ROW start_CELL italic_r start_POSTSUBSCRIPT 34 end_POSTSUBSCRIPT = italic_y start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT . end_CELL end_ROW (3)

We introduce hyperspherical coordinates in analogy with three-dimensional (3D) spherical coordinate:

y1=Rsinϕsinθy2=Rcosϕsinθy3=Rcosθsubscript𝑦1𝑅italic-ϕ𝜃subscript𝑦2𝑅italic-ϕ𝜃subscript𝑦3𝑅𝜃\begin{gathered}y_{1}=R\sin\phi\sin\theta\\ y_{2}=R\cos\phi\sin\theta\\ y_{3}=R\cos\theta\end{gathered}start_ROW start_CELL italic_y start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = italic_R roman_sin italic_ϕ roman_sin italic_θ end_CELL end_ROW start_ROW start_CELL italic_y start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = italic_R roman_cos italic_ϕ roman_sin italic_θ end_CELL end_ROW start_ROW start_CELL italic_y start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT = italic_R roman_cos italic_θ end_CELL end_ROW (4)

with the hyperradius R[0,)𝑅0R\in[0,\infty)italic_R ∈ [ 0 , ∞ ), and the hyperangles Ω{θ,ϕ}Ω𝜃italic-ϕ\Omega\equiv\{\theta,\phi\}roman_Ω ≡ { italic_θ , italic_ϕ } where ϕ[0,2π)italic-ϕ02𝜋\phi\in[0,2\pi)italic_ϕ ∈ [ 0 , 2 italic_π ) and θ[0,π]𝜃0𝜋\theta\in[0,\pi]italic_θ ∈ [ 0 , italic_π ]. Figure 1 (b) illustrates the mapping between the geometric configuration of the four particles and the hyperangles. For clarity, we only show the region ϕ[0,π/2]italic-ϕ0𝜋2\phi\in[0,\pi/2]italic_ϕ ∈ [ 0 , italic_π / 2 ] and θ[0,π]𝜃0𝜋\theta\in[0,\pi]italic_θ ∈ [ 0 , italic_π ], as the full configuration space can be reconstructed using symmetry operations (see discussion below).

The relative kinetic energy operator becomes:

T^=22μ4b1R2R2R+L^(Ω)22μ4bR2,^𝑇superscriptPlanck-constant-over-2-pi22subscript𝜇4b1𝑅superscript2superscript𝑅2𝑅^𝐿superscriptΩ22subscript𝜇4bsuperscript𝑅2\hat{T}=-\frac{\hbar^{2}}{2\mu_{{\rm 4b}}}\frac{1}{R}\frac{\partial^{2}}{% \partial R^{2}}R+\frac{\hat{L}\left(\Omega\right)^{2}}{2\mu_{{\rm 4b}}R^{2}},over^ start_ARG italic_T end_ARG = - divide start_ARG roman_ℏ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 italic_μ start_POSTSUBSCRIPT 4 roman_b end_POSTSUBSCRIPT end_ARG divide start_ARG 1 end_ARG start_ARG italic_R end_ARG divide start_ARG ∂ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG ∂ italic_R start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG italic_R + divide start_ARG over^ start_ARG italic_L end_ARG ( roman_Ω ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 italic_μ start_POSTSUBSCRIPT 4 roman_b end_POSTSUBSCRIPT italic_R start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG , (5)

with the hyperangular momentum operator defined as:

L^(Ω)2=2[1sinθθsinθθ+1sin2θ2ϕ2].^𝐿superscriptΩ2superscriptPlanck-constant-over-2-pi2delimited-[]1𝜃𝜃𝜃𝜃1superscript2𝜃superscript2superscriptitalic-ϕ2\hat{L}\left(\Omega\right)^{2}=-\hbar^{2}\left[\frac{1}{\sin\theta}\frac{% \partial}{\partial\theta}\sin\theta\frac{\partial}{\partial\theta}+\frac{1}{% \sin^{2}\theta}\frac{\partial^{2}}{\partial\phi^{2}}\right].over^ start_ARG italic_L end_ARG ( roman_Ω ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT = - roman_ℏ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT [ divide start_ARG 1 end_ARG start_ARG roman_sin italic_θ end_ARG divide start_ARG ∂ end_ARG start_ARG ∂ italic_θ end_ARG roman_sin italic_θ divide start_ARG ∂ end_ARG start_ARG ∂ italic_θ end_ARG + divide start_ARG 1 end_ARG start_ARG roman_sin start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_θ end_ARG divide start_ARG ∂ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG ∂ italic_ϕ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ] . (6)

This operator has the same form as the angular momentum operator in three-dimensional space, which is expected since the number of relative degrees of freedom in our four-body problem is also three.

The four-body Schrödinger equation (after separating the center-of-mass degree of freedom) is given by

[22μ4b1R2R2R+L^(Ω)22μ4bR2+V(R,Ω)]Ψ(R,Ω)=EΨ(R,Ω).delimited-[]superscriptPlanck-constant-over-2-pi22subscript𝜇4b1𝑅superscript2superscript𝑅2𝑅^𝐿superscriptΩ22subscript𝜇4bsuperscript𝑅2𝑉𝑅ΩΨ𝑅Ω𝐸Ψ𝑅Ω\left[-\frac{\hbar^{2}}{2\mu_{{\rm 4b}}}\frac{1}{R}\frac{\partial^{2}}{% \partial R^{2}}R+\frac{\hat{L}\left(\Omega\right)^{2}}{2\mu_{{\rm 4b}}R^{2}}+V% \left(R,\Omega\right)\right]\Psi\left(R,\Omega\right)=E\Psi\left(R,\Omega% \right).[ - divide start_ARG roman_ℏ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 italic_μ start_POSTSUBSCRIPT 4 roman_b end_POSTSUBSCRIPT end_ARG divide start_ARG 1 end_ARG start_ARG italic_R end_ARG divide start_ARG ∂ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG ∂ italic_R start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG italic_R + divide start_ARG over^ start_ARG italic_L end_ARG ( roman_Ω ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 italic_μ start_POSTSUBSCRIPT 4 roman_b end_POSTSUBSCRIPT italic_R start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG + italic_V ( italic_R , roman_Ω ) ] roman_Ψ ( italic_R , roman_Ω ) = italic_E roman_Ψ ( italic_R , roman_Ω ) . (7)

We define a rescaled wave function ψ=RΨ𝜓𝑅Ψ\psi=R\Psiitalic_ψ = italic_R roman_Ψ, which satisfies

[22μ4b2R2+L^(Ω)22μ4bR2+V(R,Ω)]ψ(R,Ω)=Eψ(R,Ω),delimited-[]superscriptPlanck-constant-over-2-pi22subscript𝜇4bsuperscript2superscript𝑅2^𝐿superscriptΩ22subscript𝜇4bsuperscript𝑅2𝑉𝑅Ω𝜓𝑅Ω𝐸𝜓𝑅Ω\left[-\frac{\hbar^{2}}{2\mu_{{\rm 4b}}}\frac{\partial^{2}}{\partial R^{2}}+% \frac{\hat{L}\left(\Omega\right)^{2}}{2\mu_{{\rm 4b}}R^{2}}+V\left(R,\Omega% \right)\right]\psi\left(R,\Omega\right)=E\psi\left(R,\Omega\right),[ - divide start_ARG roman_ℏ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 italic_μ start_POSTSUBSCRIPT 4 roman_b end_POSTSUBSCRIPT end_ARG divide start_ARG ∂ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG ∂ italic_R start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG + divide start_ARG over^ start_ARG italic_L end_ARG ( roman_Ω ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 italic_μ start_POSTSUBSCRIPT 4 roman_b end_POSTSUBSCRIPT italic_R start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG + italic_V ( italic_R , roman_Ω ) ] italic_ψ ( italic_R , roman_Ω ) = italic_E italic_ψ ( italic_R , roman_Ω ) , (8)

where V(R,Ω)𝑉𝑅ΩV\left(R,\Omega\right)italic_V ( italic_R , roman_Ω ) denotes the interaction potential among particles.

Refer to caption
Figure 1: (a) Schematic illustration of the four-particle configuration. As the system is one-dimensional, only the horizontal displacement is physically meaningful; vertical offsets are used solely for visual clarity. (b) Mapping between the geometric configuration of particles and the hyperangular coordinates θ𝜃\thetaitalic_θ and ϕitalic-ϕ\phiitalic_ϕ. (c) The four-body potential V(R,Ω)𝑉𝑅ΩV(R,\Omega)italic_V ( italic_R , roman_Ω ) plotted as a function of angles at fixed hyperradius R=10𝑅10R=10italic_R = 10. The pairwise two-body interaction is modeled using the sinh\sinhroman_sinh-cosh\coshroman_cosh potential with parameters s=s=1.5𝑠superscript𝑠1.5s=s^{\prime}=1.5italic_s = italic_s start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT = 1.5.

2.2 The sinh\sinhroman_sinh-cosh\coshroman_cosh potential

For ultracold atomic systems, the total interaction potential can usually be modeled as a sum of pairwise interactions:

V(R,Ω)=i<jvij(rij),𝑉𝑅Ωsubscript𝑖𝑗subscript𝑣𝑖𝑗subscript𝑟𝑖𝑗V\left(R,\Omega\right)=\sum_{i<j}v_{ij}\left(r_{ij}\right),italic_V ( italic_R , roman_Ω ) = ∑ start_POSTSUBSCRIPT italic_i < italic_j end_POSTSUBSCRIPT italic_v start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT ( italic_r start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT ) , (9)

where the interparticle distances rijsubscript𝑟𝑖𝑗r_{ij}italic_r start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT can then be explicitly written in terms of hyperspherical coordinates:

r12=Rsinϕsinθ,r13=R[sinθsin(ϕπ/4)+cosθ]/2,r14=R[sinθsin(ϕ+π/4)+cosθ]/2,r23=R[sinθsin(ϕ+π/4)cosθ]/2,r24=R[sinθsin(ϕπ/4)cosθ]/2,r34=Rcosϕsinθ.formulae-sequencesubscript𝑟12𝑅italic-ϕ𝜃formulae-sequencesubscript𝑟13𝑅delimited-[]𝜃italic-ϕ𝜋4𝜃2formulae-sequencesubscript𝑟14𝑅delimited-[]𝜃italic-ϕ𝜋4𝜃2formulae-sequencesubscript𝑟23𝑅delimited-[]𝜃italic-ϕ𝜋4𝜃2formulae-sequencesubscript𝑟24𝑅delimited-[]𝜃italic-ϕ𝜋4𝜃2subscript𝑟34𝑅italic-ϕ𝜃\begin{gathered}r_{12}=R\sin\phi\sin\theta,\\ r_{13}=R\left[\sin\theta\sin\left(\phi-\pi/4\right)+\cos\theta\right]/\sqrt{2}% ,\\ r_{14}=R\left[\sin\theta\sin\left(\phi+\pi/4\right)+\cos\theta\right]/\sqrt{2}% ,\\ r_{23}=-R\left[\sin\theta\sin\left(\phi+\pi/4\right)-\cos\theta\right]/\sqrt{2% },\\ r_{24}=-R\left[\sin\theta\sin\left(\phi-\pi/4\right)-\cos\theta\right]/\sqrt{2% },\\ r_{34}=R\cos\phi\sin\theta.\end{gathered}start_ROW start_CELL italic_r start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = italic_R roman_sin italic_ϕ roman_sin italic_θ , end_CELL end_ROW start_ROW start_CELL italic_r start_POSTSUBSCRIPT 13 end_POSTSUBSCRIPT = italic_R [ roman_sin italic_θ roman_sin ( italic_ϕ - italic_π / 4 ) + roman_cos italic_θ ] / square-root start_ARG 2 end_ARG , end_CELL end_ROW start_ROW start_CELL italic_r start_POSTSUBSCRIPT 14 end_POSTSUBSCRIPT = italic_R [ roman_sin italic_θ roman_sin ( italic_ϕ + italic_π / 4 ) + roman_cos italic_θ ] / square-root start_ARG 2 end_ARG , end_CELL end_ROW start_ROW start_CELL italic_r start_POSTSUBSCRIPT 23 end_POSTSUBSCRIPT = - italic_R [ roman_sin italic_θ roman_sin ( italic_ϕ + italic_π / 4 ) - roman_cos italic_θ ] / square-root start_ARG 2 end_ARG , end_CELL end_ROW start_ROW start_CELL italic_r start_POSTSUBSCRIPT 24 end_POSTSUBSCRIPT = - italic_R [ roman_sin italic_θ roman_sin ( italic_ϕ - italic_π / 4 ) - roman_cos italic_θ ] / square-root start_ARG 2 end_ARG , end_CELL end_ROW start_ROW start_CELL italic_r start_POSTSUBSCRIPT 34 end_POSTSUBSCRIPT = italic_R roman_cos italic_ϕ roman_sin italic_θ . end_CELL end_ROW (10)

The system we now wish to focus consists of two kinds of fermions with the same mass m𝑚mitalic_m, distinguished only by a label σ=±1𝜎plus-or-minus1\sigma=\pm 1italic_σ = ± 1. The pair-wise potential considered here are called the sinh\sinhroman_sinh-cosh\coshroman_cosh potential [44], which follows the form of

vσiσj(rij)=s(s+1)sinh2(rij/rc)(1+σiσj2)s(s+1)cosh2(rij/rc)(1σiσj2),subscript𝑣subscript𝜎𝑖subscript𝜎𝑗subscript𝑟𝑖𝑗superscript𝑠superscript𝑠1superscript2subscript𝑟𝑖𝑗subscript𝑟𝑐1subscript𝜎𝑖subscript𝜎𝑗2𝑠𝑠1superscript2subscript𝑟𝑖𝑗subscript𝑟𝑐1subscript𝜎𝑖subscript𝜎𝑗2v_{\sigma_{i}\sigma_{j}}(r_{ij})=\frac{s^{\prime}\left(s^{\prime}+1\right)}{% \sinh^{2}\left(r_{ij}/r_{c}\right)}\left(\frac{1+\sigma_{i}\sigma_{j}}{2}% \right)-\frac{s\left(s+1\right)}{\cosh^{2}\left(r_{ij}/r_{c}\right)}\left(% \frac{1-\sigma_{i}\sigma_{j}}{2}\right),italic_v start_POSTSUBSCRIPT italic_σ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_σ start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT end_POSTSUBSCRIPT ( italic_r start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT ) = divide start_ARG italic_s start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ( italic_s start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT + 1 ) end_ARG start_ARG roman_sinh start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_r start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT / italic_r start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT ) end_ARG ( divide start_ARG 1 + italic_σ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_σ start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG ) - divide start_ARG italic_s ( italic_s + 1 ) end_ARG start_ARG roman_cosh start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_r start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT / italic_r start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT ) end_ARG ( divide start_ARG 1 - italic_σ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_σ start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG ) , (11)

where s,s>0superscript𝑠𝑠0s^{\prime},s>0italic_s start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT , italic_s > 0. This form ensures that like particles (same σ𝜎\sigmaitalic_σ) interact repulsively, while unlike particles experience an attractive interaction. In this sense, the system is analogous to an electron-hole model, where σ𝜎\sigmaitalic_σ can be interpreted as an effective charge. At large distances, both types of interactions decay exponentially with a characteristic length scale rcsubscript𝑟𝑐r_{c}italic_r start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT, which we hereafter set to unity. The corresponding natural energy scale is then Ec=2/mrc2subscript𝐸𝑐superscriptPlanck-constant-over-2-pi2𝑚superscriptsubscript𝑟𝑐2E_{c}=\hbar^{2}/mr_{c}^{2}italic_E start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT = roman_ℏ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT / italic_m italic_r start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT. Throughout this work, we adopt rcsubscript𝑟𝑐r_{c}italic_r start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT and Ecsubscript𝐸𝑐E_{c}italic_E start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT as natural units to interpret the length and energy scales in the figures.

Figure 1 (c) shows the four-body potential at hyperradius R=10𝑅10R=10italic_R = 10 as a function of the hyperangles. The plot reveals a clear fourfold symmetry, with identical potential landscapes across four regions. Notably, a deep minimum appears near (θ,ϕ)=(π/2,π/4)𝜃italic-ϕ𝜋2𝜋4\left(\theta,\phi\right)=\left(\pi/2,\pi/4\right)( italic_θ , italic_ϕ ) = ( italic_π / 2 , italic_π / 4 ), indicated by the dark blue region, while a lighter blue figure-eight-shaped stripe and a broader light green plateau surround it. As we will see, these regions correspond to different physical configurations: the deep minimum supports dimer-dimer scattering states, the lighter stripe indicates dimer–atom–atom configurations, and the outer plateau corresponds to four free atoms.

One interesting property of the sinh\sinhroman_sinh-cosh\coshroman_cosh potential is that the two-body scattering properties are known analytically [44]. For the repulsive sinh\sinhroman_sinh-potential between like particles, with the relative momentum defined as k=k1k2𝑘subscript𝑘1subscript𝑘2k=k_{1}-k_{2}italic_k = italic_k start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - italic_k start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT, the scattering amplitude is given as

S(k)=Γ(1+ik/2)Γ(s+1ik/2)Γ(1ik/2)Γ(s+1+ik/2),𝑆𝑘Γ1𝑖𝑘2Γ𝑠1𝑖𝑘2Γ1𝑖𝑘2Γ𝑠1𝑖𝑘2S\left(k\right)=-\frac{\Gamma\left(1+ik/2\right)\Gamma\left(s+1-ik/2\right)}{% \Gamma\left(1-ik/2\right)\Gamma\left(s+1+ik/2\right)},italic_S ( italic_k ) = - divide start_ARG roman_Γ ( 1 + italic_i italic_k / 2 ) roman_Γ ( italic_s + 1 - italic_i italic_k / 2 ) end_ARG start_ARG roman_Γ ( 1 - italic_i italic_k / 2 ) roman_Γ ( italic_s + 1 + italic_i italic_k / 2 ) end_ARG , (12)

where Γ()Γ\Gamma\left(\cdot\right)roman_Γ ( ⋅ ) is the gamma function. For unlike particles, the cosh\coshroman_cosh-potential are attractive, whose reflection and transmission amplitudes are given as

R(k)=S(k)r(k),T(k)=S(k)t(k),formulae-sequence𝑅𝑘𝑆𝑘𝑟𝑘𝑇𝑘𝑆𝑘𝑡𝑘R\left(k\right)=S\left(k\right)r\left(k\right),\ T\left(k\right)=S\left(k% \right)t\left(k\right),italic_R ( italic_k ) = italic_S ( italic_k ) italic_r ( italic_k ) , italic_T ( italic_k ) = italic_S ( italic_k ) italic_t ( italic_k ) , (13)

where

r(k)=sinπ(s+1)sinπ(s+1+ik/2),t(k)=sinπik/2sinπ(s+1+ik/2).formulae-sequence𝑟𝑘𝜋𝑠1𝜋𝑠1𝑖𝑘2𝑡𝑘𝜋𝑖𝑘2𝜋𝑠1𝑖𝑘2\begin{gathered}r\left(k\right)=\frac{\sin\pi\left(s+1\right)}{\sin\pi\left(s+% 1+ik/2\right)},\\ t\left(k\right)=\frac{\sin\pi ik/2}{\sin\pi\left(s+1+ik/2\right)}.\end{gathered}start_ROW start_CELL italic_r ( italic_k ) = divide start_ARG roman_sin italic_π ( italic_s + 1 ) end_ARG start_ARG roman_sin italic_π ( italic_s + 1 + italic_i italic_k / 2 ) end_ARG , end_CELL end_ROW start_ROW start_CELL italic_t ( italic_k ) = divide start_ARG roman_sin italic_π italic_i italic_k / 2 end_ARG start_ARG roman_sin italic_π ( italic_s + 1 + italic_i italic_k / 2 ) end_ARG . end_CELL end_ROW (14)

Dimer bound states appear as poles of the reflection and transmission amplitudes, with binding energy Eτ=(s+1τ)2subscript𝐸𝜏superscript𝑠1𝜏2E_{\tau}=-\left(s+1-\tau\right)^{2}italic_E start_POSTSUBSCRIPT italic_τ end_POSTSUBSCRIPT = - ( italic_s + 1 - italic_τ ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT, τ=1,2,,s+1𝜏12𝑠1\tau=1,2,\cdots,\left\lfloor s+1\right\rflooritalic_τ = 1 , 2 , ⋯ , ⌊ italic_s + 1 ⌋, where x𝑥\left\lfloor x\right\rfloor⌊ italic_x ⌋ denote the largest integer less than x𝑥xitalic_x. The dimer can thus be described by two particles with momentums k±iκτplus-or-minus𝑘𝑖subscript𝜅𝜏k\pm i\kappa_{\tau}italic_k ± italic_i italic_κ start_POSTSUBSCRIPT italic_τ end_POSTSUBSCRIPT, where κτ=s+1τsubscript𝜅𝜏𝑠1𝜏\kappa_{\tau}=s+1-\tauitalic_κ start_POSTSUBSCRIPT italic_τ end_POSTSUBSCRIPT = italic_s + 1 - italic_τ.

For the case of s=s𝑠superscript𝑠s=s^{\prime}italic_s = italic_s start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT, the system become integrable for arbitrary number of particles. No other bound states exists other than the dimer bound states mentioned above. The atom-dimer scattering amplitude between a particle with momentum k1subscript𝑘1k_{1}italic_k start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT and a pair of bounded particles with momenta k2±iκτplus-or-minussubscript𝑘2𝑖subscript𝜅𝜏k_{2}\pm i\kappa_{\tau}italic_k start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ± italic_i italic_κ start_POSTSUBSCRIPT italic_τ end_POSTSUBSCRIPT is given by

exp[iθ0τ(k)]=S(k+iκτ)S(kiκτ)t(k+iκτ)𝑖subscript𝜃0𝜏𝑘𝑆𝑘𝑖subscript𝜅𝜏𝑆𝑘𝑖subscript𝜅𝜏𝑡𝑘𝑖subscript𝜅𝜏\exp\left[-i\theta_{0\tau}(k)\right]=S\left(k+i\kappa_{\tau}\right)S\left(k-i% \kappa_{\tau}\right)t\left(k+i\kappa_{\tau}\right)roman_exp [ - italic_i italic_θ start_POSTSUBSCRIPT 0 italic_τ end_POSTSUBSCRIPT ( italic_k ) ] = italic_S ( italic_k + italic_i italic_κ start_POSTSUBSCRIPT italic_τ end_POSTSUBSCRIPT ) italic_S ( italic_k - italic_i italic_κ start_POSTSUBSCRIPT italic_τ end_POSTSUBSCRIPT ) italic_t ( italic_k + italic_i italic_κ start_POSTSUBSCRIPT italic_τ end_POSTSUBSCRIPT ) (15)

with k=k2k1𝑘subscript𝑘2subscript𝑘1k=k_{2}-k_{1}italic_k = italic_k start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT - italic_k start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT. The dimer-dimer scattering amplitude Sττ(k)=exp[iθττ(k)]subscript𝑆𝜏superscript𝜏𝑘𝑖subscript𝜃𝜏superscript𝜏𝑘S_{\tau\tau^{\prime}}(k)=\exp\left[-i\theta_{\tau\tau^{\prime}}\left(k\right)\right]italic_S start_POSTSUBSCRIPT italic_τ italic_τ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT ( italic_k ) = roman_exp [ - italic_i italic_θ start_POSTSUBSCRIPT italic_τ italic_τ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT ( italic_k ) ] between bound atom pairs with k1±iκτplus-or-minussubscript𝑘1𝑖subscript𝜅𝜏k_{1}\pm i\kappa_{\tau}italic_k start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ± italic_i italic_κ start_POSTSUBSCRIPT italic_τ end_POSTSUBSCRIPT and k2±iκτplus-or-minussubscript𝑘2𝑖subscript𝜅superscript𝜏k_{2}\pm i\kappa_{\tau^{\prime}}italic_k start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ± italic_i italic_κ start_POSTSUBSCRIPT italic_τ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT are determined by the phase shift

θττ(k)=θ0τ(kiκτ)+θ0τ(k+iκτ),subscript𝜃𝜏superscript𝜏𝑘subscript𝜃0superscript𝜏𝑘𝑖subscript𝜅𝜏subscript𝜃0superscript𝜏𝑘𝑖subscript𝜅𝜏\theta_{\tau\tau^{\prime}}(k)=\theta_{0\tau^{\prime}}(k-i\kappa_{\tau})+\theta% _{0\tau^{\prime}}(k+i\kappa_{\tau}),italic_θ start_POSTSUBSCRIPT italic_τ italic_τ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT ( italic_k ) = italic_θ start_POSTSUBSCRIPT 0 italic_τ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT ( italic_k - italic_i italic_κ start_POSTSUBSCRIPT italic_τ end_POSTSUBSCRIPT ) + italic_θ start_POSTSUBSCRIPT 0 italic_τ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT ( italic_k + italic_i italic_κ start_POSTSUBSCRIPT italic_τ end_POSTSUBSCRIPT ) , (16)

with k=k2k1𝑘subscript𝑘2subscript𝑘1k=k_{2}-k_{1}italic_k = italic_k start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT - italic_k start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT. For scattering between the ground state dimers τ=τ=1𝜏superscript𝜏1\tau=\tau^{\prime}=1italic_τ = italic_τ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT = 1, we have κ1=ssubscript𝜅1𝑠\kappa_{1}=sitalic_κ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = italic_s, and E1=s2subscript𝐸1superscript𝑠2E_{1}=-s^{2}italic_E start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = - italic_s start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT. For s1𝑠1s\geq 1italic_s ≥ 1, the ground state dimer-dimer scattering amplitude has a simplified analytical form

𝒮(k)=(sik/2)(s+ik/2)Γ(1+ik/2)Γ(1ik/2)Γ(2s+1ik/2)Γ(2s+1+ik/2).𝒮𝑘𝑠𝑖𝑘2𝑠𝑖𝑘2Γ1𝑖𝑘2Γ1𝑖𝑘2Γ2𝑠1𝑖𝑘2Γ2𝑠1𝑖𝑘2\mathcal{S}(k)=-\frac{\left(s-ik/2\right)}{\left(s+ik/2\right)}\frac{\Gamma% \left(1+ik/2\right)}{\Gamma\left(1-ik/2\right)}\frac{\Gamma\left(2s+1-ik/2% \right)}{\Gamma\left(2s+1+ik/2\right)}.caligraphic_S ( italic_k ) = - divide start_ARG ( italic_s - italic_i italic_k / 2 ) end_ARG start_ARG ( italic_s + italic_i italic_k / 2 ) end_ARG divide start_ARG roman_Γ ( 1 + italic_i italic_k / 2 ) end_ARG start_ARG roman_Γ ( 1 - italic_i italic_k / 2 ) end_ARG divide start_ARG roman_Γ ( 2 italic_s + 1 - italic_i italic_k / 2 ) end_ARG start_ARG roman_Γ ( 2 italic_s + 1 + italic_i italic_k / 2 ) end_ARG . (17)

Since this is the lowest scattering channel, this scattering amplitude is equivalent to the S𝑆Sitalic_S-matrix. The closely related K𝐾Kitalic_K-matrix is given by

𝒦(k)=i1𝒮(k)1+𝒮(k),𝒦𝑘𝑖1𝒮𝑘1𝒮𝑘\mathcal{K}(k)=i\frac{1-\mathcal{S}\left(k\right)}{1+\mathcal{S}\left(k\right)},caligraphic_K ( italic_k ) = italic_i divide start_ARG 1 - caligraphic_S ( italic_k ) end_ARG start_ARG 1 + caligraphic_S ( italic_k ) end_ARG , (18)

and the energy-dependent dimer-dimer scattering length is given by

aD(k)=1k𝒦(k).subscript𝑎𝐷𝑘1𝑘𝒦𝑘a_{D}(k)=\frac{1}{k\mathcal{K}(k)}.italic_a start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT ( italic_k ) = divide start_ARG 1 end_ARG start_ARG italic_k caligraphic_K ( italic_k ) end_ARG . (19)

The dimer-dimer scattering length at zero scattering energy is thus given by

aDlimk0aD(k)=γ+s1+ψ(2s+1)2,subscript𝑎𝐷subscript𝑘0subscript𝑎𝐷𝑘𝛾superscript𝑠1𝜓2𝑠12a_{D}\equiv\lim_{k\rightarrow 0}a_{D}\left(k\right)=\frac{\gamma+s^{-1}+\psi% \left(2s+1\right)}{2},italic_a start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT ≡ roman_lim start_POSTSUBSCRIPT italic_k → 0 end_POSTSUBSCRIPT italic_a start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT ( italic_k ) = divide start_ARG italic_γ + italic_s start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT + italic_ψ ( 2 italic_s + 1 ) end_ARG start_ARG 2 end_ARG , (20)

where ψ()𝜓\psi\left(\cdot\right)italic_ψ ( ⋅ ) is the digamma function and γ=0.577216𝛾0.577216\gamma=0.577216italic_γ = 0.577216 is the Euler-Mascheroni constant. We will use these analytical results to benchmark our numerical calculations.

3 Results

3.1 Hyperspherical potential

To efficiently solve Eq. 8, we adopt the AHR. This involves solving the hyperangular eigenvalue problem:

[L^(Ω)22μ4bR2+V(R,Ω)]Φν(R;Ω)=Uν(R)Φν(R;Ω),delimited-[]^𝐿superscriptΩ22subscript𝜇4bsuperscript𝑅2𝑉𝑅ΩsubscriptΦ𝜈𝑅Ωsubscript𝑈𝜈𝑅subscriptΦ𝜈𝑅Ω\left[\frac{\hat{L}\left(\Omega\right)^{2}}{2\mu_{{\rm 4b}}R^{2}}+V\left(R,% \Omega\right)\right]\Phi_{\nu}(R;\Omega)=U_{\nu}\left(R\right)\Phi_{\nu}\left(% R;\Omega\right),[ divide start_ARG over^ start_ARG italic_L end_ARG ( roman_Ω ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 italic_μ start_POSTSUBSCRIPT 4 roman_b end_POSTSUBSCRIPT italic_R start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG + italic_V ( italic_R , roman_Ω ) ] roman_Φ start_POSTSUBSCRIPT italic_ν end_POSTSUBSCRIPT ( italic_R ; roman_Ω ) = italic_U start_POSTSUBSCRIPT italic_ν end_POSTSUBSCRIPT ( italic_R ) roman_Φ start_POSTSUBSCRIPT italic_ν end_POSTSUBSCRIPT ( italic_R ; roman_Ω ) , (21)

where Uν(R)subscript𝑈𝜈𝑅U_{\nu}\left(R\right)italic_U start_POSTSUBSCRIPT italic_ν end_POSTSUBSCRIPT ( italic_R ) are the hyperspherical potentials, and Φν(Ω;R)subscriptΦ𝜈Ω𝑅\Phi_{\nu}(\Omega;R)roman_Φ start_POSTSUBSCRIPT italic_ν end_POSTSUBSCRIPT ( roman_Ω ; italic_R ) are the corresponding channel functions, which depend parametrically on hyperradius R𝑅Ritalic_R. These functions form a truncated basis for expanding the total wave function:

ψνE(R,Ω)=ν=1NcFνν(R)Φν(Ω;R),superscriptsubscript𝜓superscript𝜈𝐸𝑅Ωsuperscriptsubscript𝜈1subscript𝑁𝑐subscript𝐹𝜈superscript𝜈𝑅subscriptΦ𝜈Ω𝑅\psi_{\nu^{\prime}}^{E}(R,\Omega)=\sum_{\nu=1}^{N_{c}}F_{\nu\nu^{\prime}}(R)% \Phi_{\nu}(\Omega;R),italic_ψ start_POSTSUBSCRIPT italic_ν start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_E end_POSTSUPERSCRIPT ( italic_R , roman_Ω ) = ∑ start_POSTSUBSCRIPT italic_ν = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_N start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT end_POSTSUPERSCRIPT italic_F start_POSTSUBSCRIPT italic_ν italic_ν start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT ( italic_R ) roman_Φ start_POSTSUBSCRIPT italic_ν end_POSTSUBSCRIPT ( roman_Ω ; italic_R ) , (22)

where νsuperscript𝜈\nu^{\prime}italic_ν start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT indexes the independent solutions and Ncsubscript𝑁𝑐N_{c}italic_N start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT is the number of channels retained in the expansion.

Before proceeding further, we note that the calculation can often be simplified by imposing permutation symmetry for identical particles. In what follows, we focus on the case of four fermions, where particles 1 and 2 are assigned σ=1𝜎1\sigma=1italic_σ = 1, and particles 3 and 4 are assigned σ=1𝜎1\sigma=-1italic_σ = - 1. The total wave function must be antisymmetrized using the operator

𝒜=1𝒫12𝒫34+𝒫12𝒫34,𝒜1subscript𝒫12subscript𝒫34subscript𝒫12subscript𝒫34\mathcal{A}=1-\mathcal{P}_{12}-\mathcal{P}_{34}+\mathcal{P}_{12}\mathcal{P}_{3% 4},caligraphic_A = 1 - caligraphic_P start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT - caligraphic_P start_POSTSUBSCRIPT 34 end_POSTSUBSCRIPT + caligraphic_P start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT caligraphic_P start_POSTSUBSCRIPT 34 end_POSTSUBSCRIPT , (23)

where Pijsubscript𝑃𝑖𝑗P_{ij}italic_P start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT denotes the exchange operator for particles i𝑖iitalic_i and j𝑗jitalic_j. From the expressions of the interparticle distances rijsubscript𝑟𝑖𝑗r_{ij}italic_r start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT in Eq. (10), we find that P12subscript𝑃12P_{12}italic_P start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT pas ϕ2πϕitalic-ϕ2𝜋italic-ϕ\phi\rightarrow 2\pi-\phiitalic_ϕ → 2 italic_π - italic_ϕ, while P34subscript𝑃34P_{34}italic_P start_POSTSUBSCRIPT 34 end_POSTSUBSCRIPT maps ϕπϕitalic-ϕ𝜋italic-ϕ\phi\rightarrow\pi-\phiitalic_ϕ → italic_π - italic_ϕ, revealing the fourfold symmetry observed in the potential landscape shown in Fig. 1 (c). Since the hyperangular momentum operator takes the same form as the angular momentum operator for a single particle in three-dimensional space, the wave function can be formally expanded using spherical harmonics: um(θ,ϕ)Pm(cosθ)eimϕsimilar-tosubscript𝑢subscript𝑚𝜃italic-ϕsuperscriptsubscript𝑃subscript𝑚𝜃superscript𝑒𝑖subscript𝑚italic-ϕu_{\ell m_{\ell}}\left(\theta,\phi\right)\sim P_{\ell}^{m_{\ell}}\left(\cos% \theta\right)e^{im_{\ell}\phi}italic_u start_POSTSUBSCRIPT roman_ℓ italic_m start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT end_POSTSUBSCRIPT ( italic_θ , italic_ϕ ) ∼ italic_P start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_m start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT end_POSTSUPERSCRIPT ( roman_cos italic_θ ) italic_e start_POSTSUPERSCRIPT italic_i italic_m start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT italic_ϕ end_POSTSUPERSCRIPT. Applying the antisymmetrization operator 𝒜𝒜\mathcal{A}caligraphic_A to these functions yields antisymmetrized basis functions of the form uk𝒜(θ,ϕ)P2k(cosθ)sin(2kϕ)similar-tosuperscriptsubscript𝑢subscript𝑘𝒜𝜃italic-ϕsuperscriptsubscript𝑃2subscript𝑘𝜃2subscript𝑘italic-ϕu_{\ell k_{\ell}}^{\mathcal{A}}\left(\theta,\phi\right)\sim P_{\ell}^{2k_{\ell% }}\left(\cos\theta\right)\sin\left(2k_{\ell}\phi\right)italic_u start_POSTSUBSCRIPT roman_ℓ italic_k start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT end_POSTSUBSCRIPT start_POSTSUPERSCRIPT caligraphic_A end_POSTSUPERSCRIPT ( italic_θ , italic_ϕ ) ∼ italic_P start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 italic_k start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT end_POSTSUPERSCRIPT ( roman_cos italic_θ ) roman_sin ( 2 italic_k start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT italic_ϕ ), where ksubscript𝑘k_{\ell}italic_k start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT is a positive integers.

The pair-wise potential vσiσj(r)subscript𝑣subscript𝜎𝑖subscript𝜎𝑗𝑟v_{\sigma_{i}\sigma_{j}}(r)italic_v start_POSTSUBSCRIPT italic_σ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_σ start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT end_POSTSUBSCRIPT ( italic_r ) also exhibits certain symmetries: it is even in r𝑟ritalic_r, i.e., vσiσj(r)=vσiσj(r)subscript𝑣subscript𝜎𝑖subscript𝜎𝑗𝑟subscript𝑣subscript𝜎𝑖subscript𝜎𝑗𝑟v_{\sigma_{i}\sigma_{j}}(r)=v_{\sigma_{i}\sigma_{j}}(-r)italic_v start_POSTSUBSCRIPT italic_σ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_σ start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT end_POSTSUBSCRIPT ( italic_r ) = italic_v start_POSTSUBSCRIPT italic_σ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_σ start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT end_POSTSUBSCRIPT ( - italic_r ), and identical for like-particle pairs, i.e., v11(r)=v11(r)subscript𝑣11𝑟subscript𝑣11𝑟v_{11}\left(r\right)=v_{-1-1}\left(r\right)italic_v start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT ( italic_r ) = italic_v start_POSTSUBSCRIPT - 1 - 1 end_POSTSUBSCRIPT ( italic_r ). As a result, the full Hamiltonian is invariant under two commuting operations 𝒫θ:θπθ:subscript𝒫𝜃𝜃𝜋𝜃\mathcal{P}_{\theta}:\theta\rightarrow\pi-\thetacaligraphic_P start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT : italic_θ → italic_π - italic_θ and 𝒫ϕ:ϕπ/2ϕ:subscript𝒫italic-ϕitalic-ϕ𝜋2italic-ϕ\mathcal{P}_{\phi}:\phi\rightarrow\pi/2-\phicaligraphic_P start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT : italic_ϕ → italic_π / 2 - italic_ϕ. These symmetries allow us to define two good quantum numbers: pθ=()subscript𝑝𝜃superscriptp_{\theta}=\left(-\right)^{\ell}italic_p start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT = ( - ) start_POSTSUPERSCRIPT roman_ℓ end_POSTSUPERSCRIPT, pϕ=()k+1subscript𝑝italic-ϕsuperscriptsubscript𝑘1p_{\phi}=\left(-\right)^{k_{\ell}+1}italic_p start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT = ( - ) start_POSTSUPERSCRIPT italic_k start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT + 1 end_POSTSUPERSCRIPT. Consequently, the calculation can be restricted to the reduced domain θ[0,π/2]𝜃0𝜋2\theta\in[0,\pi/2]italic_θ ∈ [ 0 , italic_π / 2 ] and ϕ[0,π/4]italic-ϕ0𝜋4\phi\in[0,\pi/4]italic_ϕ ∈ [ 0 , italic_π / 4 ], with numerical basis ui(θ,ϕ)subscript𝑢𝑖𝜃italic-ϕu_{i}\left(\theta,\phi\right)italic_u start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( italic_θ , italic_ϕ ) subject to the appropreate boundary conditions. We find that the ground-state dimer-dimer channel lies in the symmetry sector with pθ=+1subscript𝑝𝜃1p_{\theta}=+1italic_p start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT = + 1 and pϕ=+1subscript𝑝italic-ϕ1p_{\phi}=+1italic_p start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT = + 1, corresponding to even \ellroman_ℓ and odd ksubscript𝑘k_{\ell}italic_k start_POSTSUBSCRIPT roman_ℓ end_POSTSUBSCRIPT, which implies ui|=ϕ=00u_{i}\left|{}_{\phi=0}\right.=0italic_u start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT | start_FLOATSUBSCRIPT italic_ϕ = 0 end_FLOATSUBSCRIPT = 0, ϕui|=ϕ=π/40\partial_{\phi}u_{i}\left|{}_{\phi=\pi/4}\right.=0∂ start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT italic_u start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT | start_FLOATSUBSCRIPT italic_ϕ = italic_π / 4 end_FLOATSUBSCRIPT = 0, θui|=θ=π/20\partial_{\theta}u_{i}\left|{}_{\theta=\pi/2}\right.=0∂ start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT italic_u start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT | start_FLOATSUBSCRIPT italic_θ = italic_π / 2 end_FLOATSUBSCRIPT = 0 and ui|θ=0u_{i}\left|{}_{\theta=0}\right.italic_u start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT | start_FLOATSUBSCRIPT italic_θ = 0 end_FLOATSUBSCRIPT remains finite. (Boundary conditions for other symmetry sector are summarized in Table 1.) We employ B-spline basis functions, which are piecewise polynomials with local support. This choice facilitates the implementation of boundary conditions and leads to a sparse Hamiltonian matrix that can be efficiently diagonalized using eigensolvers such as ARPACK. In practice, we typically use approximately 100100100100 B-splines for the θ𝜃\thetaitalic_θ coordinate and 50 for ϕitalic-ϕ\phiitalic_ϕ, resulting in a sparse matrix of dimension around 5,00050005,0005 , 000.

Table 1: Boundary conditions for the numerical basis functions in different symmetry sectors, characterized by pθsubscript𝑝𝜃p_{\theta}italic_p start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT and pϕsubscript𝑝italic-ϕp_{\phi}italic_p start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT.
\toprulepθsubscript𝑝𝜃p_{\theta}italic_p start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT pϕsubscript𝑝italic-ϕp_{\phi}italic_p start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT θ=0𝜃0\theta=0italic_θ = 0 θ=π/2𝜃𝜋2\theta=\pi/2italic_θ = italic_π / 2 ϕ=π/4italic-ϕ𝜋4\phi=\pi/4italic_ϕ = italic_π / 4 ϕ=π/4italic-ϕ𝜋4\phi=\pi/4italic_ϕ = italic_π / 4
\midrule+11+1+ 1 +11+1+ 1 Finite θui=0subscript𝜃subscript𝑢𝑖0\partial_{\theta}u_{i}=0∂ start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT italic_u start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = 0 ui=0subscript𝑢𝑖0u_{i}=0italic_u start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = 0 ϕui=0subscriptitalic-ϕsubscript𝑢𝑖0\partial_{\phi}u_{i}=0∂ start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT italic_u start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = 0
+11+1+ 1 11-1- 1 Finite θui=0subscript𝜃subscript𝑢𝑖0\partial_{\theta}u_{i}=0∂ start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT italic_u start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = 0 ui=0subscript𝑢𝑖0u_{i}=0italic_u start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = 0 ui=0subscript𝑢𝑖0u_{i}=0italic_u start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = 0
11-1- 1 +11+1+ 1 Finite ui=0subscript𝑢𝑖0u_{i}=0italic_u start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = 0 ui=0subscript𝑢𝑖0u_{i}=0italic_u start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = 0 ϕui=0subscriptitalic-ϕsubscript𝑢𝑖0\partial_{\phi}u_{i}=0∂ start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT italic_u start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = 0
11-1- 1 11-1- 1 Finite ui=0subscript𝑢𝑖0u_{i}=0italic_u start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = 0 ui=0subscript𝑢𝑖0u_{i}=0italic_u start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = 0 ui=0subscript𝑢𝑖0u_{i}=0italic_u start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = 0
\botrule

Figure 2 and 3 present the hyperspherical potentials for s=s=1.5𝑠superscript𝑠1.5s=s^{\prime}=1.5italic_s = italic_s start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT = 1.5 across different symmetry sectors, while Figure 4 shows the hyperspherical potentials for fixed s=1.5𝑠1.5s=1.5italic_s = 1.5 and varying ssuperscript𝑠s^{\prime}italic_s start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT. A common feature among these curves is their convergence toward specific asymptotic thresholds at large hyperradii. For s=1.5𝑠1.5s=1.5italic_s = 1.5, the attractive interaction between unlike particles supports two dimer bound states with binding energies E1=2.25subscript𝐸12.25E_{1}=-2.25italic_E start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = - 2.25 and E2=0.25subscript𝐸20.25E_{2}=-0.25italic_E start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = - 0.25. As a result, the dimer–dimer thresholds appear at E1+E1=4.5subscript𝐸1subscript𝐸14.5E_{1}+E_{1}=-4.5italic_E start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT + italic_E start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = - 4.5, E1+E2=2.5subscript𝐸1subscript𝐸22.5E_{1}+E_{2}=-2.5italic_E start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT + italic_E start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = - 2.5, and E2+E2=0.5subscript𝐸2subscript𝐸20.5E_{2}+E_{2}=-0.5italic_E start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT + italic_E start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = - 0.5, indicated by solid horizontal lines in the figures. The dimer–atom–atom thresholds, corresponding to single dimer binding energies E1subscript𝐸1E_{1}italic_E start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT and E2subscript𝐸2E_{2}italic_E start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT, are shown as dashed horizontal lines. The atom–atom–atom–atom threshold lies at zero energy and is indicated by dash-dotted lines. An exception occurs in the second-lowest channel of Fig. 4 (b), which does not approach any of the aforementioned thresholds. This behavior arises because the reduced repulsion between like particles at s=0.9superscript𝑠0.9s^{\prime}=0.9italic_s start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT = 0.9 permits the formation of trimer bound states, causing this channel to asymptotically approach a trimer–atom threshold. Additionally, many sharp avoided crossings are observed between potential curves that approach different thresholds, which may pose numerical challenges for accurate scattering calculations.

Figure 5 displays the channel functions |Φν(Ω;R)|2superscriptsubscriptΦ𝜈Ω𝑅2|\Phi_{\nu}(\Omega;R)|^{2}| roman_Φ start_POSTSUBSCRIPT italic_ν end_POSTSUBSCRIPT ( roman_Ω ; italic_R ) | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT for several representative configurations. In panel (a), we show the dimer–dimer channel function corresponding to the lowest hyperspherical potential at R=10𝑅10R=10italic_R = 10 in Fig. 2 (a). Panel (b) shows the dimer–atom–atom channel function from the second-lowest channel at the same R=10𝑅10R=10italic_R = 10 in Fig. 2 (a). Panel (c) presents the atom–atom–atom–atom channel function from the 41st channel at R=40𝑅40R=40italic_R = 40 in Fig. 2 (a). Finally, panel (d) shows the trimer–atom channel function corresponding to the lowest potential at R=10𝑅10R=10italic_R = 10 in Fig. 4 (b). These channel functions reveal the spatial structure of the associated scattering configurations. The dimer–dimer channel is primarily localized in the deep blue minimum of the four-body potential landscape [see Fig. 1 (c)], while the atom–atom–atom–atom configuration occupies the broader light green plateau. The dimer–atom–atom and trimer–atom configurations are mostly concentrated along the figure-eight-shaped light blue region, highlighting the distinct geometrical signatures of different few-body scattering processes.

Refer to caption
Figure 2: Hyperspherical potential curves for s=s=1.5𝑠superscript𝑠1.5s=s^{\prime}=1.5italic_s = italic_s start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT = 1.5 in symmetry sectors: (a) (pθ,pϕ)=(+1,+1)subscript𝑝𝜃subscript𝑝italic-ϕ11\left(p_{\theta},p_{\phi}\right)=\left(+1,+1\right)( italic_p start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT , italic_p start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT ) = ( + 1 , + 1 ) and (b) (pθ,pϕ)=(+1,1)subscript𝑝𝜃subscript𝑝italic-ϕ11\left(p_{\theta},p_{\phi}\right)=\left(+1,-1\right)( italic_p start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT , italic_p start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT ) = ( + 1 , - 1 ). Solid horizontal lines indicate dimer-dimer thresholds; dashed lines mark dimer–atom–atom thresholds; dash-dotted lines represent the atom-atom-atom-atom threshold at zero energy.
Refer to caption
Figure 3: Same as Fig. 2, but for symmetry sectors: (a) (pθ,pϕ)=(1,+1)subscript𝑝𝜃subscript𝑝italic-ϕ11\left(p_{\theta},p_{\phi}\right)=\left(-1,+1\right)( italic_p start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT , italic_p start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT ) = ( - 1 , + 1 ) and (b) (pθ,pϕ)=(1,1)subscript𝑝𝜃subscript𝑝italic-ϕ11\left(p_{\theta},p_{\phi}\right)=\left(-1,-1\right)( italic_p start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT , italic_p start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT ) = ( - 1 , - 1 ).
Refer to caption
Figure 4: Hyperspherical potential curves for (pθ,pϕ)=(+1,+1)subscript𝑝𝜃subscript𝑝italic-ϕ11\left(p_{\theta},p_{\phi}\right)=\left(+1,+1\right)( italic_p start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT , italic_p start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT ) = ( + 1 , + 1 ) and s=1.5𝑠1.5s=1.5italic_s = 1.5 with different ssuperscript𝑠s^{\prime}italic_s start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT: (a) s=1.9superscript𝑠1.9s^{\prime}=1.9italic_s start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT = 1.9 and (b) s=0.9superscript𝑠0.9s^{\prime}=0.9italic_s start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT = 0.9. In panel (b), the second-lowest channel approaches a trimer–atom threshold due to weakened repulsion between like particles, indicating trimer formation.
Refer to caption
Figure 5: Channel functions for (a) dimer-dimer channel (b) dimer-atom-atom channel (c) Four-atom channel and (d) trimer-atom channel. The color scale indicates the probability density of the wave function, shown in arbitrary units.

3.2 Scattering

With the adiabatic potential curves and diabatic corrections obtained from the channel functions, we solve the radial Schrödinger equation as a set of coupled differential equations. However, near avoided crossings, the diabatic corrections can exhibit rapid variations that lead to numerical instability. To address this challenge, we use the slow variable discretization (SVD) method with a radial discrete variable representation (DVR) basis πn(R)subscript𝜋𝑛𝑅\pi_{n}(R)italic_π start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT ( italic_R ) to expand Fνν(R)subscript𝐹𝜈superscript𝜈𝑅F_{\nu\nu^{\prime}}(R)italic_F start_POSTSUBSCRIPT italic_ν italic_ν start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT ( italic_R ), which gives

ψν(R,Ω)=ν,ncnν,νπn(R)Φν(Ω;Rn).subscript𝜓superscript𝜈𝑅Ωsubscript𝜈𝑛subscript𝑐𝑛𝜈superscript𝜈subscript𝜋𝑛𝑅subscriptΦ𝜈Ωsubscript𝑅𝑛\psi_{\nu^{\prime}}(R,\Omega)=\sum_{\nu,n}c_{n\nu,\nu^{\prime}}\pi_{n}(R)\Phi_% {\nu}\left(\Omega;R_{n}\right).italic_ψ start_POSTSUBSCRIPT italic_ν start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT ( italic_R , roman_Ω ) = ∑ start_POSTSUBSCRIPT italic_ν , italic_n end_POSTSUBSCRIPT italic_c start_POSTSUBSCRIPT italic_n italic_ν , italic_ν start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT italic_π start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT ( italic_R ) roman_Φ start_POSTSUBSCRIPT italic_ν end_POSTSUBSCRIPT ( roman_Ω ; italic_R start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT ) . (24)

The DVR approximation allows integrals of the form:

a1a2πi(R)H(R)πj(R)𝑑RH(Ri)δijsuperscriptsubscriptsubscript𝑎1subscript𝑎2subscript𝜋𝑖𝑅𝐻𝑅subscript𝜋𝑗𝑅differential-d𝑅𝐻subscript𝑅𝑖subscript𝛿𝑖𝑗\int_{a_{1}}^{a_{2}}\pi_{i}(R)H(R)\pi_{j}(R)dR\cong H\left(R_{i}\right)\delta_% {ij}∫ start_POSTSUBSCRIPT italic_a start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_a start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT end_POSTSUPERSCRIPT italic_π start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( italic_R ) italic_H ( italic_R ) italic_π start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT ( italic_R ) italic_d italic_R ≅ italic_H ( italic_R start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) italic_δ start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT (25)

for integrating smooth function H(R)𝐻𝑅H(R)italic_H ( italic_R ) with quadrature method. Under this approximation, the expansion coefficients cnν,νsubscript𝑐𝑛𝜈superscript𝜈c_{n\nu,\nu^{\prime}}italic_c start_POSTSUBSCRIPT italic_n italic_ν , italic_ν start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT satisfy the following equation:

n,μTnnOnν,nμcnμ,ν+[Uν(Rn)E]cnν,ν=0,subscript𝑛𝜇subscript𝑇𝑛superscript𝑛subscript𝑂𝑛𝜈superscript𝑛𝜇subscript𝑐superscript𝑛𝜇superscript𝜈delimited-[]subscript𝑈𝜈subscript𝑅𝑛𝐸subscript𝑐𝑛𝜈superscript𝜈0\sum_{n,\mu}T_{nn^{\prime}}O_{n\nu,n^{\prime}\mu}c_{n^{\prime}\mu,\nu^{\prime}% }+\left[U_{\nu}\left(R_{n}\right)-E\right]c_{n\nu,\nu^{\prime}}=0,∑ start_POSTSUBSCRIPT italic_n , italic_μ end_POSTSUBSCRIPT italic_T start_POSTSUBSCRIPT italic_n italic_n start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT italic_O start_POSTSUBSCRIPT italic_n italic_ν , italic_n start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT italic_μ end_POSTSUBSCRIPT italic_c start_POSTSUBSCRIPT italic_n start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT italic_μ , italic_ν start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT + [ italic_U start_POSTSUBSCRIPT italic_ν end_POSTSUBSCRIPT ( italic_R start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT ) - italic_E ] italic_c start_POSTSUBSCRIPT italic_n italic_ν , italic_ν start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT = 0 , (26)

with kinetic energy matrix elements:

Tnn=𝑑Rπn(R)[12μ4b2R2]πn(R)𝑑Rsubscript𝑇𝑛superscript𝑛differential-d𝑅subscript𝜋𝑛𝑅delimited-[]12subscript𝜇4bsuperscript2superscript𝑅2subscript𝜋superscript𝑛𝑅differential-d𝑅T_{nn^{\prime}}=\int dR\pi_{n}(R)\left[-\frac{1}{2\mu_{{\rm 4b}}}\frac{% \partial^{2}}{\partial R^{2}}\right]\pi_{n^{\prime}}(R)dRitalic_T start_POSTSUBSCRIPT italic_n italic_n start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT = ∫ italic_d italic_R italic_π start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT ( italic_R ) [ - divide start_ARG 1 end_ARG start_ARG 2 italic_μ start_POSTSUBSCRIPT 4 roman_b end_POSTSUBSCRIPT end_ARG divide start_ARG ∂ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG ∂ italic_R start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ] italic_π start_POSTSUBSCRIPT italic_n start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT ( italic_R ) italic_d italic_R (27)

and overlap matrix:

Onν,nμ=𝑑ΩΦν(Ω;Rn)Φμ(Ω;Rn).subscript𝑂𝑛𝜈superscript𝑛𝜇differential-dΩsubscriptΦ𝜈superscriptΩsubscript𝑅𝑛subscriptΦ𝜇Ωsubscript𝑅superscript𝑛O_{n\nu,n^{\prime}\mu}=\int d\Omega\Phi_{\nu}\left(\Omega;R_{n}\right)^{*}\Phi% _{\mu}\left(\Omega;R_{n^{\prime}}\right).italic_O start_POSTSUBSCRIPT italic_n italic_ν , italic_n start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT italic_μ end_POSTSUBSCRIPT = ∫ italic_d roman_Ω roman_Φ start_POSTSUBSCRIPT italic_ν end_POSTSUBSCRIPT ( roman_Ω ; italic_R start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT roman_Φ start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT ( roman_Ω ; italic_R start_POSTSUBSCRIPT italic_n start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT ) . (28)

Once the expansion coefficients cnν,νsubscript𝑐𝑛𝜈superscript𝜈c_{n\nu,\nu^{\prime}}italic_c start_POSTSUBSCRIPT italic_n italic_ν , italic_ν start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT are determined, the total wave function is fully specified, allowing access to all relevant information for quantum scattering analysis.

In practice, we compute the R𝑅Ritalic_R-matrix, which serves a similar role to the inverse logarithmic derivative of the wave function in single-channel scattering problems. It is defined as

(R)=F(R)[F~(R)]1,𝑅𝐹𝑅superscriptdelimited-[]~𝐹𝑅1\mathcal{R}(R)=F(R)[\tilde{F}(R)]^{-1},caligraphic_R ( italic_R ) = italic_F ( italic_R ) [ over~ start_ARG italic_F end_ARG ( italic_R ) ] start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT , (29)

where F(R)𝐹𝑅F(R)italic_F ( italic_R ) and F~(R)~𝐹𝑅\tilde{F}(R)over~ start_ARG italic_F end_ARG ( italic_R ) denote the wave function and its derivative, respectively, evaluated at a given hyperradius R𝑅Ritalic_R. At large distances, the physical scattering matrix 𝒮𝒮\mathcal{S}caligraphic_S (and the related reaction matrix 𝒦𝒦\mathcal{K}caligraphic_K) can be obtained by applying asymptotic boundary conditions. For dimer-dimer scattering in the lowest channel, we use:

𝒦𝒦\displaystyle\mathcal{K}caligraphic_K =(ff)(gg)1,absent𝑓superscript𝑓superscript𝑔superscript𝑔1\displaystyle=\left(f-f^{\prime}\mathcal{R}\right)\left(g-g^{\prime}\mathcal{R% }\right)^{-1},= ( italic_f - italic_f start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT caligraphic_R ) ( italic_g - italic_g start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT caligraphic_R ) start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT , (30)
𝒮=(1+i𝒦)(1i𝒦)1,𝒮1𝑖𝒦superscript1𝑖𝒦1\mathcal{S}=(1+i\mathcal{K})(1-i\mathcal{K})^{-1},caligraphic_S = ( 1 + italic_i caligraphic_K ) ( 1 - italic_i caligraphic_K ) start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT , (31)

where f𝑓fitalic_f and g𝑔gitalic_gare the regular and irregular solutions in the lowest symmetric channel:

f=2μ4bπkHcos(kHR)𝑓2subscript𝜇4𝑏𝜋subscript𝑘𝐻subscript𝑘𝐻𝑅f=\sqrt{\frac{2\mu_{4b}}{\pi k_{H}}}\cos\left(k_{H}R\right)italic_f = square-root start_ARG divide start_ARG 2 italic_μ start_POSTSUBSCRIPT 4 italic_b end_POSTSUBSCRIPT end_ARG start_ARG italic_π italic_k start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT end_ARG end_ARG roman_cos ( italic_k start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT italic_R ) (32)
g=2μ4bπkHsin(kHR)𝑔2subscript𝜇4𝑏𝜋subscript𝑘𝐻subscript𝑘𝐻𝑅g=\sqrt{\frac{2\mu_{4b}}{\pi k_{H}}}\sin\left(k_{H}R\right)italic_g = square-root start_ARG divide start_ARG 2 italic_μ start_POSTSUBSCRIPT 4 italic_b end_POSTSUBSCRIPT end_ARG start_ARG italic_π italic_k start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT end_ARG end_ARG roman_sin ( italic_k start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT italic_R ) (33)

and kH=2μ4b(EE11)subscript𝑘𝐻2subscript𝜇4𝑏𝐸subscript𝐸11k_{H}=\sqrt{2\mu_{4b}\left(E-E_{11}\right)}italic_k start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT = square-root start_ARG 2 italic_μ start_POSTSUBSCRIPT 4 italic_b end_POSTSUBSCRIPT ( italic_E - italic_E start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT ) end_ARG, which is the hyperradial momentum with respect to the dimer–dimer threshold energy E11=E1+E1subscript𝐸11subscript𝐸1subscript𝐸1E_{11}=E_{1}+E_{1}italic_E start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT = italic_E start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT + italic_E start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT. The scattering energy between two dimers is then Es=EE11subscript𝐸𝑠𝐸subscript𝐸11E_{s}=E-E_{11}italic_E start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT = italic_E - italic_E start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT, and the dimer-dimer relative momentum is k=2kH𝑘2subscript𝑘𝐻k=\sqrt{2}k_{H}italic_k = square-root start_ARG 2 end_ARG italic_k start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT. The energy-dependent scattering length is computed using Eq. (19).

Figures 6 and 7 present the computed dimer–dimer scattering lengths between ground-state dimers of the sinh\sinhroman_sinh-potential. In Fig. 6 (a), the energy-dependent scattering length aD(k)subscript𝑎𝐷𝑘a_{D}(k)italic_a start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT ( italic_k ) is shown as a function of scattering energy Essubscript𝐸𝑠E_{s}italic_E start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT for s=s=1.5𝑠superscript𝑠1.5s=s^{\prime}=1.5italic_s = italic_s start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT = 1.5 (red squares) and s=s=4𝑠superscript𝑠4s=s^{\prime}=4italic_s = italic_s start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT = 4 (blue circles). The solid curves represent the corresponding analytical solutions, with excellent agreement across all energies. Panel (b) shows the zero-energy scattering length as a function of s𝑠sitalic_s, numerical results (red diamonds) again match the analytical curve closely. In contrast, Fig. 7 explores the more general case where ss𝑠superscript𝑠s\neq s^{\prime}italic_s ≠ italic_s start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT, for a fixed s𝑠sitalic_s and varying ssuperscript𝑠s^{\prime}italic_s start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT. In panel (a), the energy-dependent scattering length is plotted s=0.9superscript𝑠0.9s^{\prime}=0.9italic_s start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT = 0.9 (blue circles), 1.51.51.51.5 (red squares), and 1.91.91.91.9 (purple triangles). As ssuperscript𝑠s^{\prime}italic_s start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT decreases, the repulsion between like particles weakens, allowing the formation of a tetramer. When the tetramer becomes bound, a resonance appears in the scattering length, which diverges as shown in panel (b). This resonance occurs near s0.9superscript𝑠0.9s^{\prime}\approx 0.9italic_s start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ≈ 0.9, signaling the onset of tetramer binding.

Refer to caption
Figure 6: Dimer–dimer scattering length between ground-state dimers for the case s=s𝑠superscript𝑠s=s^{\prime}italic_s = italic_s start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT. (a) Numerical results for the energy-dependent scattering length aD(k)subscript𝑎𝐷𝑘a_{D}(k)italic_a start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT ( italic_k ) as a function of scattering energy Essubscript𝐸𝑠E_{s}italic_E start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT, shown for s=1.5𝑠1.5s=1.5italic_s = 1.5 (red squares) and s=4𝑠4s=4italic_s = 4 (blue circles). Solid curves represent the corresponding analytical solutions. (b) Zero-energy scattering length aDsubscript𝑎𝐷a_{D}italic_a start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT as a function of s𝑠sitalic_s. Red diamonds indicate numerical results; the solid curve shows the analytical prediction.
Refer to caption
Figure 7: Dimer–dimer scattering length between ground-state dimers for s=1.5𝑠1.5s=1.5italic_s = 1.5 and varying ssuperscript𝑠s^{\prime}italic_s start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT. Solid curves are provided as visual guides; analytical results are not available for ss𝑠superscript𝑠s\neq s^{\prime}italic_s ≠ italic_s start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT. (a) Energy-dependent scattering length aD(k)subscript𝑎𝐷𝑘a_{D}(k)italic_a start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT ( italic_k ) as a function of scattering energy Essubscript𝐸𝑠E_{s}italic_E start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT, for s=0.9superscript𝑠0.9s^{\prime}=0.9italic_s start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT = 0.9 (blue circles), 1.51.51.51.5 (red squares), and 1.91.91.91.9 (purple triangles). (b) Zero-energy scattering length aDsubscript𝑎𝐷a_{D}italic_a start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT as a function of ssuperscript𝑠s^{\prime}italic_s start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT, computed numerically.

4 Summary

In this work, we developed a robust and numerically accurate framework for studying four-body quantum scattering in one-dimensional systems using the adiabatic hyperspherical representation combined with the slow variable discretization method. This approach enables precise resolution of dimer-dimer scattering properties, even in the presence of sharp avoided crossings in hyperspherical potentials.

We validated our method against known analytical results for sinh\sinhroman_sinh-cosh\coshroman_cosh model interactions, finding excellent agreement across energy scales. Extending beyond analytically tractable regimes, we explored scattering in non-integrable systems and identified a resonance in the dimer-dimer scattering length indicative of tetramer formation. These results not only affirm the accuracy of our method but also demonstrate its potential to uncover emergent few-body phenomena in low-dimensional systems.

Our findings pave the way for future investigations into more complex configurations, such as mass-imbalanced mixtures, excited-state collisions, and recombination processes. The hyperspherical approach, enhanced by SVD, offers a powerful platform for quantitatively probing universal physics in ultracold atomic gases.

\bmhead

Acknowledgements J. W. acknowledges Nirav Mehta for useful discussions. This research was supported by the Australian Research Council’s (ARC) Discovery Program, Grants No. FT230100229 (J. W.), No. DP240101590 (H.H.), and No. DP240100248 (X.-J. L.).

References

  • \bibcommenthead
  • Metcalf and Straten [1999] Metcalf, H.J., Straten, P.: Laser Cooling and Trapping. Springer, New York (1999)
  • Dalfovo et al. [1999] Dalfovo, F., Giorgini, S., Pitaevskii, L.P., Stringari, S.: Theory of bose-einstein condensation in trapped gases. Rev. Mod. Phys. 71, 463–512 (1999)
  • Leggett [2001] Leggett, A.J.: Bose-einstein condensation in the alkali gases: Some fundamental concepts. Rev. Mod. Phys. 73, 307–356 (2001)
  • Giorgini et al. [2008] Giorgini, S., Pitaevskii, L.P., Stringari, S.: Theory of ultracold atomic fermi gases. Rev. Mod. Phys. 80, 1215–1274 (2008)
  • Anderson et al. [1995] Anderson, M.H., Ensher, J.R., Matthews, M.R., Wieman, C.E., Cornell, E.A.: Observation of bose-einstein condensation in a dilute atomic vapor. Science 269, 198 (1995)
  • DeMarco and Jin [1999] DeMarco, B., Jin, D.S.: Onset of fermi degeneracy in a trapped atomic gas. Science 285, 1703 (1999)
  • Bloch et al. [2008] Bloch, I., Dalibard, J., Zwerger, W.: Many-body physics with ultracold gases. Rev. Mod. Phys. 80, 885–964 (2008)
  • Chin et al. [2010] Chin, C., Grimm, R., Julienne, P., Tiesinga, E.: Feshbach resonances in ultracold gases. Rev. Mod. Phys. 82, 1225–1286 (2010)
  • DeMarco and Jin [1970] DeMarco, B., Jin, D.S.: Weakly bound states of three resonantly interacting particles. Yad. Fiz. 12, 1080 (1970)
  • Kraemer et al. [2006] Kraemer, T., Mark, M., Waldburger, P., Danzl, J.G., Chin, C., Engeser, B., Lange, A.D., Pilch, K., Jaakkola, A., Nägerl, H.-C., Grimm, R.: Evidence for efimov quantum states in an ultracold gas of caesium atoms. Nature 440, 315–318 (2006)
  • Ottenstein et al. [2008] Ottenstein, T.B., Lompe, T., Kohnen, M., Wenz, A.N., Jochim, S.: Collisional stability of a three-component degenerate fermi gas. Phys. Rev. Lett. 101, 203202 (2008)
  • Zaccanti et al. [2009] Zaccanti, M., Deissler, B., D’Errico, C., Fattori, M., JonaLasinio, M., Müller, S., Roati, G., Inguscio, M., Modugno, G.: Observation of an efimov spectrum in an atomic system. Nat. Phys. 5, 586–591 (2009)
  • Pollack et al. [2009] Pollack, S.E., Dries, D., Hulet, R.G.: Universality in three- and four-body bound states of ultracold atoms. Science 326, 1683 (2009)
  • Gross et al. [2009] Gross, N., Shotan, Z., Kokkelmans, S., Khaykovich, L.: Observation of universality in ultracold Li7superscriptLi7{}^{7}\mathrm{Li}start_FLOATSUPERSCRIPT 7 end_FLOATSUPERSCRIPT roman_Li three-body recombination. Phys. Rev. Lett. 103, 163202 (2009)
  • Huckans et al. [2009] Huckans, J.H., Williams, J.R., Hazlett, E.L., Stites, R.W., O’Hara, K.M.: Three-body recombination in a three-state fermi gas with widely tunable interactions. Phys. Rev. Lett. 102, 165302 (2009)
  • Williams et al. [2009] Williams, J.R., Hazlett, E.L., Huckans, J.H., Stites, R.W., Zhang, Y., O’Hara, K.M.: Evidence for an excited-state efimov trimer in a three-component fermi gas. Phys. Rev. Lett. 103, 130404 (2009)
  • Nakajima et al. [2010] Nakajima, S., Horikoshi, M., Mukaiyama, T., Naidon, P., Ueda, M.: Nonuniversal efimov atom-dimer resonances in a three-component mixture of Li6superscriptLi6{}^{6}\mathrm{Li}start_FLOATSUPERSCRIPT 6 end_FLOATSUPERSCRIPT roman_Li. Phys. Rev. Lett. 105, 023201 (2010)
  • Barontini et al. [2009] Barontini, G., Weber, C., Rabatti, F., Catani, J., Thalhammer, G., Inguscio, M., Minardi, F.: Observation of heteronuclear atomic efimov resonances. Phys. Rev. Lett. 103, 043201 (2009)
  • Gross et al. [2010] Gross, N., Shotan, Z., Kokkelmans, S., Khaykovich, L.: Nuclear-spin-independent short-range three-body physics in ultracold atoms. Phys. Rev. Lett. 105, 103203 (2010)
  • Lompe et al. [2010] Lompe, T., Ottenstein, T.B., Serwane, F., Viering, K., Wenz, A.N., Zürn, G., Jochim, S.: Atom-dimer scattering in a three-component fermi gas. Phys. Rev. Lett. 105, 103201 (2010)
  • Berninger et al. [2011] Berninger, M., Zenesini, A., Huang, B., Harm, W., Nägerl, H.-C., Ferlaino, F., Grimm, R., Julienne, P.S., Hutson, J.M.: Universality of the three-body parameter for efimov states in ultracold cesium. Phys. Rev. Lett. 107, 120401 (2011)
  • Wild et al. [2012] Wild, R.J., Makotyn, P., Pino, J.M., Cornell, E.A., Jin, D.S.: Measurements of tan’s contact in an atomic bose-einstein condensate. Phys. Rev. Lett. 108, 145305 (2012)
  • Wang et al. [2012a] Wang, J., D’Incao, J.P., Esry, B.D., Greene, C.H.: Origin of the three-body parameter universality in efimov physics. Phys. Rev. Lett. 108, 263001 (2012)
  • Wang et al. [2012b] Wang, Y., Wang, J., D’Incao, J.P., Greene, C.H.: Universal three-body parameter in heteronuclear atomic systems. Phys. Rev. Lett. 109, 243201 (2012)
  • D’Incao et al. [2013] D’Incao, J.P., Wang, J., Esry, B.D., Greene, C.H.: The universality of the efimov three-body parameter. Few-body Syst 54, 1523–1527 (2013)
  • Nishida et al. [2013] Nishida, Y., Moroz, S., Son, D.T.: Super efimov effect of resonantly interacting fermions in two dimensions. Phys. Rev. Lett. 110, 235301 (2013)
  • Gridnev [2014] Gridnev, D.K.: Three resonating fermions in flatland: proof of the super efimov effect and the exact discrete spectrum asymptotics. J. Phys. A 47, 505204 (2014)
  • Gao et al. [2015] Gao, C., Wang, J., Yu, Z.: Revealing the origin of super efimov states in the hyperspherical formalism. Phys. Rev. A 92, 020504 (2015)
  • Nishida [2017] Nishida, Y.: Semisuper efimov effect of two-dimensional bosons at a three-body resonance. Phys. Rev. Lett. 118, 230601 (2017)
  • Esry et al. [1996] Esry, B.D., Lin, C.D., Greene, C.H.: Adiabatic hyperspherical study of the helium trimer. Phys. Rev. A 54, 394–401 (1996)
  • Suno et al. [2002] Suno, H., Esry, B.D., Greene, C.H., Burke, J.P.: Three-body recombination of cold helium atoms. Phys. Rev. A 65, 042725 (2002)
  • D’Incao et al. [2009] D’Incao, J.P., Greene, C.H., Esry, B.D.: Adiabatic hyperspherical study of the helium trimer. J. Phys. B. 42, 044016 (2009)
  • Wang et al. [2011] Wang, J., D’Incao, J.P., Greene, C.H.: Numerical study of three-body recombination for systems with many bound states. Phys. Rev. A 84, 052721 (2011)
  • Wang et al. [2012] Wang, J., D’Incao, J.P., Wang, Y., Greene, C.H.: Universal three-body recombination via resonant d𝑑ditalic_d-wave interactions. Phys. Rev. A 86, 062511 (2012)
  • Chen and Greene [2022] Chen, Y.-H., Greene, C.H.: Efimov physics implications at p𝑝pitalic_p-wave fermionic unitarity. Phys. Rev. A 105, 013308 (2022)
  • Higgins and Greene [2025] Higgins, M.D., Greene, C.H.: Resonances and collisional properties of neutron-rich helium isotopes in the adiabatic hyperspherical representation. Phys. Rev. C 111, 014001 (2025)
  • Esry et al. [1999] Esry, B.D., Greene, C.H., Burke, J.P.: Recombination of three atoms in the ultracold limit. Phys. Rev. Lett. 83, 1751–1754 (1999)
  • Wang et al. [2011] Wang, Y., D’Incao, J.P., Greene, C.H.: Efimov effect for three interacting bosonic dipoles. Phys. Rev. Lett. 106, 233201 (2011)
  • Mehta et al. [2007] Mehta, N.P., Rittenhouse, S.T., D’Incao, J.P., Greene, C.H.: Hyperspherical Approach to the Four-Body Problem. arXiv:0706.1296 (2007)
  • D’Incao et al. [2009] D’Incao, J.P., Rittenhouse, S.T., Mehta, N.P., Greene, C.H.: Dimer-dimer collisions at finite energies in two-component fermi gases. Phys. Rev. A 79, 030501 (2009)
  • Rittenhouse et al. [2011] Rittenhouse, S.T., Stecher, J., D’Incao, J.P., Mehta, N.P., Greene, C.H.: The hyperspherical four-fermion problem. J. Phys. B 44, 172001 (2011)
  • Tolstikhin et al. [1996] Tolstikhin, O.I., Watanabe, S., Matsuzawa, M.: "slow" variable discretization: a novel approach for hamiltonians allowing adiabatic separation of variables. J. Phys. B 29, 389 (1996)
  • Wang [2004] Wang, J.: Hyperspherical approach to quantal three-bodytheory. Phd thesis, University of Colorado, Boulder (2004)
  • Sutherland [2004] Sutherland, B.: Beautiful Models: 70 Years of Exactly Solved Quantum Many-Body Problems. World Scientific Publishing Company, Singapore (2004)