Three-dimensional non-approximate Coulomb interaction between two trapped quantum particles

The two trapped quantum particles interacting problem is generalized to three dimensions, and the exact Coulomb potential is used. The system is solved by expanding the wavefunction in terms of the isotropic harmonic oscillator eigenfunctions and Hydrogen atom eigenfunctions independently, showing that each one results in a prime approximation for different domains of the normalized coupling constant of the relative interactions, suggesting that the combination of the basis is enough to build a well-suited base for the non-approximate problem. The results are compared to previous works that use a model of approximate finite-rage soft-core interaction model of the problem to give insights into the many-body states of strongly correlated ultracold bosons and fermions. We conclude that the proposed three-dimensional approach facilitates the distinction between bosons and fermions while the solutions given by the expansions better define the behavior of the particles for repulsive potentials. In addition, we discuss the substantial differences between our work and the previous approximate model.

In addition, the study of quantum particles trapped in a harmonic oscillator represents a useful approximation of the microsystems trapped in an optical lattice, so a precise understanding of the interactions in this elementary model may help comprehend the physics occurring in more complex and realistic systems.For instance, stationary solutions are helpful to study the non-equilibrium dynamic of two atoms [21][22][23][24] .This is true with all other simplified models previously cited and those before Kościk and Sowiński 7 , whose work is used as the primary reference in this paper.In the spirit of keep developing new simplified models from which we could build more complex physics, Kościk et al. proposed a finite-rage soft-core interaction model of two trapped interacting quantum particles as an approximation of the Coulomb potential and many properties of the system were derived.Similar to the comparison they made with Busch et al. 9 , mainly focusing on the spectral properties of the system, we will aim to enlarge the already vast pool of studies made in this area but focusing not on the quasi-exact solution, but on devising a proper base to numerically construct solutions to the non-approximate interaction between two particles trapped in a harmonic oscillator.
Hence, we present the generalization of this problem by considering the numerical solution to the exact Coulomb interaction potential in the diagonalization of the relative motion Hamiltonian between the two trapped particles.To do so, we use the three-dimensional isotropic harmonic oscillator eigenfunctions to expand the Hamiltonian and solve the eigenproblem, studying the system behavior when considering the combined harmonic and hydrogen-atom systems.Two distinguishable domains are identified and a clear suggestion on how to build the Hilbert space for a more in detail study is deduced: a combination of the harmonic eigenfunctions, which work nicely for the limit of a repulsive interaction; and Hydrogen atom eigenfunctions, which provide a good approximation for an attractive interaction.A proper combination might be enough to build a well-suited base for the exact Coulomb potential problem between these limits.

Model
The first step we take to generalize the work of Kościk and Sowiński 7 is to consider three dimensions, so we have a system described by the following Hamiltonian: where V corresponds to the Coulomb potential.We use boldface to indicate vectors and normal font to refer to their norm.Using the standard transformation to the coordinates of the center of mass: it is possible to separate Eq.(1) as the sum of two independent Hamiltonians H = H R + H r , one corresponding to the center of mass ( R ) and the other to the relative distance between particles ( r ), which have the following form: where µ is the reduced mass.While H R has the known form of the isotropic quantum harmonic oscillator (IHO), H r has an extra interaction given by the Coulomb potential.The main focus of this report is to find a well-suited base to solve the eigenproblem given by Eq. (3) and, therefore, have a better understanding of the exact Coulomb interaction between quantum particles in a trap.To do so, we consider the dimensionless operators: with the momentum operator p = −i ∇ r , to transform Eq. (3) into the dimensionless Hamiltonian where L 2 is the magnitude of the angular momentum operator, such that the spherical harmonic functions are eigenfunctions of L 2 , namely, L 2 Y ℓ,m (θ, φ) = −ℓ(ℓ + 1)Y ℓ,m (θ, φ) .Notice that the model only has one free parameter γ = Z 1 Z 2 e 2 µ/h 3 ω , unlike the finite-range interaction of Ref. 7 , where a strength interaction param- eter V and a range parameter a that describes the distance in which the particles interact with each other are needed.In our case, as we do not approximate the Coulomb interaction, we consider the range parameter a → ∞ and only need γ to help us study the intensity of attraction ( γ < 0 ) or repulsion ( γ > 0 ) between the two particles.Furthermore, using γ = 0 and γ → −∞ recovers the IHO and a Hydrogen-like atom, respectively, which we use to compare later. (1)

