Exactly solvable model of two trapped quantum particles interacting via finite-range soft-core interactions

The exactly solvable model of two indistinguishable quantum particles (bosons or fermions) confined in a one-dimensional harmonic trap and interacting via finite-range soft-core interaction is presented and many properties of the system are examined. Particularly, it is shown that independently on the potential range, in the strong interaction limit bosonic and fermionic solutions become degenerate. For sufficiently large ranges a specific crystallization appears in the system. The results are compared to predictions of the celebrated Busch et al. model and those obtained in the Tonks-Girardeau limit. The assumed inter-particle potential is very similar to the potential between ultra-cold dressed Rydberg atoms. Therefore, the model can be examined experimentally.


The Model
In the following we study properties of the system of two identical quantum particles of mass m confined in a one-dimensional harmonic trap of frequency Ω and interacting via soft-core finite-range rectangular potential. The Hamiltonian of the system reads , if 0, if (2) i.e., the interaction energy is constant and it is non-zero only when the distance between particles is not larger than a. Similar model was considered previously in three dimensions 34 . Depending on the situation we consider symmetric (for bosons) or antisymmetric (for fermions) wave functions with respect to an exchange of particles' positions It is quite natural that in the limit of a → 0 with constrain 2aV = const one restores the Busch et al. model with delta-like contact interaction 1 , while in the limit V → ∞ an extensively studied model of hard spheres is obtained 14 .
Our aim is to give a straightforward and analytical prescription for the eigenstates of the Hamiltonian (1) as a function of the potential depth V and its range a. With these solutions, we examine different properties of a few of the lowest eigenstates in the bosonic and fermionic cases. In particular, we consider different single-particle system characteristics (density profile, momentum distribution) as well as inter-particle correlations reflected in a reduced single-particle density matrix.

