Exactly solvable model of two interacting Rydberg-dressed atoms confined in a two-dimensional harmonic trap

Exactly solvable model of two Rydberg-dressed atoms moving in a quasi-two-dimensional harmonic trap is introduced and its properties are investigated. Depending on the strength of inter-particle interactions and the critical range of the potential, the two-particle eigenstates are classified with respect to the excitations of the center-of-mass motion, relative angular momentum, and relative distance variable. Having these solutions in hand, we discuss inter-particle correlations as functions of interaction parameters. We also present a straightforward prescription of how to generalize obtained solutions to higher dimensions.

exactly solvable model of two Rydberg-dressed atoms moving in a quasi-two-dimensional harmonic trap is introduced and its properties are investigated. Depending on the strength of inter-particle interactions and the critical range of the potential, the two-particle eigenstates are classified with respect to the excitations of the center-of-mass motion, relative angular momentum, and relative distance variable. Having these solutions in hand, we discuss inter-particle correlations as functions of interaction parameters. We also present a straightforward prescription of how to generalize obtained solutions to higher dimensions.
Few-body systems of ultra-cold atoms provide a very comprehensive toolbox for exploring fundamental properties of quantum systems containing a mesoscopic number of particles [1][2][3] . Due to accessible tunability of their different parameters they can serve as quantum simulators of strongly correlated quantum systems of a few particles described by models being far beyond computational facilities of computers nowadays 4 . Typically, in the context of ultra-cold physics one assumes that mutual interactions between atoms are dominated by short-range forces and may be represented by simple s-wave (for bosons) or p-wave (for fermions) scattering processes. However, mainly due to the experimental progress with polar atoms and molecules, also dipolar long-range and anisotropic interactions are widely considered and they were found to have very interesting consequences for the system's properties [5][6][7] . Alternatively, long-range interactions between ultra-cold atoms can be achieved when their excitations to the Rydberg states are considered, i.e., when atoms become excited to large principal numbers 8 . In such a case, mutual interactions are not only long-range and strong but also have multipolar properties. In fact, the resulting interaction potential can be viewed as a combination of the standard long-range van der Waals interaction acting over large distances (measured recently with a very sensitive experimental scheme 9 ) and soft-core, almost constant, potential when atoms are close enough [10][11][12][13][14] . Since coherent excitations to Rydberg states on demand were recently announced by many experimental groups (for review see 15 ) Rydberg atoms become one of the candidates for fundamental blocks of future quantum simulations. In such cases, a deep understanding of their spatial correlations may be fundamentally important.
Inspired with the above experimental motivations, in this work we study the problem of two interacting Rydberg atoms confined in a two-dimensional isotropic parabolic trap. The discussion is carried out in the framework of the simplified but exactly solvable model. The simplification is based on the assumption that the main contributions to the spatial properties of the system come from the soft-core part of interactions. This part we model simply by a flat potential of a finite range and strength. At the same time, we assume that the long-range part of the interactions is adequately less important and can be safely omitted. In this way, we end up with the interaction potential modeled by a step function in the relative distance between particles. It turns out that assuming this simplified shape of inter-particle interactions one can fully solve the corresponding two-particle Schrödinger equation in terms of special functions. Along with the discussion, we argue that the simplified approach is justified provided that the critical radius is not very large when compared to the natural length of the trapping potential.
Our work is organized as follows. In Sec. 2 we introduce the simplified model of interacting Rydberg-dressed atoms confined in a two-dimensional parabolic trap and we explain its origin. In Sec. 3 we present a full solution of the corresponding two-particle eigenproblem in terms of the hypergeometric confluent functions. Then in Sec. 4 we perform classification and discussion of the spectrum of the system and its lowest eigenstates. At this point, we also validate our simplified model by comparing its predictions with numerical results obtained for the same system but with interactions modeled by the realistic potential. In Sec. 5 we focus on inter-particle correlations for different parameters of interactions, while in Sec. 6 we shortly explain how one can generalize our analytical solutions to the case of two Rydberg-dressed atoms in a three-dimensional harmonic trap. Finally, we conclude in Sec. 7.

