2+1 system in 1D

This example reproduces the results in the article [2]. It studies a one-dimensional 2+1 system of two identical particles interacting with a third particle via a central potential. Here, we consider the interaction to be either a contact interaction or a Gaussian potential, which supports a weakly-bound ground state. The two identical particles do not interact.

Setup

using Printf, FewBodyToolkit

This function returns the universal energy ratios of Table I in the article for a given mass ratio:

function exfun(mr)
    if mr == 2.2
        return [-2.1966, -1.0520]
    elseif mr == 12.4
        return [-2.5963, -1.4818, -1.1970, -1.0377, -1.0002]
    elseif mr == 22.2
        return [-2.7515, -1.6904, -1.3604, -1.1479, -1.0525, -1.0040]
    else
        error("Unknown mass ratio: $mr")
    end
end;

Two-body inverse problem

First we define the 2+1 system via the mass ratio. Here we choose 22.2 since it results in the most bound states. Feel free to change it to either 2.2 or 12.4.

massratio = 22.2 # other values are 2.2 and 12.4
mass_arr = [1.0, massratio, massratio]
mur = 1/(1/mass_arr[1]+1/mass_arr[2]) # reduced mass
println("Mass ratio: ", massratio, ", reduced mass: ", round(mur,digits=4))

Mass ratio: 22.2, reduced mass: 0.9569

We need a potential whose ground state is weakly bound. Here, we use a Gaussian potential and set the target energy to -10^-3. Then, we find the required potential strength via the GEM2B code.

v0 = -1.0; mu_g = 1.0;
vg = GaussianPotential(v0,mu_g)

phys_params2B = make_phys_params2B(;mur,vint_arr=[vg],dim=1)
num_params2B = make_num_params2B(;gem_params=(;nmax=16, r1=1.0, rnmax=120.0))

stateindex = 1; target_e2 = -1e-3;
println("1. Two-body inverse problem")
pps,nps,vscale = GEM2B.v0GEMOptim(phys_params2B,num_params2B,stateindex,target_e2)
println("Found potential parameters: v0 = ", vscale*v0, ", mu_g = ", mu_g)

1. Two-body inverse problem
Found potential parameters: v0 = -0.026696322347653095, mu_g = 1.0

We define the rescaled potential and corresponding physical parameters

vgscaled = GaussianPotential(v0*vscale,mu_g)
pps = make_phys_params2B(;mur,vint_arr=[vgscaled],dim=1)

(hbar = 1.0, mur = 0.9568965517241379, vint_arr = GaussianPotential[GaussianPotential(-0.026696322347653095, 1.0)], lmin = 0, lmax = 0, dim = 1)

Check if the two-body system indeed has the desired binding energy:

e2s = GEM2B.GEM2B_solve(pps,nps)
println("Two-body binding energy: ", e2s[1], " (target: $(target_e2) )")

Two-body binding energy: -0.0010000000110296809 (target: -0.001 )

We can also perform the calculation with a contact potential:

vc = ContactPotential1D(-sqrt(-2*target_e2),0.0)
ppc = make_phys_params2B(;mur,vint_arr=[vc],dim=1)
npc = make_num_params2B(;gem_params=(;nmax=16, r1=1.0, rnmax=120.0))

(gem_params = (nmax = 16, r1 = 1.0, rnmax = 120.0), theta_csm = 0.0, omega_cr = 0.9, threshold = 1.0e-8)

For a contact interaction we don't need to find the potential strength, but parameters should be optimized

r1cs,rnmaxcs,e2copt=GEM_Optim_2B(ppc,npc,stateindex)

1×3 Matrix{Float64}:
 0.00283203  83.5395  -0.000956816

Three-body problem

Two identical bosons

Inputs

Having found the potential parameters, we can now set up the three-body problem with the scaled potential. For bosons we can use their symmetry with the argument svals=["x","b","b"], where x is the different particle and b denotes two identical bosons.

vint_arr=[[],[vgscaled],[vgscaled]] #[[v23],[v31],[v12]]
phys_params3B = make_phys_params3B1D(;mass_arr=mass_arr,svals=["x","b","b"],vint_arr)

(hbar = 1.0, mass_arr = [1.0, 22.2, 22.2], svals = ["x", "b", "b"], vint_arr = Vector{Any}[[], [GaussianPotential(-0.026696322347653095, 1.0)], [GaussianPotential(-0.026696322347653095, 1.0)]], parity = 1)

For the numerical parameters nmax, r1, and rnmax we use the optimized ones, found by the two-body inverse problem. The parameters for the other Jacobi coordinates are set manually.

nmax = nps.gem_params.nmax;
r1 = nps.gem_params.r1; rnmax = nps.gem_params.rnmax;
num_params3B = make_num_params3B1D(;gem_params=(;nmax, r1, rnmax, Nmax=16, R1=1.5, RNmax=250.0))