The Eigenproblem
To find eigenstates and corresponding eigenenergies of the Hamiltonian (1) it is very convenient to perform standard transformation to the coordinates of the center-of-mass frame: In these new variables the Hamiltonian (1) can be written as a sum of two independent single-particle Hamiltonians  =  R +  r : The Hamiltonian  R has a textbook form of the harmonic oscillator and it can be diagonalized straightforwardly. The Hamiltonian  r has an additional term related to the interactions (2) and the corresponding eigenequation, when written in the natural units of an external harmonic oscillator, has the form: Since the Hamiltonian  r commutes with the operator of the parity inversion, r → −r, the eigenstates of  r can be chosen as either even or odd functions of the relative position r. They directly correspond to bosonic and fermionic statistics, respectively.
The eigenequation (7) has the form of the Weber differential equation 35 : with u equal to −E + V − 1/2 and −E − 1/2 for |r| < a and |r| ≥ a, respectively. The Weber equation was originally studied to solve Laplace equation expressed in parabolic coordinates 35 but it appears in different problems of mathematical physics and many of its properties are well known 36,37 . In the case studied, when the problem is not reduced to the ordinary harmonic oscillator problem (V ≠ 0), it is very convenient to consider two different pairs of the solutions ϕ ϕ having appropriate symmetry under parity transformation  but different properties on the boundaries. The first pair can be expressed in terms of the confluent hypergeometric function 1 F 1 as: These functions may be consider as appropriate solutions only in the region |r| < a because they are divergent in the infinity, r → ±∞. The second pair of solutions is expressed in terms of other confluent hypergeometric function U: In contrast to the first pair, these functions decay appropriately in the limit r → ±∞ but they do not have appropriate behavior at r = 0, i.e., odd functions φ − r ( ) u ( ) are discontinuous, whereas even ones φ + r ( ) u ( ) have discontinuous first derivative. It means that functions φ ± r ( ) u ( ) may be considered as appropriate solutions only in the region |r| ≥ a. Consequently, any solution of the eigenequation (7) with energy E i (which is continuous and has a continuous first derivative in a whole space) may be constructed as following is a normalization coefficient of the resulting function. An additional coefficient ν ± A ( ) together with the eigenenergy E i are determined by matching conditions at |r| = a, ensuring that the function (13) and their first derivatives are continuous in the whole space: As typical for such problems, the conditions (14) and (15) can be fulfilled only for some particular values of an eigenenergy E i leading directly to the quantization of the physical spectrum. In the case studied, the matching conditions (14) and (15) are fulfilled when eigenenergy E i is a solution of the following transcendental equation and μ ν μ determining even and odd solutions, respectively. Equations (16) and (17) are quite complicated but they can be solved straightforwardly with simple numerical methods. After determining eigenenergies one finds the corresponding coefficients and thus corresponding wave functions of the relative motion (13).
It is worth mentioning that in the limiting case of noninteracting particles (V = 0) standard solutions of a one-dimensional harmonic oscillator are restored. Indeed, in this particular case one finds μ = ν = −E i − 1/2 and the matching conditions (14) and (15) reduce to a simple demanding that eigenenergies are half-integer numbers, E i = i + 1/2. In consequence, appropriate functions ϕ ν ± r ( ) ( ) and φ μ ± r ( ) ( ) become equivalent and they are expressed in terms of Hermite polynomials H i : In the opposite limit of infinite repulsions (V → ∞) situation is also simplified. In this case the relative wave functions (13) must vanish in a whole range |r| < a, i.e., all amplitudes = ν ± A 0 ( ) . It immediately leads to the simplified quantization condition and in consequence to the typical for hard-core limit degeneracy between neighboring even and odd solutions.
Having analytical solutions (13) one can express any eigenstate of the Hamiltonian (1) as a simple product of two wave functions and Φ j (r) are appropriate eigenstates of the center-of-mass and relative motion Hamiltonians, respectively. Although these two particular coordinates (R and r) are completely decoupled, the wave function (21) cannot be written (for any finite V) as a product (for bosons) or single Slater determinant (for fermions) of wave functions of independent particles. This observation leads directly to non trivial quantum correlations between particles which are discussed in the following.
Spectral properties. Many properties of the system studied can be extracted directly from the spectrum of the relative motion Hamiltonian (6). In Fig. 1 we show several the lowest eigenenergies of  r as functions of the interaction strength V for different potential ranges a. Solid lines correspond to even wave functions (bosons) while dashed lines to odd cases (fermions). As suspected, for V = 0 the spectrum of noninteracting particles is restored, i.e., alternating bosonic (symmetric) and fermionic (antisymmetric) states of the relative motion have equally distributed energies ħ/2, 3ħ/2, 5ħ/2, etc. For non-vanishing interactions V ≠ 0, depending on the potential range a, eigenenergies vary. A rapidity of these changes crucially depends on a sign of interactions-for the attraction is much higher than for the repulsion. It is interesting to note that for any finite range (a ≠ 0) and sufficiently large attractions each eigenstate may have arbitrary large negative energy. Note also that, independently on the potential range a, energies never cross. It means that for any finite a and V any eigenstate of the system is not degenerated. In fact, it is a direct consequence of one-dimensionality of the relative motion Hamiltonian H r 38 .
In the particular limit of strong repulsions, the neighboring states of opposite symmetry become degenerate independently on the interaction range a. This observation is a direct consequence of the Bose-Fermi mapping 14,16 .
It is clearly seen in Fig. 1 that properties of even and odd solutions of the relative coordinate eigenproblem are essentially different when the potential range a becomes smaller than the natural length scale of the problem (in dimensionless units equal to 1). As it is seen, in contrast to bosonic states, energies of fermionic states (dashed lines) become more horizontal, i.e., they are less sensitive to the interaction energy strength. This behavior is a consequence of the fact that odd functions always vanish at r = 0, i.e., along with decreasing a the interaction energy rapidly decreases independently on interaction strength V. In the limit of vanishing a the interaction is completely described in terms of the s-wave scattering which is not present between indistinguishable fermionic particles.