the Model
In this work we consider the system of two interacting quantum particles of mass m confined in an external quasi-two-dimensional harmonic isotropic trap of frequency Ω. The Hamiltonian of the system reads where r i = (x i , y i ) are positions of particles. We assume that the dimensional reduction of the problem to two spatial dimensions is granted by very deep confinement in the remaining third spatial dimension. In such a case any spatial excitations in this direction are strongly suppressed. Consequently, one can safely assume that the dynamics is frozen and particles occupy only the lowest single-particle orbital. It is well-known that in the case of highly-excited ultra-cold off-resonantly dressed Rydberg atoms the effective interaction potential have a very characteristic form [10][11][12][13][14]16 . It can be viewed as the natural van der Waals interaction acting on large distances (~1/r 6 ) with significant modification when the inter-particle distance is comparable with so-called critical Rydberg radius R c . It is argued that the potential can be written as: Here, the two independent parameters g and R c describe characteristic scales of the potential and they are related to the interaction strength and aforementioned critical distance at which interactions change their character.
Since these parameters directly depend on experimentally accessible quantities, namely the effective Rabi frequency ω R and the detuning Δ, as ω = Δ g 2 / R 4 3 and R c = (C 6 /2Δ) (1/6) (C 6 is the dispersion coefficient being constant for chosen atom), they can be treated as parameters which can be tuned on demand (for details see for example 14 ).
In fact, the potential (2) belongs to the large class of potentials of the form 1 having an almost flat soft-core in the center ( ⪅ r R c ) and a long-range tail decaying algebraically as ~r −α . In the limit of very large powers α this class of potentials is exactly equivalent to the simplified potential of the form provided that V = g and a = πR c /[α sin (π/α)]. It means that in the case studied (α = 6) the potential (2) can be modeled by approximate potential (3) by fixing a = πR c /3 (see Fig. 1). Since the power α is quite large the approximation can be treated as reasonable.
Having this argumentation in mind, in our work we model mutual interactions between particles with the finite-range soft-core potential (3). In this approximation, the interaction energy vanishes whenever the distance between particles is larger than the potential range a and it has non-vanishing constant value g at short distances. In the following, we show how to find all eigenstates of the Hamiltonian (1) with corresponding eigenenergies. In this way, we generalize recent results obtained for infinite interaction strength (V → ∞) 17 , as well as for the corresponding one-dimensional problem 18,19 , and we extend the list of exactly (or almost exactly) solvable models Figure 1. Schematic comparison of the realistic shape of interaction potential V R (r) between Rydberg-dressed atoms (2) and the simplified rectangular box potential V(r) given by (3). In the latter case, the long-range part is completely neglected while the short-range part is replaced by a constant energy shift. Substantial differences between the two models are visible only in the vicinity of the critical radius R c . eigenproblem The first step towards diagonalization of the Hamiltonian (1) is to separate the center-of-mass motion. Since particles are confined in a harmonic trap the separation is done by performing the following transformation of coordinates to the center-of-mass and relative motion positions Indeed, in this coordinates the Hamiltonian separates into two independent parts R    = + ξˆˆ of the form: where M = 2m and μ = m/2. Consequently, every eigenstate of (1) can be written as Υ(r 1 , r 2 ) = Φ(R)Ψ(ξ). The center-of-mass Hamiltonian (5a) describes the single particle of mass M confined in a two-dimensional harmonic trap and its eigenfunctions Φ NL (R) with corresponding eigenenergies can be found straightforwardly. They are enumerated with two quantum numbers (N, L) related to excitations in the radial direction and the angular momentum of the center of mass, respectively. The relative motion Hamiltonian (5b) is equivalent to the Hamiltonian of a single particle of mass μ confined in a two-dimensional harmonic trap imposed in the center by the additional rectangular potential V( ) |ξ| . Our aim is to analyze all the properties of the Hamiltonian ξ  . For convenience we express all quantities in the natural units of the harmonic oscillator, i.e., energies, positions, and momenta are measured in Ω  , m / Ω  , and m  Ω, respectively. Note that the Hamiltonian ˆξ  commutes with the relative angular momentum operator L i / = − ξ × ∂ ∂ξ. Therefore, to reduce complexity of the problem we rewrite it to the polar coordinates ξ = (ρ, φ) and we represent all its eigen wave functions Ψ(ξ) in the standard angular momentum representation The wave function Ψ(ρ, φ) obeys the Schrödinger equation and the radial part f(ρ) fulfills the one-dimensional radial Schrödinger equation of the form Due to the specific form of the interaction (3), the eigenequation (7) has simplified form in the two disjoint regions ρ < a and ρ > a: It is a matter of fact that for any E the Equation (8) has two independent solutions (for V ≠ 0) which can be expressed in terms of confluent hypergeometric functions U and 1 F 1 as follows: Since the function f 1 (ρ) is divergent in the limit ρ → ∞, any physically acceptable solution of (8) can be constructed only as the composition www.nature.com/scientificreports www.nature.com/scientificreports/ with appropriately chosen energy  and coefficient A to match both parts at ρ = a and to assure that the wave function f(ρ) is continuous and differentiable in a whole space where the quantum number n = 0, 1, … enumerates successive roots of the Equation (12). Finally, the eigenstates of the relative motion Hamiltonian  ξ are determined by two quantum numbers n ( , )  and the angular momentum orientation and they have the form where f ( ) n ρ is given by (10) provided that the eigenenergy   n is the n-th root of the transcendental Equation (12). At this point, it is interesting to point out that for some particular potential parameters V and a the exact solutions may be significantly simplified. It may happen when the first arguments of the confluent hypergeometric functions U and 1 F 1 are integers and the functions are expressed in terms of simple algebraic expressions 31 Table 1. Such simple algebraic solutions play a very important role, since they may serve as benchmarks for the accuracy of different numerical techniques.
Finally, we want to emphasize that up to now the solutions are in fact obtained for distinguishable particles described by the Hamiltonian (1). In the case of indistinguishable atoms (fermions or bosons), one should impose additional requirements to the wave functions of the relative motion under exchange of particles' positions, ξ → −ξ. In consequence, the bosonic (fermionic) states have even (odd) angular momentum quantum numbers . It means that the radial distribution of fermionic relative motion must necessarily vanish at ρ = 0. This fact can  www.nature.com/scientificreports www.nature.com/scientificreports/ be viewed as a direct manifestation of the Pauli exclusion principle forbidding any two identical fermions to occupy the same position.

