API

CentralPotential(f::Function)

A concrete implementation of PotentialFunction that represents a central potential. It takes a function f that defines the potential as a function of the radial distance r.

ContactPotential1D(v0::Float64, z0::Float64)

A concrete implementation of PotentialFunction that represents a 1D Contact (Dirac) potential:

\[V(z) = v_0 \delta(z - z_0)\]

where z is the 1D coordinate.

Arguments:

• v0::Float64: The strength of the potential.
• z0::Float64: The position of the delta function.

function (gp::ContactPotential1D)(z::Float64)

Evaluates the contact potential at a given 1D coordinate z.

GaussianPotential(v0::Float64, mu_g::Float64)

A concrete implementation of PotentialFunction that represents a Gaussian potential:

\[V(r) = v_0 e^{-\mu_g r^2}\]

where r is the radial distance.

Arguments:

• v0::Float64: The strength of the potential.
• mu_g::Float64: The width parameter of the Gaussian potential.

function (gp::GaussianPotential)(r::Float64)

Evaluates the Gaussian potential at a given radial distance r.

PotentialFunction

Abstract type for potential functions used in few-body calculations. This type serves as a parent type for more specific potential implementations

PowerLawPotential(v0::Float64, p::Float64)

A concrete implementation of PotentialFunction that represents a power-law potential:

\[V(r) = v_0 |r|^{p}\]

Note the absolute value: the potential is defined via |r|, i.e. it is an even function. This matters only in 1D, where the coordinate can be negative; there, PowerLawPotential(v0,p) describes v0*|x|^p. Odd potentials are not supported. In 2D and 3D, r >= 0 anyway and |r| = r.

Matrix elements are evaluated analytically (closed-form radial integration), so neither numerical integration nor range-interpolation is needed. Prominent cases are the Coulomb interaction (p=-1) and the harmonic oscillator (p=2).

The matrix elements are finite only if p is not too negative. The bound depends on the dimension and on the angular momentum, and is therefore not checked here but in the individual modules:

• GEM2B: requires p > -2*lmax - dim, e.g. p > -1 for lmax=0 in 1D, so a pure 1/|x| is already divergent in 1D.
• ISGL: requires p > -3.

Arguments:

• v0::Float64: The strength of the potential.
• p::Float64: The exponent of the power law.

function (pp::PowerLawPotential)(r)

Evaluates the power-law potential at a given radial distance r. Uses abs(r), see the note on 1D in PowerLawPotential.

check_cr_csm_sector(complex_ranged, complex_scaling, complex_range_freq, complex_scaling_angle)

Guards the one combination in which complex ranges and complex scaling do not simply coexist.

Complex scaling is applied here by rotating the integration contour back onto the real axis, so that a user-supplied potential is only ever called with a real argument. The price is that the rotation lands on the Gaussian instead: the radial integrals are taken at alpha*csmfac^2, with csmfac = exp(-i*theta_csm) and alpha anywhere in the sector |arg(alpha)| <= atan(omega). The factor exp(-alpha*csmfac^2*r^2) decays only while atan(omega) + 2*theta_csm < pi/2, and the corner of the sector where this is tightest (arg(alpha) = -atan(omega)) is attained exactly.

This is a limit of the representation, not of the matrix element: the same quantity written without the contour rotation, int v(r*exp(i*theta_csm)) r^n exp(-alpha*r^2) dr, converges throughout (both forms agree to full precision wherever the rotated one is computable). Lifting the limit would mean evaluating the potential at complex arguments, which not every user-supplied function supports.

The tail of the potential does not rescue the rotated form. For a slowly decaying potential the integral genuinely diverges; for a fast decaying one it converges mathematically but still fails numerically, since v(r) and exp(-alpha*csmfac^2*r^2) are separate factors and underflow times overflow gives NaN. Either way quadgk reports an unrelated-looking DomainError, so the combination is caught here instead.

eigen2step(e_arr, H, S; threshold=1e-13)

Solves the generalized eigenvalue problem H*x = E*S*x and writes the energies into e_arr.

Similar to eigvals(H,S), but the basis is first orthonormalized by reduce_basis, which cuts the small eigenvalues of S. The number of energies can therefore be smaller than the size of H.

eigen2step_valvec(e_arr, v_arr, H, S; threshold=1e-13)

Same as eigen2step, but also writes the coefficient vectors into the columns of v_arr, normalized as x' * S * x == 1.

inverse_solve(T, V, S; target_energy=0.0, threshold=1e-13, return_vectors=false)

Solves the inverse problem: instead of the energies at a given interaction strength, it returns the interaction strengths lambda at which a state of energy target_energy exists.

For a Hamiltonian which is linear in a strength parameter, H = T + lambda*V, the condition that target_energy is an eigenvalue of H is the generalized eigenvalue problem

(T - target_energy*S) x = -lambda * V x

in the non-orthogonal basis with norm-overlap S. One matrix build and one eigen solve therefore replace a whole scan over the strength with bracketing of the energy.

Arguments

• T: kinetic energy matrix (not the Hamiltonian T+V), as returned by GEM2B_matrices, GEM3B1D_matrices, ISGL_matrices.
• V: interaction matrix at unit strength. lambda is the factor multiplying all interactions, so pass the potential at unit strength.
• S: norm-overlap matrix.

Keywords

• target_energy=0.0: the energy that should be reached. Must lie below the lowest continuum threshold of the basis, such that T - target_energy*S is positive definite; otherwise an error is raised. target_energy=0.0 gives the critical strengths at which a state becomes bound.
• threshold=1e-13: cut-off for the eigenvalues of S, see reduce_basis. Results near target_energy=0.0 are sensitive to this value: too small a cut-off can produce spurious, nearly linearly dependent solutions, too large a one cuts the diffuse basis functions needed close to threshold.
• return_vectors=false: whether to also return the coefficient vectors of the corresponding states.