Method
As H r may suggest, we approach this problem using two different known solutions as a base: the three-dimen- sional isotropic harmonic oscillator and the Hydrogen-like atom (HA).We take one of these basis functions and expand the "extra" term in Eq. ( 4) to build the complete matrix we later diagonalize.For example, if we use the IHO eigenfunctions, we end up with a diagonal matrix given by the eigenenergies of the Harmonic problem plus a matrix corresponding to the Coulomb term, which we build expanding γ /r on the IHO base.Considering the Hydrogen eigenfunctions, we would have a diagonal matrix given by the Hydrogen energies plus a harmonic matrix expanded into the HA base.Having this in mind, we only need to calculate the extra term in both scenarios to diagonalize and find the system's spectra.
Starting with the known solution to the three-dimensional IHO: where f nℓ is a polynomial with coefficients given by the recurrent relationship This solution already lets us distinguish the symmetric and anti-symmetric spatial wave function with respect to an exchange of particles positions, given by the symmetry of the spherical harmonics: This is a simple distinction: an even value of ℓ corresponds to bosons and an odd value of ℓ to fermions.We then have a natural separation between particle types because of the three-dimensional solution we have chosen.
Continuing with the eigenproblem, we end up with the following equation using ( 5) and ( 4) and since the energy levels of the harmonic oscillator are we only need to calculate the matrix elements of the Coulomb interaction, which can be reduced to expressions involving the Gamma function with β = µω/ .
Conversely, using the radial part of the HA eigenfunctions we end up with a new equation similar to Eq. ( 8), namely, Equations ( 8) and ( 12) summarize the eigenproblem we will focus on since they provide the non-diagonal matrices to test the efficiency and accuracy when using either the IHO or the HA eigenbasis in the solution via expansion of two trapped quantum particles interacting via a non-approximate Coulomb potential.The spectral properties and the shape of the orbitals obtained using this method will be compared to the numerical solution of the Schrödinger equation with the relative motion Hamiltonian (4) using the finite elements method (FEM), in which we discretize Hψ(r) = E ℓ ψ(r) into a mesh using divided differences with a spatial discretization r to transform it into a matrix eigenvalue problem.Using standard techniques we obtain the eigenstates ψ n,ℓ (r) with eigenenergies E n,ℓ .