Eigenstates Classification
It is very instructive to start the analysis of the system's properties focusing on the ground-state energy of the relative motion Hamiltonian (5b) in individual subspaces of given relative angular momentum, i.e., states with n ( , ) (0, ) =  . In Fig. 2 we present resulting spectra for several potential ranges a ∈ {0.5, 0.75, 1.0, 1.25}. For completeness, we compare eigenenergies obtained in our model of interaction potential V (solid black lines) with those obtained numerically when the realistic model of mutual interactions V R is considered (red dots). It is clear that for not too large ranges of the potential, predictions of both approaches agree in a wide range of interaction strengths. Deviations are clearly visible for large potential ranges a 1 ⪆ and/or adequately strong interactions ⪆ V 7. Note, that even for quite large potential range, a = 0.5, and quite large strengths V, only the s-state with  0 = is influenced by interactions. It is clear that states which are the most sensitive to interactions are characterized by the smallest relative angular momentum quantum numbers . Along with increasing , radial distributions of relative motion are pushed out from the center due to the additional centrifugal term in the Hamiltonian. In consequence, they are less sensitive to the interaction core. This effect is clearly seen when the radial density distribution of the relative motion is considered, F(ρ) = f 2 (ρ)/ρ. In Fig. 3 we plot this distribution for the bosonic  = n ( , ) (0, 0) and the fermionic  n ( , ) (0, 1) = ground states of the relative Hamiltonian (5b) and different potential strengths V (here we set a = 1). For increasing repulsions the density probability is suppressed in the (labels on solid lines) for different values of the potential range a. Solid lines are obtained for simplified potential V while red dots for realistic interaction potential  V with R c = 3a/π. It is clearly seen that for ranges smaller than the natural length of the harmonic oscillator (a 1 ⪅ ) both approaches give very similar results in a wide range of interaction strengths. www.nature.com/scientificreports www.nature.com/scientificreports/ center, while for attractions it is enhanced. In fact, this effect almost does not depend on long-range tails of the potential and it is an exclusive consequence of the potential core. It is clear when we compare the distributions obtained in the simplified model of interaction V (solid lines) with predictions of the realistic model V R (dotted lines). In both cases, the resulting distributions are almost identical (see Fig. 3). In the limiting case, a → 0, only the s-states ( 0 =  ) are affected by interactions since only these states have non-vanishing distributions at ρ = 0. The situation becomes even more interesting when excitations of the relative motion are considered. In Fig. 4 we present the energy spectrum of the relative motion Hamiltonian (5b) for two different potential ranges classified accordingly with their quantum numbers  n ( , ). As it is seen, for interactions having larger potential ranges, an energetic order of eigenstates can be changed. Moreover, almost perfect degeneracies between different states visible for smaller ranges are lifted (compare behavior of states {(0, 4), (1, 2)} or {(1, 1), (0, 3)} for a = 0.5 and a = 1.0). Finally, let us draw some attention to the effect of decreasing splitting between the two lowest eigenstates  (2). Note that in the case of fermionic particles, due to the Pauli exclusion principle, the radial density distribution F(ρ) must necessarily vanish at ρ = 0. (0, 0) and (0, 1) (the latter is in fact doubly degenerated due to the orientation of the angular momentum) when larger ranges are considered. Vanishing of this particular gap may have crucial experimental consequences for distinguishable particles since then even very small but a finite temperature of the system may lead to the statistical mixing of these states and significantly change measurable properties of the system. On the other hand, in the situation of a small gap the states can be easily coupled by some additional well-controlled but non-conserving angular momentum interactions (for example spin-orbit coupling).