Contact interactions limit.
Mentioned above limiting case of contact forces can be explored more precisely by considering a formal limit in which the interactions  (r) become identical with δ-like potential of the form gδ(r), i.e., in the limit a → 0 with fixed product 2aV = g = const. In this limit, the problem reduces to the celebrated model of two quantum particles interacting via contact forces for which exact analytical solutions are known 1 . To show how the limiting spectrum is restored we fix the potential range a and for given limiting interaction g we calculate a rescaled value of potential strength V = g/(2a). In this way one obtains the spectrum of the relative motion Hamiltonian H r for different ranges a rescaled to the interaction strength g of the Busch et al. model. Results of this procedure adopted to even (bosonic) eigenstates of the relative motion Hamiltonian are presented in Fig. 2. As it is seen, with decreasing potential range a corresponding eigenenergies approach the results for contact interactions (thick black line). For a = 0.15 (red solid line) an agreement is almost perfect in a wide range of interactions. These results are in qualitative agreement with previously obtained finite range corrections obtained within the Green's function approach for higher dimensionality 39 .
At this point, it should be noted that for any finite a, in contrast to the contact interactions limit, all eigenenergies become negative for sufficiently strong attractions (see Fig. 1). Only in the case of contact interactions there exists exactly one bound state of the Hamiltonian (6)-the ground-state of the bosonic system.
Single-particle quantities. The simplest quantities which can be measured experimentally quite easily are related to single-particle properties. All of them are fully captured by the reduced single-particle density matrix which for the model studied has a form: where Ψ(x 1 , x 2 ) is a chosen two-particle state of the system. In the following we focus on the properties of the bosonic and the fermionic ground-states, i.e., according to the notation of Eq. (21) the states Ψ 00 (x 1 , x 2 ) and Ψ 01 (x 1 , x 2 ), respectively. Typically, we are mostly interested not in the whole reduced single-particle density matrix Γ(x,x′) but only in its diagonal part which represents a density profile of particles. Analogously, a diagonal part of its Fourier transform ip x x ( )/ encodes distribution of a single-particle momentum. Properties of these two simple quantities crucially depend on the range of the potential a. These differences are especially manifested in the cases which are beyond applicability of the Busch et al. model. In Fig. 3 and Fig. 4 we show density and momentum distributions for bosonic (red solid line) and fermionic (dashed blue line) ground-states obtained for a few representative potential ranges (a = 1, a = 1.5, and a = 2) and different potential strengths V, including hard-core limit case V → ∞. Let us recall that in the hard-core limit, the bosonic wave functions necessarily satisfy the condition Ψ 00 (x 1 , x 2 ) = 0 on the line x 1 = x 2 , regardless of a since corresponding wave functions of relative motion Φ(r) vanish at r = 0. This observation, usually called fermionization, enables one to map the bosonic ground-state wave function to the fermionic one via the following relation Ψ 00 (x 1 , x 2 ) = |Ψ 01 (x 1 , x 2 )|. In consequence, in the hard-core limit, the bosonic and fermionic ground states share not only the same energy but also have the same spatial density profiles (bottom row in Fig. 3). Note however, that this particular mapping (forced by infinite repulsions) does not necessarily mean that also momentum distributions for bosonic and fermionic ground-states are the same (bottom row in Fig. 4).
The situation becomes essentially different for potential ranges > ∼ a 2. In these cases, not only the density profiles n(x) but also the momentum distributions π(p) of bosonic and fermionic ground-states become identical in the hard-core limit (right bottom plots in Figs 3 and 4). At the same time the density profiles exhibit spatial  separation into two independent peaks indicating localization of particles at opposite sides of the trap. This result is clearly understandable when two-particle density profile ρ(x 1 , x 2 ) = |Ψ(x 1 , x 2 )| 2 is considered. As it is seen in Fig. 5 the probability of finding both particles at the same side of the trap (x 1 , x 2 0 or x 1 , x 2 < 0) vanishes when > ∼ a 2. It is often said that the system enters the crystallization regime where individual particles are spatially separated 40 and behave as distinguishable parties. In consequence, any physical property of the system does not depend on a quantum statistics. However, as explained in the next section, spatial separation does not mean that individual particles can be treated as independent since non-classical correlations between them are still present.
For finite interaction strengths V < ∞ (see Figs 3 and 4) one can observe how the quasi-fermionization is built along with increasing V. It is quite instructive to note that the fermionization regime is reached earlier if the interaction range a is bigger. It is consistent with previous results concerning quasi-degeneracy in the spectrum of the relative motion Hamiltonian (compare to Fig. 1).
Inter-particle correlations. As was mentioned before, the eigenstates (21) are always separable when written with respect to coordinates (4). However, for any finite interaction, it cannot be written as a product with respect to positions of particles x 1 and x 2 . This simple observation means that the eigenstates of the interacting system studied encode nonclassical correlations between particles. These inter-particle correlations are nicely captured by spectral properties of the reduced single-particle density matrix Γ(x, x′). It is known that the matrix Γ(x, x′) can be decomposed to its natural Schmidt orbitals where u i (x) and λ i are eigenvectors and eigenvalues of the reduced density matrix Γ(x, x′), respectively. Coefficients λ i have a direct interpretation of probabilities of finding a single particle in quantum states described by the corresponding orbitals and they are normalized to unity, ∑ i λ i = 1. Let us note that in the case of fermions, due to Pauli exclusion principle (3), all non-zero eigenvalues are doubly degenerated. If both bosons occupy exactly the same orbital u 0 (x) then the reduced density matrix simply projects to the orbital u 0 (x) and particles are trivially correlated. What is less intuitive, correlations between particles are also trivial whenever particles occupy two different orbitals u 0 (x) and u 1 (x) in the way that the two-particle wave function is represented by their Slater determinant (for fermions) or permanent (for bosons). In all these cases, the state is regarded as non-entangled since correlations originate only in the quantum statistics of indistinguishable particles 41,42 . In other cases, additional correlations are forced by mutual interactions and they are directly reflected in increasing number of non-vanishing occupations λ i in the decomposition (25). In the Fig. 6a we present four the largest eigenvalues λ i of the reduced density matrix Γ(x, x′) calculated for bosonic (solid black lines) and fermionic (dashed blue line) ground-states, Ψ 00 (x 1 , x 2 ) and Ψ 01 (x 1 , x 2 ), as functions of potential range a in the hard-core limit V → ∞. As noted above, in the fermionic case eigenvalues are doubly degenerated. In the contact interaction limit (a → 0), fermions, in contrast to bosons, become noninteracting and therefore there is no additional correlation beyond that induced by quantum statistics (λ 0 = λ 1 = 0.5). In opposite limit of large ranges > ∼ a 2, bosonic occupations λ i become doubly degenerated and equal to fermionic ones. Simultaneously, reduced single-particle matrices of bosonic and fermionic ground-states become identical. By a direct inspection of the two-particle state we found that the ground-state of the system can be written as (± sign for bosons and fermions, respectively): where  i (x) and  i (x) are single-particle orbitals localized in left and right side of the trap constructed from the corresponding even and odd single-particle orbitals of the reduced density matrix Γ(x, x′) 43 : The construction is possible, since in this case the appropriate eigenorbitals are degenerated and any of their linear combination remains as an eigenvector of Γ(x, x′). The amplitudes κ i are related directly to the occupations This observation is one of quite spectacular manifestations of the crystallization mentioned before. Note that the values of the dominant occupations λ 0 and λ 1 decreases with potential range a. It means that even in the crystallization regime particles cannot be treated as trivially correlated parties and therefore they cannot be locally described with individual well-defined orbitals.
Non-classical correlations between particles are quite well quantified by the von Neumann entropy. When occupancies λ i are known, it can be calculated straightforwardly as This measure exactly vanishes in the non-entangled product state of bosons (λ 0 = 1) and it is equal to 1 in the non-entangled bosonic and fermionic states (λ 0 = λ 1 = 0.5). Note however that in the bosonic case, an opposite implication does not hold, i.e., the condition S = 1 does not necessarily means that the state in non-entangled 42 . In Fig. 6b (for bosons) and Fig. 6c (for fermions) we show the dependence of the von Neumann entropy on the potential strength V for different ranges a calculated in the ground-state of interacting particles. It is clear that in a considered range of parameters the von Neumann entropy S is a monotonic function of the interaction V and its growth crucially depends on the potential range a. It is also worth noticing that in the fermionic case and vanishing range (a → 0) the von Neumann entropy remains unchanged independently on the interaction strength V. This is a direct consequence of vanishing contact interaction for fermionic species.