Results
For simplicity, we will analyze the case of particles with equal mass, so µ = m/2 .We use energies in dimension- less units and normalize the distance with respect to the Bohr radius ( a 0 ).We obtain an efficient convergence of the energies of a repulsive ( 0 < γ ) system H r using the IHO eigenfunctions with ℓ = 0 when solving the eigen- problem given by Eq. ( 8).In fact, we obtain a good approximation to the ground state with just five elements in the base.This tendency continues as we calculate higher energies: we only need n + 3 to n + 5 elements in the IHO base for the n-th energy value to converge.For example, when going from n = 5 to n = 10 , with ℓ = 0 and γ = 5 , the relative error in the 5-th eigenvalue decreases by 38%, from 17.4 to 12.6.
We then contrast these results using second-order perturbation theory to check that the approximation converges to the proper result.A further check is made by rearranging the Hamiltonian to obtain the known Hydrogen atom energy levels as follows: we have www.nature.com/scientificreports/if we consider γ = −1 such that V 1 = −1/r , then the HA Hamiltonian is with V 2 as the harmonic potential.Solving Eq. ( 13) with the same expansion used for the IHO eigenfunctions results in a suitable approximation for the known Hydrogen atom energies E HA n = − 1 2n 2 in atomic units, where different values for the amplitude of oscillation gives us additional precision.Therefore we confirm that the IHO eigenfunctions build a valuable basis for this system.
On the other hand, calculating the matrix elements with HA eigenfunctions to solve the diagonalization problem with the expansion of the harmonic potential, produces, as expected, a poor result as HA eigenfunctions do not form a complete basis if we do not consider the scattering states: Hydrogen atom eigenfunctions on their own, despite converging, does not work and results in wrong values for the system energy levels.
A further study of the behavior of the system is made by analyzing how the energy levels change by the coupling factor γ in Eq. ( 4) and comparing them to those obtained by Kościk and Sowiński, and the direct solution to the differential equation via the finite elements method, which can be seen in Fig. 1.The finite element solution was calculated using divided differences up to 2nd order for the derivatives, using r = 0.01 .For the calculation, we write the eigenvalue equation in terms of finite elements up to r max = 10 in normalized units with a r reso- lution.In Fig. 2, we note that this is a reasonable value of r max from the trapped system.In addition, we studied the eigenvalue convergence for other values of r max and reached the conclusion that r max = 10 was a reasonable compromise between calculation time and accuracy.We also checked with r = 0.001 , and we obtained the same result as in Fig. 1.Additionally, the use of isotropic harmonic oscillator eigenfunctions with ℓ = 0 results in an excellent approximation for the repulsive potential 0 < γ while failing for γ ≪ 0 , which represents a highly attractive potential between the particles.For a broader analysis of this behavior, we contrast the radial probability density distribution of the IHO, HA, the system solved using 18 IHO eigenfunctions and finite elements method for different values of γ , as shown in Fig. 2. We can see that, opposite to the IHO, the HA base reasonably describes the attractive system γ < 0 while failing to approximate the repulsive case 0 < γ .For this last case, we used γ = −0.1 and γ = −1 to calculate the Hydrogen-like density distribution while using repulsive potentials for the other three curves.An analogous analysis was made by plotting a set of three-dimensional orbitals of the Hydrogen atom and comparing them to those obtained by solving the relative motion Hamiltonian (4), which are shown in Table 1.Despite the differences in scale, which are informed by the expectation value of the radius, it can be seen that the solution of the relative motion Hamiltonian for a highly attractive potential results in a shape identical to that of the Hydrogen atom.www.nature.com/scientificreports/Finally, as we mention in Eq. ( 7), the distinction between bosons and fermions is reduced to the parity of ℓ , so using ℓ = 1 helps us describe the interaction between two fermions, whose contrast with a bosonic (n, ℓ) = (0, 0) system can be seen in Fig. 3. Despite getting closer for positive values of γ , the two systems are not equal as previous works 7 have suggested.Moreover, the exact solution using FEM shows that there is also degeneracy between neighboring bosonic and fermionic energy levels for the highly attractive potential.Another meaningful difference between symmetric and antisymmetric particles is that the fermion eigenfunction converges faster to the expected result due to the Pauli exclusion principle.