Real-symmetric and hermitian matrices are supported, i.e. both the usual and the complex-ranged basis functions. The complex scaling method is not: it makes T and V complex-symmetric rather than hermitian, and is rejected.

Returns

• lambdas: the positive strengths, in ascending order. The first entry is the smallest strength at which the system supports a state of energy target_energy. Negative strengths are discarded; if they are of interest, call the function with -V instead.
• vectors: (optional) the corresponding coefficient vectors, normalized as x' * S * x == 1.

Example

T,V,S = GEM2B_matrices(phys_params, num_params)
lambdas = inverse_solve(T,V,S; target_energy=0.0, threshold=1e-8) # critical strengths for binding

parse_complex_ranged(complex_ranged) -> (complex_ranged_r, complex_ranged_R)

Translates the user-facing complex_ranged option into two booleans, one per Jacobi coordinate.

Accepted values are the symbols :none, :r, :R, :both, and a Bool (true meaning :both, false meaning :none) for consistency with GEM2B_solve, whose basis has only one coordinate.

Used by both three-body modules (ISGL and GEM3B1D).

quadgk_scaled(f, domain::Tuple; segbuf=nothing)

quadgk over domain (e.g. (-Inf,0,Inf)), with a fallback for integrals that vanish or nearly cancel (e.g. odd moments of even potentials, needed for parity=0 in 1D). With its purely relative default tolerance, quadgk cannot converge on those and runs until maxevals (10^7 evaluations); the fallback measures the error relative to ∫|f| instead. Integrals that converge normally are computed as before.

reduce_basis(S; threshold=1e-13)

Transformation to an orthonormalized, possibly reduced basis: returns the matrix L (size n x m, m <= n) with L' * S * L == I.

L is built from the eigenvectors of the norm-overlap matrix S, dropping all eigenvalues smaller than threshold times the largest one. This removes the near-linear dependence of the non-orthogonal basis; m < n whenever some eigenvalues are cut. Both real-symmetric S and hermitian S (complex-ranged basis functions) are supported.