Summary
In this paper, we present properties of the exactly solvable model of two interacting particles confined in a harmonic trap. Inter-particle forces are modeled by a square wall controlled by two independent parameters: potential range and its strength. The results enabled us to investigate and discuss different properties of the system in a whole range of parameters between limiting cases of well known Busch et al. and hard-core models. Obtained results suggest that any finite range of the inter-particle forces is directly reflected in simple quantities which can be measured experimentally. The prominent example is the many-body spectrum where, in contrast to the zero-range case, all eigenstates become unbounded from below for attractive forces. We show that also density and momenta distributions maybe strongly affected by the finite range of the potential. All these deviations from the Busch et al. model maybe examined in recent experiments on a few ultra-cold particles.
The model presented belongs to the specific class of quite realistic quantum many-body problems having exact analytical solutions [44][45][46][47] . From this point of view it is not only an interesting academic example. In fact, it may serve as a first building block for constructions of the many-body ground states of larger number of interacting particles. For example, it can be used as an input for the variational Jastrow-like ansatz based on analytical solutions of two-body problems 48 . Figure 6. (a) First four of the largest eigenvalues of the reduced single-particle density matrix Γ(x, x′) calculated for the bosonic (solid black lines) and fermionic (dashed blue lines) ground-state in the hard-core limit V → ∞ as functions of the potential range a. In the limit of contact interactions (a → 0) fermionic groundstate does not manifest any non-trivial correlations. For large ranges ( > ∼ a 2) bosonic eigenvalues become degenerate and equal to appropriate fermionic ones indicating crystallization regime. Note that fermionic eigenvalues are always doubly degenerated. (b,c) The von Neumann entropy S as a function of interaction strength V for different potential ranges a. Note that in the fermionic case its growth is strongly suppressed in the limit of vanishing range (a → 0) as a consequence of vanishing contact interaction in this limit. In all plots the potential range a and strength V are measured in units of Ω ħ m /( ) and ħΩ respectively.