(lmin = 0, Lmin = 0, lmax = 0, Lmax = 0, gem_params = (nmax = 16, r1 = 0.7375, rnmax = 150.0125, Nmax = 16, R1 = 1.5, RNmax = 250.0), theta_csm = 0.0, omega_cr = 0.9, kmax_interpol = 1000, threshold = 1.0e-8)

Solving the three-body problem

println("\n2. Solving the three-body problem and comparing to the article's results")
e3 = GEM3B1D.GEM3B1D_solve(phys_params3B,num_params3B);

2. Solving the three-body problem and comparing to the article's results

We compute the ratio of three-body to two-body binding energies.

println("Results for the bosonic case, mass ratio = $massratio")
epsilon = e3 /abs(e2s[1])

ex_arr = exfun(massratio)[1:2:end]
comparison(epsilon, ex_arr, min(length(epsilon),length(ex_arr)); s1="Gaussian", s2="Contact (Ref.)")

Results for the bosonic case, mass ratio = 22.2
Index   Gaussian        Contact (Ref.)  Difference
1       -2.742642       -2.751500       -0.008858
2       -1.360521       -1.360400       0.000121
3       -1.052605       -1.052500       0.000105

We can also do the calculation with a contact potential:

vint_arrC=[[],[vc],[vc]]
pp3BC = make_phys_params3B1D(;mass_arr=mass_arr,svals=["x","b","b"],vint_arr = vint_arrC)
nmax = nps.gem_params.nmax;
np3BC = make_num_params3B1D(;gem_params=(;nmax, r1=r1cs, rnmax=rnmaxcs, Nmax=16, R1=1.5, RNmax=250.0))
println("\nThree-body problem with contact potential")
e3c = GEM3B1D.GEM3B1D_solve(pp3BC,np3BC);
epsilonC = e3c /abs(e2copt);

comparison(epsilonC, ex_arr, min(length(epsilonC),length(ex_arr)); s1="Contact (FBTK)", s2="Contact (Ref.)")

Three-body problem with contact potential
Index   Contact (FBTK)  Contact (Ref.)  Difference
1       -2.751575       -2.751500       0.000075
2       -1.360447       -1.360400       0.000047
3       -1.052559       -1.052500       0.000059

Two identical fermions

Inputs

For fermions we can use the same potential. To account for their different statistics and parity, we use svals=["x","f","f"] and parity=-1. To allow for basis functions that obey these requirements, we need to set lmax, Lmax=1.

println("\n3. Results for the fermionic case, mass ratio = $massratio")
phys_params3B_F = make_phys_params3B1D(;mass_arr=mass_arr,svals=["x","f","f"],vint_arr,parity=-1)
num_params3B_F = make_num_params3B1D(;gem_params=(;nmax, r1, rnmax, Nmax=16, R1=1.5, RNmax=250.0), lmin=0, Lmin=0, lmax=1, Lmax=1)
e3_F = GEM3B1D.GEM3B1D_solve(phys_params3B_F,num_params3B_F);

3. Results for the fermionic case, mass ratio = 22.2

We compute the ratio of three-body to two-body binding energies.

epsilon_F = e3_F /abs(e2s[1])

ex_arr_F = exfun(massratio)[2:2:end]
comparison(epsilon_F, ex_arr_F, min(length(epsilon_F),length(ex_arr_F)); s1="Gaussian", s2="Contact (Ref.)")

Index   Gaussian        Contact (Ref.)  Difference
1       -1.694973       -1.690400       0.004573
2       -1.149204       -1.147900       0.001304
3       -1.004203       -1.004000       0.000203

Again, we can also do the calculation with a contact potential:

pp3BC = make_phys_params3B1D(;mass_arr=mass_arr,svals=["x","f","f"],vint_arr = vint_arrC, parity=-1)
np3BC = make_num_params3B1D(;gem_params=(;nmax, r1=r1cs, rnmax=rnmaxcs, Nmax=16, R1=1.5, RNmax=250.0),lmin=0,Lmin=0,lmax=1,Lmax=1)
println("\nThree-body problem with contact potential")
e3c_F = GEM3B1D.GEM3B1D_solve(pp3BC,np3BC);
epsilonC_F = e3c_F /abs(e2copt);

comparison(epsilonC_F, ex_arr_F, min(length(epsilonC_F),length(ex_arr_F)); s1="Contact (FBTK)", s2="Contact (Ref.)")

Three-body problem with contact potential
Index   Contact (FBTK)  Contact (Ref.)  Difference
1       -1.690480       -1.690400       0.000080
2       -1.147950       -1.147900       0.000050
3       -1.004036       -1.004000       0.000036

Overall, we can reproduce the article's results quite well for both bosonic and fermionic systems. For the contact interaction, the results match within the provided accuracy. For the Gaussian interaction, better results could be obtained with more basis functions and/or optimized basis parameters. Note, however, that perfect agreement between the finite-range Gaussian interaction and the contact interaction cannot be reached. Only in the limit of vanishing two-body binding energy, the two potentials should yield the same results.

Page References

See also the full bibliography for further references cited throughout this documentation.

This page was generated using Literate.jl.