Conclusions and Discussion
We focused on the non-approximate Coulomb interaction between two trapped particles using the isotropic harmonic oscillator and Hydrogen atom eigenfunctions independently to expand the relative motion Hamiltonian between two trapped quantum particles and numerically solve the eigenproblem by diagonalization, showing that each of these functions has complementary domains in which they result in a good approximation to the system.The IHO eigenfunctions are an ideal base for the repulsive potential but also work for weak attraction, resulting in adequate approximations for values up to γ ≈ −5 , as shown in Fig. 1.In contrast, the HA eigen- functions work best for highly attractive potentials.This behavior was visualized by studying the radial density distribution of our model and comparing it to the isotropic harmonic oscillator and the Hydrogen atom, as well as the numerical solution obtained by the finite elements method.The system's inclination towards the two different domains is shown in Fig. 2, where it is clear that the IHO basis fails to approximate the system when it approaches the Hydrogen-like atom.Furthermore, the similarity to the Hydrogen atom for highly attractive potentials γ ≪ 0 is checked in Table 1 at the γ = −10 column; despite the differences in scale due to the har- monic trap, it is identical to the Hydrogen atom orbitals.This suggests that a combination of these functions may be enough to build a base that helps us solve the non-approximate Coulomb interaction between trapped Figure 2. Radial probability density distribution of the relative motion in bosonic ground state (n, ℓ) = (0, 0) for different values of the attraction coefficient γ , where a 0 is the Bohr radius.As a repulsive potential does not allow the Hydrogen-like eigenfunctions to be localized, we use attractive potentials when γ = 0.1 and γ = 1 .It can be seen that our solution (blue) using 18 elements of the IHO base best approximates the numerical solution (dotted) for 0 < γ , where the system approaches the IHO (soft orange).For γ ≪ 0 , where the system resembles the hydrogen-like (red) behavior, our method fails to model the system, just like Fig. 1  www.nature.com/scientificreports/particles using the expansion in eigenfunctions.The wide domain in which IHO functions work and the low number of elements that are needed in the base for the approximation tells us that maybe a low number of HA eigenfunctions are required to be added for the solution to be generalized to positive and negative values of γ , thus solving the differences observed in Figs. 1 and 2 between our results and the numerical result using FEM.
As we think that the expansion in the IHO basis may be useful for further studies, we include the coefficients of the IHO expansion for a number of γ values as supplementary data.
Our model, which has just one variable coefficient to explore the attraction and repulsion between trapped particles, captures qualitatively the results from previous works and deepens the understanding of what features are specific to the election of a particular potential.Thus it expands on the problem that started with Busch et al.In addition, our analysis clearly distinguishes between the behavior of bosons and fermions, a characterization that comes naturally with a three-dimensional approach.With the exact potential, substantial differences are observed in degeneracy for bosonic and fermionic energy levels: degeneracy is observed between neighboring eigenenergies for highly attractive systems, as Fig. 3 shows; and degeneracy for coefficients 0 < γ more accurately approach the hard-core limit.This makes us confident that with few additions, a robust model can be built to later study different types of interaction and other numbers of particles, thus adding more theoretical understanding to an area at the forefront of experimental studies.Figure 3. Spectrum of the relative motion Hamiltonian H r as a function of the coupling factor γ for bosonic (n, ℓ) = (0, 0) (teal) and fermionic (n, ℓ) = (0, 1) (red) ground state obtained using 18 elements in the expansion in IHO eigenfunctions (solid).These results are compared to the FEM solution (triangles) using r = 0.01 .There is a more accurate behavior between neighboring eigenenergies for the repulsive non-approximated interaction, as the degeneracy between them should be obtained in the hard-core limit γ → ∞.Table 1.Three-dimensional visualization of the orbitals of the Hydrogen atom compared to those obtained by solving the relative motion Hamiltonian (4) for attractive and repulsive coupling constants γ (columns) and different atomic numbers (rows).As Eq. ( 7) states, odd (even) values of ℓ correspond to fermionic (bosonic) particles.The expectation value for radius is displayed under every figure to show the difference in scale.This is mainly due to the effect of the harmonic trap but also because of the attractive potential in the case of γ = −10, −1 .Despite the similarities at first glance, details in the geometry and size of the orbitals change, where it can be seen that the highly attractive system approaches the Hydrogen atom orbital best.

Figure 1 .
Figure 1.Spectrum of the normalized Hamiltonian (4) as a function of the coupling constant γ using the first 6 (red) and 15 (teal) eigenstates of the IHO base to expand and solve the eigenproblem.These results are compared to the direct solution of the differential equation via the finite elements method (dotted) with r = 0.01 , showing that, as the number of elements used to solve the diagonalization problem increases, the energy gets closer to those obtained by the finite elements method.