inter-particle correlations
Having analytical expressions for two-particle eigenstates of the interacting system it is very easy to perform the inverse transformation and obtain two-particle wave functions ( , ; , ) NL n ; 1 1 2 2  ρ ϕ ρ ϕ ϒ expressed by particles' real-space positions r 1 = (ρ 1 , φ 1 ) and r 2 = (ρ 2 , φ 2 ). Then, one can straightforwardly analyze different interesting features of inter-particle correlations. Here we focus on the two simplest quantities which directly encode information about relative spatial correlations between particles. The first is the two-particle radial distribution n(ρ 1 , ρ 2 ) defined as which can be directly interpreted as the probability density that simultaneously observed particles are at distances ρ 1 and ρ 2 from the center of the trap (see Fig. 5a). It is evident that along with increasing interaction strength V particles are pushed out from the center of the trap and their radial distances become correlated, i.e., the probability of finding particles at the same distance from the center increases. The second quantity is the two-particle azimuthal distribution Γ(φ 1 , φ 2 ) defined as  www.nature.com/scientificreports www.nature.com/scientificreports/ It is related to the probability that in a simultaneous measurement of particles' positions the position vectors will be oriented at angles φ 1 and φ 2 , respectively. If the two-particle quantum state is rotationally invariant (for example the two-boson ground-state of the system has this property) then the distribution (15) depends only on a difference γ = φ 1 − φ 2 and then one can consider the simplified distribution Γ(γ) = 2πΓ(γ + φ 0 , φ 0 ) (with arbitrary chosen φ 0 ) encoding probability density for the relative angle between position vectors.
We plot the distribution Γ(γ) for different interaction strengths V and a = 1 in Fig. 5b. Obviously, for V = 0 the distribution is flat and equal to (2π) −1 . When repulsive interactions are switched on, the probability that particles occupy opposed sides of the trap (γ = ±π) is strongly enhanced. Contrary, for attractive interactions (V < 0), particles are more likely to be located on the same side of the trap (γ = 0). One can quantify an uncertainty that particles are exactly on the opposite (same) side by calculating the standard deviation defined as with γ π = and 0 γ = for repulsive and attractive interactions, respectively. We display this quantity in Fig. 5c. It is clear that along with increasing interactions uncertainty decreases. Surprisingly, it decreases also when the range of the potential a increases. It means that larger ranges try to force particles to form a line with the center of the trap. This observation is also supported by results obtained for infinite potential strength V → ∞ (inset in Fig. 5c). As it is seen, in this limit the standard deviation σ Γ (∞) decreases with the potential range a. Interestingly, the effect of the potential range is not so obvious for attractive interactions. In this case, the behavior of the system strongly depends on potential strength. All these three quantities together [n(ρ 1 , ρ 2 ), Γ(φ 1 , φ 2 ), and σ Γ (V)] give a quite nice view on the spatial correlations induced by interactions build in the system. It can be summarized as follows. When repulsions are increased, positions of particles are forced to arrange exactly on opposite sides of the trap in a quite well-established distance from the center which is determined by the potential range. The effect is stronger for larger a. Contrary, when interactions are attractive, the most probable situation is that particles are found on the same side of the trap.