Used by eigen2step and inverse_solve; useful on its own to transform a generalized eigenvalue problem H*x = E*S*x into the ordinary one (L'*H*L)*y = E*y, with x = L*y.

PreallocStruct2B{TTV, TS, TE}

A structure that holds preallocated arrays for the GEM-2B solver, used to store basis parameters, matrices, and results for two-body calculations.

Fields

• nu_arr::Vector{TTV}: Array of nonlinear variational parameters (basis exponents).
• S::Matrix{TS}: Overlap matrix between basis functions.
• T::Matrix{TTV}: Kinetic energy matrix.
• V::Matrix{TTV}: Potential energy matrix.
• energies::Vector{TE}: Array to store computed eigenvalues (energies).
• wavefunctions::Matrix{TTV}: Matrix to store eigenvectors (wavefunctions).

Type Parameters

• TTV: Element type for kinetic, potential, and basis parameter arrays (e.g., Float64 or ComplexF64).
• TS: Element type for the overlap matrix (e.g., Float64 or ComplexF64).
• TE: Element type for the energies array (e.g., Float64 or ComplexF64).

Description

This struct is designed to minimize memory allocations and improve performance by reusing arrays during repeated GEM-2B calculations. The types and sizes of the arrays are determined by the numerical parameters and whether complex rotation or complex scaling is used.

Keyword arguments

• num_params: A named tuple containing numerical parameters, including the maximum number of basis functions (nmax).
• complex_ranged: Boolean indicating whether complex range basis functions are used.
• complex_scaling: Boolean indicating whether complex scaling is used.

Example

PreallocStruct2B(num_params, complex_ranged=false, complex_scaling=false) # for real basis functions and no complex scaling
PreallocStruct2B(num_params, complex_ranged=true, complex_scaling=false) # for complex range basis functions and no complex scaling
PreallocStruct2B(num_params, complex_ranged=false, complex_scaling=true) # for real basis functions with complex scaling

GEM2B_matrices(phys_params, num_params; complex_ranged=false, complex_scaling=false)

Builds and returns the matrices of a two-body problem without solving it: (; T, V, S), where T is the kinetic energy (not the Hamiltonian T+V), V the interaction, and S the norm-overlap of the basis functions.

Arguments and keywords are the same as for GEM2B_solve. Intended for problems which are more naturally posed in terms of the matrices than of the energies, in particular the inverse problem, see inverse_solve.

With complex_scaling=true, T already carries the factor exp(-2*i*theta) and V the rotated potential.

Example

T,V,S = GEM2B_matrices(phys_params, num_params)

GEM2B_solve(phys_params, num_params; return_wavefunctions=false, complex_ranged=false, complex_scaling=false)

Solves two-body quantum mechanical problems using the Gaussian Expansion Method (GEM).

Arguments

• phys_params: Physical parameters describing the two-body system:
  • hbar::Float64: reduced Planck constant
  • mur::Float64: reduced mass
  • interactions=Vector{Any}: a vector of interactions
  • lmax::Int: power of r^lmax in the basis functions; indicator for the angular momentum in 3D
  • dim::Int: the spatial dimension
• num_params: Numerical parameters struct containing information on the set of basis functions:
  • gem_params::NamedTuple: (number of basis functions, smallest and largest range parameters).
  • complex_scaling_angle::Float64: Complex scaling angle (in radians) for the Complex Scaling Method.
  • complex_range_freq::Float64: Parameter controlling the frequency for complex-ranged basis functions.
  • threshold::Float64: Numerical threshold generalized eigenvalue solver.

Keywords

• return_wavefunctions=false: Whether to return wavefunctions (false: energies only, true: energies and wavefunctions)
• complex_ranged=false: Whether to use complex-ranged basis functions.
• complex_scaling=false: Whether to use the complex scaling method.

Returns

• energies: Array of energy eigenvalues
• wavefunctions: (Optional) Array of eigenvectors if return_wavefunctions=true

Example

phys_params = make_phys_params2B(hbar=1.0, mur=1.0, interactions=[GaussianPotential(-1.0, 0.5)], lmax=0, dim=3)
num_params = make_num_params2B(gem_params=(nmax=5, r1=1.0, rnmax=10.0), complex_range_freq=0.5, complex_scaling_angle=0.0, threshold=1e-10)
energies = GEM2B_solve(phys_params, num_params; return_wavefunctions=false, complex_ranged=false, complex_scaling=false)
# or with wavefunctions:
energies, wavefunctions = GEM2B_solve(phys_params, num_params; return_wavefunctions=true)
# Note: The function can handle 1D, 2D, or 3D problems based on the `dim` parameter in `phys_params`.

GEM_Optim_2B(phys_params, num_params, stateindex; complex_ranged=0, g_tol=1e-9)

Optimize the ranges used in the Gaussian Expansion Method (GEM) for a specific 2-body state.

Arguments

• phys_params::NamedTuple: Physical parameters such as hbar, mur, interactions, lmax, lmin, and dim.
• num_params::NamedTuple: Numerical parameters such as gem_params, complex_range_freq, complex_scaling_angle, and threshold.
• stateindex::Int: An integer specifying the index of the state to optimize (e.g., 1 for the ground state).

Keyword Arguments

• complex_ranged::Int=0: Indicates whether to use complex rotation.
• g_tol::Float64=1e-9: A value specifying the tolerance for the optimizer.

Returns

• A vector containing:
  • Optimized GEM parameters: [r1, rnmax].
  • The energy value of the optimized state.

Example

r1opt,rnmaxopt,energy = GEM_Optim_2B(phys_params, num_params, 1) # optimize for the ground state
r1opt,rnmaxopt,energy = GEM_Optim_2B(phys_params, num_params, 3; complex_ranged = 1) # optimize for the third state (2nd excited) using complex-ranged basis functions

eigen2step(e_arr, H, S; threshold=1e-13)

Solves the generalized eigenvalue problem H*x = E*S*x and writes the energies into e_arr.

Similar to eigvals(H,S), but the basis is first orthonormalized by reduce_basis, which cuts the small eigenvalues of S. The number of energies can therefore be smaller than the size of H.

eigen2step_valvec(e_arr, v_arr, H, S; threshold=1e-13)

Same as eigen2step, but also writes the coefficient vectors into the columns of v_arr, normalized as x' * S * x == 1.

inverse_solve(T, V, S; target_energy=0.0, threshold=1e-13, return_vectors=false)

Solves the inverse problem: instead of the energies at a given interaction strength, it returns the interaction strengths lambda at which a state of energy target_energy exists.

For a Hamiltonian which is linear in a strength parameter, H = T + lambda*V, the condition that target_energy is an eigenvalue of H is the generalized eigenvalue problem

(T - target_energy*S) x = -lambda * V x

in the non-orthogonal basis with norm-overlap S. One matrix build and one eigen solve therefore replace a whole scan over the strength with bracketing of the energy.

Arguments

• T: kinetic energy matrix (not the Hamiltonian T+V), as returned by GEM2B_matrices, GEM3B1D_matrices, ISGL_matrices.
• V: interaction matrix at unit strength. lambda is the factor multiplying all interactions, so pass the potential at unit strength.
• S: norm-overlap matrix.

Keywords

• target_energy=0.0: the energy that should be reached. Must lie below the lowest continuum threshold of the basis, such that T - target_energy*S is positive definite; otherwise an error is raised. target_energy=0.0 gives the critical strengths at which a state becomes bound.
• threshold=1e-13: cut-off for the eigenvalues of S, see reduce_basis. Results near target_energy=0.0 are sensitive to this value: too small a cut-off can produce spurious, nearly linearly dependent solutions, too large a one cuts the diffuse basis functions needed close to threshold.
• return_vectors=false: whether to also return the coefficient vectors of the corresponding states.

Real-symmetric and hermitian matrices are supported, i.e. both the usual and the complex-ranged basis functions. The complex scaling method is not: it makes T and V complex-symmetric rather than hermitian, and is rejected.

Returns

• lambdas: the positive strengths, in ascending order. The first entry is the smallest strength at which the system supports a state of energy target_energy. Negative strengths are discarded; if they are of interest, call the function with -V instead.
• vectors: (optional) the corresponding coefficient vectors, normalized as x' * S * x == 1.

Example

T,V,S = GEM2B_matrices(phys_params, num_params)
lambdas = inverse_solve(T,V,S; target_energy=0.0, threshold=1e-8) # critical strengths for binding

make_num_params2B(; gem_params=(nmax=5, r1=1.0, rnmax=10.0), complex_scaling_angle=0.0, complex_range_freq=0.9, threshold=1e-8)

Create and return a named tuple containing the numerical parameters for a two-body GEM calculation.

Keyword arguments

• gem_params::NamedTuple = (nmax=5, r1=1.0, rnmax=10.0): Parameters for the Gaussian Expansion Method (number of basis functions, smallest and largest range parameters).
• complex_scaling_angle::Float64 = 0.0: Complex scaling angle (in radians) for the Complex Scaling Method.
• complex_range_freq::Float64 = 0.9: Parameter controlling the frequency for complex-ranged basis functions.
• threshold::Float64 = 1e-8: Numerical threshold generalized eigenvalue solver.

Returns

• NamedTuple: Named tuple with the specified numerical parameters.

Example

make_num_params2B(gem_params=(nmax=20, r1=0.1, rnmax=50.0)) # for a larger basis set
make_num_params2B(gem_params=(nmax=20, r1=0.1, rnmax=50.0), complex_scaling_angle = 10.0) # non-zero rotation angle for complex scaling method (CSM)

make_phys_params2B(; hbar=1.0, mur=1.0, interactions=[[GaussianPotential(-1.0, 1.0)]], lmin=0, lmax=0, dim=3, parity=nothing)

Create and return a named tuple containing the physical parameters for a two-body system.

Keyword arguments

• hbar::Float64 = 1.0: Reduced Planck constant used in calculations.
• mur::Float64 = 1.0: Reduced mass of the two-body system.
• interactions::Vector{Any} = [[GaussianPotential(-1.0, 1.0)]]: Array of interaction potentials or related parameters.
• lmin::Int = 0: Minimum orbital angular momentum quantum number.
• lmax::Int = 0: Maximum orbital angular momentum quantum number.
• dim::Int = 3: Spatial dimension of the system.
• parity = nothing: With parity=0 (only for dim=1), the basis contains both even (l=0) and odd (l=1) functions, as needed for potentials with V(r) != V(-r); lmin, lmax are then not used.

Returns

• NamedTuple: Named tuple with the specified physical parameters.

Example

make_phys_params2B(interactions=[GaussianPotential(-1.0, 0.5)], dim=1) # 1D system with Gaussian potential
make_phys_params2B(mur=0.5, interactions=[r -> -1/r], lmax=2)       # 3D Coulomb potential in d-wave (l=2) with reduced mass 0.5

reduce_basis(S; threshold=1e-13)

Transformation to an orthonormalized, possibly reduced basis: returns the matrix L (size n x m, m <= n) with L' * S * L == I.

L is built from the eigenvectors of the norm-overlap matrix S, dropping all eigenvalues smaller than threshold times the largest one. This removes the near-linear dependence of the non-orthogonal basis; m < n whenever some eigenvalues are cut. Both real-symmetric S and hermitian S (complex-ranged basis functions) are supported.

Used by eigen2step and inverse_solve; useful on its own to transform a generalized eigenvalue problem H*x = E*S*x into the ordinary one (L'*H*L)*y = E*y, with x = L*y.

v0GEMOptim(phys_params, num_params, stateindex, target_e2; complex_ranged=0, rtol=1e-4, atol=10*eps(), g_tol=1e-9, output=false)

Finds a value v0crit to globally scale the potential in order to achieve the target energy target_e2 for the state specified by stateindex. Additionally, performs intermediate optimization of Gaussian Expansion Method (GEM) parameters.

Arguments

• phys_params::NamedTuple: Physical parameters including hbar, mur, interactions, lmax, lmin, and dim.
• num_params::NamedTuple: Numerical parameters including gem_params, complex_range_freq, complex_scaling_angle, and threshold.
• stateindex::Int: Index of the state to optimize, where 1 indicates the ground state.
• target_e2::Float64: Target energy value.

Keyword Arguments

• complex_ranged::Int: Use complex rotation (default: 0).
• rtol::Float64: Relative tolerance for energy convergence (default: 1e-4).
• atol::Float64: Absolute tolerance for energy convergence (default: 10*eps()).
• g_tol::Float64: Tolerance for the optimizer (default: 1e-9).
• output::Bool: If true, prints intermediate optimization results (default: false).

Returns

• phys_params::NamedTuple: Updated physical parameters with optimized potential.
• num_params::NamedTuple: Updated numerical parameters with optimized GEM ranges.
• v0crit::Float64: The value to scale the potential with, to achieve the target energy. This is not the overall value of v0 but rather a scaling factor for the potential.

Example

phys_params_scaled,num_params_optimized,scalingfactor = v0GEMOptim(phys_params, num_params, 1, -2.5)

wavefun_arr(r_arr, phys_params, num_params, wf_arr; complex_ranged=false)

Calculates the 2-body wavefunction at specified positions based on input parameters and coefficients.

Arguments

• r_arr::Vector{Float64}: Array of positions where the wavefunction is evaluated.
• phys_params::NamedTuple: Physical parameters, here we only need lmaxand dim
• num_params::NamedTuple: Numerical parameters containing the information about the Gaussian ranges
• wf_arr::Vector: Eigenvector from the diagonalization routine. Can be complex-valued.

Keyword Arguments

• complex_ranged::Bool=false: Determines whether to use complex-ranged Gaussians.

Returns

• Vector{Float64}: The wavefunction values at the specified positions r_arr.

Example

wavefun_arr(0.0:0.1:10.0, phys_params, num_params, wf_arr; complex_ranged=false)

wavefun_point(r, nu_arr, wf_arr, ll, dim)

Compute the value of the 2-body wavefunction at a given position r.

Arguments

• r::Float64: The radial position where the wavefunction is evaluated.
• nu_arr::Vector: Array of Gaussian widths.
• wf_arr::Vector: Coefficients of the wavefunction for the given state.
• ll::Int: Orbital angular momentum quantum number.
• dim::Int: Dimensionality of the system.

Returns

• Float64: The computed value of the wavefunction at the specified position r.

Example

nu_arr = [1.0, 0.0025]
wf_arr = [0.8,0.6]
psi_value = wavefun_point(1.0, nu_arr, wf_arr, 0, 3)

GEM3B1D_matrices(phys_params, num_params; complex_scaling=false, complex_ranged=:none)

Builds and returns the matrices of a 1D three-body problem without solving it: (; T, V, S), where T is the kinetic energy (not the Hamiltonian T+V), V the interaction, and S the norm-overlap of the basis functions.

Arguments and keywords are the same as for GEM3B1D_solve. Intended for problems which are more naturally posed in terms of the matrices than of the energies, in particular the inverse problem, see inverse_solve.

With complex_scaling=true, T already carries the factor exp(-2*i*theta) and V the rotated potential.

Example

T,V,S = GEM3B1D_matrices(phys_params, num_params)

GEM3B1D_solve(phys_params, num_params; return_wavefunctions=false, complex_scaling=false, complex_ranged=:none, observ_params=(;stateindices=[],centobs_arr=[[],[],[]],R2_arr=[0,0,0]))

Solves the 1D three-body problem using the Gaussian Expansion Method (GEM).

Arguments

• phys_params: Physical parameters for the three-body system (e.g., masses, interaction potentials, etc.).
• num_params: Numerical parameters for the GEM calculation (e.g., basis size, grid parameters, etc.).
• return_wavefunctions: (optional) If true, also returns wavefunction-related observables. Default is false.
• complex_scaling: (optional) If true, uses complex scaling method. Default is false.
• complex_ranged: (optional) Selects complex-ranged basis functions nu -> nu*(1 + i*omega) per Jacobi coordinate, with omega = complex_range_freq from num_params. One of :none (default), :r (only the internal coordinate $r$), :R (only $R$), or :both. A Bool is also accepted, with true meaning :both, for consistency with GEM2B_solve. Each complex-ranged coordinate doubles its number of basis functions, since the ranges enter as conjugate pairs; choosing :both therefore multiplies the total basis size by four. All interaction types are supported, including generic central potentials: their radial integrals are interpolated over the sector of the complex plane in which the effective Gaussian range lives, with kmax_theta angular nodes (see make_num_params3B1D).
• observ_params: (optional) Parameters for observable calculations. Currently not supported for 1D.

Returns

• If return_wavefunctions=false: Returns an array of computed energies.
• If return_wavefunctions=true: Returns a tuple (energies, wavefunctions).

Example

phys_params = make_phys_params3B1D()
num_params = make_num_params3B1D()
energies = GEM3B3D_solve(phys_params, num_params) #solving with default parameters: three particles with the same mass and gaussian interaction

eigen2step(e_arr, H, S; threshold=1e-13)

Solves the generalized eigenvalue problem H*x = E*S*x and writes the energies into e_arr.

Similar to eigvals(H,S), but the basis is first orthonormalized by reduce_basis, which cuts the small eigenvalues of S. The number of energies can therefore be smaller than the size of H.

eigen2step_valvec(e_arr, v_arr, H, S; threshold=1e-13)

Same as eigen2step, but also writes the coefficient vectors into the columns of v_arr, normalized as x' * S * x == 1.

inverse_solve(T, V, S; target_energy=0.0, threshold=1e-13, return_vectors=false)

Solves the inverse problem: instead of the energies at a given interaction strength, it returns the interaction strengths lambda at which a state of energy target_energy exists.

For a Hamiltonian which is linear in a strength parameter, H = T + lambda*V, the condition that target_energy is an eigenvalue of H is the generalized eigenvalue problem

(T - target_energy*S) x = -lambda * V x

in the non-orthogonal basis with norm-overlap S. One matrix build and one eigen solve therefore replace a whole scan over the strength with bracketing of the energy.

Arguments

• T: kinetic energy matrix (not the Hamiltonian T+V), as returned by GEM2B_matrices, GEM3B1D_matrices, ISGL_matrices.
• V: interaction matrix at unit strength. lambda is the factor multiplying all interactions, so pass the potential at unit strength.
• S: norm-overlap matrix.

Keywords

• target_energy=0.0: the energy that should be reached. Must lie below the lowest continuum threshold of the basis, such that T - target_energy*S is positive definite; otherwise an error is raised. target_energy=0.0 gives the critical strengths at which a state becomes bound.
• threshold=1e-13: cut-off for the eigenvalues of S, see reduce_basis. Results near target_energy=0.0 are sensitive to this value: too small a cut-off can produce spurious, nearly linearly dependent solutions, too large a one cuts the diffuse basis functions needed close to threshold.
• return_vectors=false: whether to also return the coefficient vectors of the corresponding states.

Real-symmetric and hermitian matrices are supported, i.e. both the usual and the complex-ranged basis functions. The complex scaling method is not: it makes T and V complex-symmetric rather than hermitian, and is rejected.

Returns

• lambdas: the positive strengths, in ascending order. The first entry is the smallest strength at which the system supports a state of energy target_energy. Negative strengths are discarded; if they are of interest, call the function with -V instead.
• vectors: (optional) the corresponding coefficient vectors, normalized as x' * S * x == 1.

Example

T,V,S = GEM2B_matrices(phys_params, num_params)
lambdas = inverse_solve(T,V,S; target_energy=0.0, threshold=1e-8) # critical strengths for binding

make_num_params3B1D(;lmin=0, Lmin=0, lmax=0, Lmax=0, gem_params=(nmax=5, r1=1.0, rnmax=10.0, Nmax=5, R1=1.0, RNmax=10.0), complex_scaling_angle=0.0, complex_range_freq=0.9, kmax_interpol=1000, kmax_theta=25, threshold=10^-8)

Create and return a named tuple containing the numerical parameters for a three-body GEM calculation in 1D.

Keyword arguments

• lmin::Int = 0: Minimum power r^l used in the basis functions of the $r$ Jacobi coordinate.
• Lmin::Int = 0: Minimum power r^L used in the basis functions of the $R$ Jacobi coordinate.
• lmax::Int = 0: Maximum power r^l used in the basis functions of the $r$ Jacobi coordinate.
• Lmax::Int = 0: Maximum power r^L used in the basis functions of the $R$ Jacobi coordinate.
• gem_params::NamedTuple = (nmax=5, r1=1.0, rnmax=10.0, Nmax=5, R1=1.0, RNmax=10.0): Parameters for the Gaussian Expansion Method (number of basis functions, smallest and largest range parameters for both Jacobi coordinates).
• complex_scaling_angle::Float64 = 0.0: Complex scaling angle (in degrees) for the Complex Scaling Method.
• complex_range_freq::Float64 = 0.9: Parameter $\omega$ controlling the frequency for complex-ranged basis functions, $\nu \to \nu(1 + i\omega)$. Only used when GEM3B1D_solve is called with complex_ranged set to something other than :none. The same $\omega$ is used for both Jacobi coordinates.
• kmax_interpol::Int = 1000: Number of numerical integration with effective Gaussian ranges used for interpolation.
• kmax_theta::Int = 25: Number of angular nodes of the range-interpolation. Only used with complex-ranged basis functions, where the effective Gaussian range becomes complex and the interpolation runs over the sector $|\arg\alpha| \le \arctan\omega$ instead of over a real interval. The setup cost of the interpolation grows in proportion, so kmax_interpol can be reduced in compensation. Increase it to check convergence; complex_ranged=:both is the most demanding case, since the four-fold enlarged basis makes the overlap matrix nearly singular and amplifies any error in the matrix elements.
• threshold::Float64 = 1e-8: Numerical threshold for the generalized eigenvalue solver.

Returns

• NamedTuple: Named tuple with the specified numerical parameters.

Example

make_num_params3B1D() # for the default basis set
make_num_params3B1D(gem_params=(nmax=10, r1=0.5, rnmax=20.0, Nmax=15, R1=1.0, RNmax=100.0), complex_scaling_angle=10.0) for a larger basis set and non-zero rotation angle for complex scaling method

make_phys_params3B1D(; hbar=1.0, masses=[1.0,1.0,1.0], species=[:x,:y,:z], interactions=[[GaussianPotential(-1.0, 1.0)],[GaussianPotential(-1.0, 1.0)],[GaussianPotential(-1.0, 1.0)]], parity=+1)

Create and return a named tuple containing the physical parameters for a three-body system in 1D.

Keyword arguments

• hbar::Float64 = 1.0: Reduced Planck constant used in calculations.
• masses::Vector{Float64} = [1.0,1.0,1.0]: 3-element vector containing the masses of the three particles.
• species::Vector = [:x,:y,:z]: 3-element vector containing particle labels used for automatic (anti-)symmetrization. Use Symbols :b, (:f) for identical bosons or fermions
• interactions::Vector{Vector{Any}} = [[GaussianPotential(-1.0, 1.0)],[GaussianPotential(-1.0, 1.0)],[GaussianPotential(-1.0, 1.0)]]: Vector of Vector of interaction potentials for each pair of particles: [[v231,v232,...],[v311,v312,...],[v121,v122,...]].
• parity::Int = 1: Parity parity=(-1)^(l+L) of the wave function. Possible values: +1,-1 for positive/negative parity, 0 for parity-violating potentials.

Returns

• NamedTuple: Named tuple with the specified physical parameters.

Example

make_phys_params3B1D() # default: system of three different particles with the same mass and Gaussian interactions
make_phys_params3B1D(masses=[1.0,10.0,20.0], species=[:i,:j,:k], interactions=[[v23],[v31],[v12_1,v12_2]]) # system of three different particles with interactions v23 (between particles 2 and 3), v31 (between particles 3 and 1), and v12_1, v12_2 (between particles 1 and 2). The interactions need to be defined before this call.

reduce_basis(S; threshold=1e-13)

Transformation to an orthonormalized, possibly reduced basis: returns the matrix L (size n x m, m <= n) with L' * S * L == I.

L is built from the eigenvectors of the norm-overlap matrix S, dropping all eigenvalues smaller than threshold times the largest one. This removes the near-linear dependence of the non-orthogonal basis; m < n whenever some eigenvalues are cut. Both real-symmetric S and hermitian S (complex-ranged basis functions) are supported.

Used by eigen2step and inverse_solve; useful on its own to transform a generalized eigenvalue problem H*x = E*S*x into the ordinary one (L'*H*L)*y = E*y, with x = L*y.

ISGL_matrices(phys_params, num_params; complex_scaling=false, complex_ranged=:none)

Builds and returns the matrices of a 3D three-body problem without solving it: (; T, V, S), where T is the kinetic energy (not the Hamiltonian T+V), V the interaction, and S the norm-overlap of the basis functions.

Arguments and keywords are the same as for ISGL_solve. Intended for problems which are more naturally posed in terms of the matrices than of the energies, in particular the inverse problem, see inverse_solve.

With complex_scaling=true, T already carries the factor exp(-2*i*theta) and V the rotated potential.

Example

T,V,S = ISGL_matrices(phys_params, num_params)

ISGL_solve(phys_params, num_params; return_wavefunctions=false, complex_scaling=false, complex_ranged=:none, observ_params=(;stateindices=[],centobs_arr=[[],[],[]],R2_arr=[0,0,0]))

Solves the 3D three-body problem using the Gaussian Expansion Method (GEM).

Arguments

• phys_params: Physical parameters for the three-body system (e.g., masses, interaction potentials, etc.).
• num_params: Numerical parameters for the GEM calculation (e.g., basis size, grid parameters, etc.).
• return_wavefunctions: (optional) If true, also returns wavefunction-related observables. Default is false.
• complex_scaling: (optional) If true, uses complex scaling method. Default is false.
• complex_ranged: (optional) Selects complex-ranged basis functions nu -> nu*(1 + i*omega) per Jacobi coordinate, with omega = complex_range_freq from num_params. One of :none (default), :r (only the internal coordinate $r$), :R (only $R$), or :both. A Bool is also accepted, with true meaning :both, for consistency with GEM2B_solve. Each complex-ranged coordinate doubles its number of basis functions, since the ranges enter as conjugate pairs; choosing :both therefore multiplies the total basis size by four. All interaction types are supported, including generic central potentials: their radial integrals are interpolated over the sector of the complex plane in which the effective Gaussian range lives, with kmax_theta angular nodes (see make_num_params3B3D). Observables are not supported with complex ranges.
• observ_params: (optional) Parameters for observable calculations.
  • stateindices: Indices of states for which observables are calculated.
  • centobs_arr: Array of central (only dependent on $r$; must be defined as functions) observables, for each Jacobi set (similar to interactions in phys_params).
  • R2_arr: Array which indicates whether the observable $\langle R^2 \rangle$ should be calculated (1) for any of the three Jacobi sets, or not (0).

Returns

• If return_wavefunctions=false: Returns an array of computed energies.
• If return_wavefunctions=true: Returns a tuple (energies, wavefunctions, centobs_output, R2_output).
  • energies: Vector of computed energies.
  • wavefunctions: Matrix of eigenvectors (column-wise) which contain the coefficients of the basis functions.
  • centobs_output: Mean values of central observables for the specified states. The first dimension corresponds to the Jacobi sets, the second to the observables, and the third to the states.
  • R2_output: Mean squared radii for the R-coordinate. Rows corresond to the Jacobi sets, columns to the states.

Example

phys_params = make_phys_params3B3D()
num_params = make_num_params3B3D()
energies = ISGL_solve(phys_params, num_params) #solving with default parameters: three particles with the same mass and gaussian interaction

csmgaussopt(gaussopt, complex_scaling, complex_scaling_angle)

If complex_scaling=true, returns a new Vector{Vector{Tuple{Float64,ComplexF64}}} where each mu has been scaled by exp(2im * complex_scaling_angle * pi/180). Otherwise returns the original gaussopt unchanged.

eigen2step(e_arr, H, S; threshold=1e-13)

Solves the generalized eigenvalue problem H*x = E*S*x and writes the energies into e_arr.

Similar to eigvals(H,S), but the basis is first orthonormalized by reduce_basis, which cuts the small eigenvalues of S. The number of energies can therefore be smaller than the size of H.

eigen2step_valvec(e_arr, v_arr, H, S; threshold=1e-13)

Same as eigen2step, but also writes the coefficient vectors into the columns of v_arr, normalized as x' * S * x == 1.

inverse_solve(T, V, S; target_energy=0.0, threshold=1e-13, return_vectors=false)

Solves the inverse problem: instead of the energies at a given interaction strength, it returns the interaction strengths lambda at which a state of energy target_energy exists.

For a Hamiltonian which is linear in a strength parameter, H = T + lambda*V, the condition that target_energy is an eigenvalue of H is the generalized eigenvalue problem

(T - target_energy*S) x = -lambda * V x

in the non-orthogonal basis with norm-overlap S. One matrix build and one eigen solve therefore replace a whole scan over the strength with bracketing of the energy.

Arguments

• T: kinetic energy matrix (not the Hamiltonian T+V), as returned by GEM2B_matrices, GEM3B1D_matrices, ISGL_matrices.
• V: interaction matrix at unit strength. lambda is the factor multiplying all interactions, so pass the potential at unit strength.
• S: norm-overlap matrix.

Keywords

• target_energy=0.0: the energy that should be reached. Must lie below the lowest continuum threshold of the basis, such that T - target_energy*S is positive definite; otherwise an error is raised. target_energy=0.0 gives the critical strengths at which a state becomes bound.
• threshold=1e-13: cut-off for the eigenvalues of S, see reduce_basis. Results near target_energy=0.0 are sensitive to this value: too small a cut-off can produce spurious, nearly linearly dependent solutions, too large a one cuts the diffuse basis functions needed close to threshold.
• return_vectors=false: whether to also return the coefficient vectors of the corresponding states.

Real-symmetric and hermitian matrices are supported, i.e. both the usual and the complex-ranged basis functions. The complex scaling method is not: it makes T and V complex-symmetric rather than hermitian, and is rejected.

Returns

• lambdas: the positive strengths, in ascending order. The first entry is the smallest strength at which the system supports a state of energy target_energy. Negative strengths are discarded; if they are of interest, call the function with -V instead.
• vectors: (optional) the corresponding coefficient vectors, normalized as x' * S * x == 1.

Example

T,V,S = GEM2B_matrices(phys_params, num_params)
lambdas = inverse_solve(T,V,S; target_energy=0.0, threshold=1e-8) # critical strengths for binding

make_num_params3B3D(;lmin=0, Lmin=0, lmax=0, Lmax=0, gem_params=(nmax=10, r1=0.2, rnmax=10.0, Nmax=10, R1=0.2, RNmax=20.0), complex_scaling_angle=0.0, complex_range_freq=0.9, mu0=0.08, c_shoulder=1.6, kmax_interpol=1000, kmax_theta=25, threshold=10^-8)

Create and return a named tuple containing the numerical parameters for a three-body GEM calculation in 3D.

Keyword arguments

• lmin::Int = 0: Minimum power r^l used in the basis functions of the $r$ Jacobi coordinate.
• Lmin::Int = 0: Minimum power r^L used in the basis functions of the $R$ Jacobi coordinate.
• lmax::Int = 0: Maximum power r^l used in the basis functions of the $r$ Jacobi coordinate.
• Lmax::Int = 0: Maximum power r^L used in the basis functions of the $R$ Jacobi coordinate.
• gem_params::NamedTuple = (nmax=5, r1=1.0, rnmax=10.0, Nmax=5, R1=1.0, RNmax=10.0): Parameters for the Gaussian Expansion Method (number of basis functions, smallest and largest range parameters for both Jacobi coordinates).
• complex_scaling_angle::Float64 = 0.0: Complex scaling angle (in degrees) for the Complex Scaling Method.
• complex_range_freq::Float64 = 0.9: Parameter $\omega$ controlling the frequency for complex-ranged basis functions, $\nu \to \nu(1 + i\omega)$. Only used when ISGL_solve is called with complex_ranged set to something other than :none. The same $\omega$ is used for both Jacobi coordinates.
• mu0::Float64 = 0.08: Parameter (prefactor) for the ISGL method.
• c_shoulder::Float64 = 1.6: Parameter (base) for the ISGL method.
• kmax_interpol::Int = 1000: Number of numerical integration with effective Gaussian ranges used for interpolation.
• kmax_theta::Int = 25: Number of angular nodes of the range-interpolation. Only used with complex-ranged basis functions, where the effective Gaussian range becomes complex and the interpolation runs over the sector $|\arg\alpha| \le \arctan\omega$ instead of over a real interval. The setup cost of the interpolation grows in proportion, so kmax_interpol can be reduced in compensation. Increase it to check convergence; complex_ranged=:both is the most demanding case, since the four-fold enlarged basis makes the overlap matrix nearly singular and amplifies any error in the matrix elements.
• threshold::Float64 = 1e-8: Numerical threshold for the generalized eigenvalue solver.

Returns

• NamedTuple: Named tuple with the specified numerical parameters.

Example

make_num_params3B3D() # for the default basis set
make_num_params3B3D(gem_params=(nmax=10, r1=0.5, rnmax=20.0, Nmax=15, R1=1.0, RNmax=100.0), kmax_interpol=5000) for a larger basis set, and a finer interpolation grid for matrix-element calculations of interactions with many features.

make_phys_params3B3D(; hbar=1.0, masses=[1.0,1.0,1.0], species=[:x,:y,:z], interactions=[[GaussianPotential(-1.0, 1.0)],[GaussianPotential(-1.0, 1.0)],[GaussianPotential(-1.0, 1.0)]], J_tot=0, parity=+1, spins=[0,0,0])

Create and return a named tuple containing the physical parameters for a three-body system in 3D.

Keyword arguments

• hbar::Float64 = 1.0: Reduced Planck constant used in calculations.
• masses::Vector{Float64} = [1.0,1.0,1.0]: 3-element vector containing the masses of the three particles.
• species::Vector = [:x,:y,:z]: 3-element vector containing particle labels used for automatic (anti-)symmetrization. Use Symbols (e.g. :b, :f).
• interactions::Vector{Vector{Any}} = [[GaussianPotential(-1.0, 1.0)],[GaussianPotential(-1.0, 1.0)],[GaussianPotential(-1.0, 1.0)]]: Vector of Vector of interaction potentials for each pair of particles: [[v231,v232,...],[v311,v312,...],[v121,v122,...]].
• J_tot::Int = 0: Total angular momentum of the system. Can take on half-integer values when spins are involved.
• parity::Int = 1: Parity parity=(-1)^(l+L) of the wave function. Possible values: +1,-1 for positive/negative parity, 0 for parity-violating potentials.
• spins::Vector{Float64} = [0,0,0]: 3-element vector containing the spins [z1,z2,z_3] of the three particles. Used for automatic (anti-)symmetrization.

Returns

• NamedTuple: Named tuple with the specified physical parameters.

Example

make_phys_params3B3D() # default: system of three different particles with the same mass and Gaussian interactions
make_phys_params3B3D(masses=[10.0,10.0,20.0], species=[:i,:j,:k], interactions=[[v23],[v31],[v12_1,v12_2]], J_tot=3/2, spins=[1/2,1/2,1/2],parity=-1, lmax=1, Lmax=1) # system of three different particles with half-integer spin, and interactions v23 (between particles 2 and 3), v31 (between particles 3 and 1), and v12_1, v12_2 (between particles 1 and 2). The interactions need to be defined before this call.

reduce_basis(S; threshold=1e-13)

Transformation to an orthonormalized, possibly reduced basis: returns the matrix L (size n x m, m <= n) with L' * S * L == I.

L is built from the eigenvectors of the norm-overlap matrix S, dropping all eigenvalues smaller than threshold times the largest one. This removes the near-linear dependence of the non-orthogonal basis; m < n whenever some eigenvalues are cut. Both real-symmetric S and hermitian S (complex-ranged basis functions) are supported.

Used by eigen2step and inverse_solve; useful on its own to transform a generalized eigenvalue problem H*x = E*S*x into the ordinary one (L'*H*L)*y = E*y, with x = L*y.