Generalization to three Dimensions
Finally, let us also mention that the presented solutions may be easily generalized to the problem of two-particles confined in the isotropic three-dimensional harmonic trap. In this case, all eigenfunctions of corresponding three-dimensional relative motion Hamiltonian are classified by three quantum numbers n m ( , , )  and have a form are three-dimensional spherical harmonic functions. The radial part of the eigenfunction ρ f ( ) n fulfills the following single-particle Schrödinger equation This eigenproblem is exactly equivalent to the previous two-dimensional problem (7) provided that one perform appropriate substitution 1/2 + in (7). Consequently, by applying this substitution in (10) and (12), one obtains three-dimensional eigenstates and transcendental equation for eigenenergies, respectively.

Summary
To conclude, in our work we introduced a simplified model of two interacting Rydberg-dressed atoms confined in a harmonic trap. The model is a consequence of replacing the realistic shape of interaction potential by the soft-core finite-range forces modeled by a step function. The main advantage of the model proposed is its exact solvability in terms of special functions. This gives a route for analytical analysis of different properties of the system which are crucially important when inter-particle correlation are considered. By performing detailed numerical analysis, we show that the eigenstates obtained in the simplified model are very close to those obtained in a realistic model in a wide range of interaction parameters. Although our work is devoted only to two interacting atoms, the solutions presented can be used as building-blocks for approximate methods dedicated to a larger number of particles. For example, similarly as it was done in different one-32-34 , two- 35,36 , or three-dimensional 37 , one of the possible extensions is to use these two-particle solutions when variational ansatz of pair-correlated Jastrow wave functions are constructed 38 . Typically for this construction of the variational family, the corresponding two-body exact solution may serve as a prescription for trial wave functions.
Finally, we want to point out that in cases of more than two interacting Rydberg atoms, the simplified model of interactions should be used carefully since neglecting the long-range part of interactions may lead to false conclusions. For sure, the simplified soft-core finite-range potential captures the most important part of realistic interactions between Rydberg atoms. As long as long-range tails do not affect (or affect much weaker) a third particle being far from considered pair, the simplified model should appropriately describe properties of the many-body system. It means that the simplified model of interactions gives an appropriate description only when two-body interactions are not substantially affected by long-range interactions with other particles. In other cases, the long-range tails may substantially change properties of the system and introduce additional correlations.