Quantum Entanglement of Two Harmonically Trapped Dipolar Particles

We study systems of two identical dipolar particles confined in quasi one-dimensional harmonic traps. Numerical results for the dependencies of the entanglement on the control parameters of the systems are provided and discussed in detail. In the limit of a strong interaction between the particles, the occupancies and the von Neumann entropies of the bosonic and fermionic ground states are derived in closed analytic forms by applying the harmonic approximation. The strong correlation regimes of the system with the dipolar bosons and the system with the charged ones are compared with each other in regard to aspects of their entanglement.


Introduction
In recent years, systems of quasi one-dimensional (1D) systems of cold atoms with a short-range interaction have attracted considerable research attention. In particular, an observation of the Tonks-Girardeau (TG) systems [1], where bosonic systems behave as gases of spinless non-interacting fermions [2], has inspired a great interest in exploring their properties. The recent technological developments have also opened up perspectives for the experimental realization of systems of spatially confined particles with long-range dipole-dipole interactions (DDI) [3][4][5], providing, among other things, new possibilities for studying quantum correlation effects in many-body systems. Since then, there has been a remarkable increase of interest in understanding the properties of such systems [6][7][8][9][10][11][12][13]. For recent developments in studies of Dipolar Quantum Gases, see the overview [14].
Also, the study of entanglement has attracted much attention within the last few years, as entangled systems play an important role in many different research areas such as quantum information technology [15] and quantum phase transitions [16]. Particularly the research activity has expanded towards investigating the entanglement properties in various systems composed of interacting particles. For instance, the recent studies include model systems such as the Moshinsky atom [17][18][19][20][21][22], helium atoms and helium-like atoms [23][24][25][26], quantum dot systems [27][28][29][30], and 1D systems of atoms interacting via a short-range contact interaction [31][32][33]. For details concerning the recent progress in entanglement studies of quantum composite systems, see [34]. However, according to our best knowledge, studies on the entanglement of dipolar particles have not yet drawn much attention.
In this paper, we address this issue and gain some insight into the quantum entanglement properties of systems composed of two particles confined in a harmonic trap with the DDI modeled by where θ rd is the angle between r and d, cos 2 θ rd = r · d/rd, and d 2 being the strength of the DDI. For the sake of simplicity, we focus here on the regime of strong anisotropy, = ω ⊥ /ω 1. The theoretical description of such quasi-1D systems can be simplified, since one can assume that the particles stay in the lowest transverse confinement mode and the single-mode approximation (SMA) can be applied. Within this approximation, the system of four dipolar bosons in a trap of anisotropy of = 50 has recently been considered in [7], wherein the effects of the interaction strength on various characteristics (such as the density, the momentum distribution, and the occupation number distribution) have been determined. For the two-particle system under consideration here, the Hamiltonian in the SMA is with where we have omitted the short-range contact interaction in order to explore only the pure dipolar effects (for details on this point, see, for example, Deuretzbacher et al. [7]). The coordinates and the energies are measured in terms of √ mω/h,hω, respectively, and the dimensionless coupling g is related to the control parameters of the system by g = d 2 √ ωm 3 2 (1 + 3cos2θ)/8h 5 2 . The system (3) has the convenient feature that the center of mass (cm.) and relative motion can be decoupled. In terms of the Hamiltonian (3) separates into H = H x + H X , where the cm. Hamiltonian H X = −1/2d 2 X + X 2 /2 is exactly solvable and the relative motion Hamiltonian is given by The Taylor expansion of U ( √ 2|x|) around = ∞ gives √ 2|x| −3 , so that the relative motion is governed in the strictly 1D limit by In this paper, we discuss the effect of both the anisotropy parameter and the interaction strength g on the entanglement in the bosonic and fermionic ground states. Moreover, we derive, within the framework of the harmonic approximation (HA), closed form expressions for the natural orbitals and their occupancies in the g → ∞ limit. Furthermore, using these results, we analytically calculate the corresponding asymptotic bosonic and fermionic von Neumann (vN) entropies and compare their values with the ones obtained numerically. This paper is structured as follows. In Sect. 2, we discuss the entanglement characteristics of quasi-1D systems of two spinless identical particles. Section 3 is devoted to the results, and some concluding remarks are made in Sect. 4.

The Measure of Entanglement
The tool that is usually used to characterize quantum entanglement is the spectrum of the one-particle reduced density matrix (RDM) [36]. For systems of a quasi-1D geometry composed of two spinless particles, the RDM takes the form The eigenvectors and eigenvalues of the RDM are sometimes called the natural orbitals and the occupancies, respectively, and we shall do the same. The normalization conditions for the occupancies of the bosonic (+) and fermionic (−) states give l λ (+) l = 1, and 2 l λ where the factor 2 comes from the fact that the occupancies of the fermionic state are doubly degenerate. As discussed in [35], the state of two identical bosons is non-entangled only in two cases, i.e., when its Schmidt number (Sn) is 1 or 2, which corresponds to λ = 0.5, respectively. As to a fermion state, it is non-entangled if, and only if, its total wavefunction can be expressed as one single determinant [35], which in turn corresponds to λ (−) i = 0.5. We will measure the amount of the entanglement via the vN entropy, where S = −Tr[ρLog 2 ρ] is the ordinary vN entropy [36] and S 0 = 0 and S 0 = −1 stand for the bosonic (B) and fermionic (F) states, respectively. In terms of the occupancies, we have l . The measure (8) vanishes for the non-entangled points discussed above, except for the bosonic states with Sn = 2, (S B = 1). As noted in the already cited [35], the vN entropy alone is insufficient to distinguish whether the bosonic state is entangled or not, since the following may happen: S B = 1 for the state with Sn different than 2. For example, such a situation was observed in the system of bosons interacting via the short-range contact potential confined in a split trap [31].

The Weak-Interaction Limit
In the limit as → ∞ (the strictly 1D case) the system of bosons gets fermionized for any g = 0 [2,7], which is due to the singular behavior of the interaction potential |x 2 − x 1 | −3 at x 1 = x 2 . In this limit, the two bosonic ground-state wavefunctions approach, as g → 0, the modulus of the Slater determinant ψ where ϕ n are the single-particle orbitals of the ideal system (g = 0). The above function is nothing else but a TG wavefunction, which means that the 1D dipolar bosons form a TG gas in the weak interaction regime. We have determined the value of the vN entropy associated with ψ g→0 B by calculating the occupancies numerically through a discretization technique (see, for example, Kościk and Okopińska [29]). The value obtained by us S g→0 B ≈ 0.9851 agrees well with that reported in [32] (0.984). On the other hand, in the non-interacting case g = 0, the bosonic and fermionic ground-state wavefunctions are given by a simple product ψ , respectively, and these states must be regarded as non-entangled (their corresponding vN entropies vanish). One can importantly conclude that the vN entropy of the fermionic ground-state is always continuous at g = 0, which is in contrast to the vN entropy of the bosonic ground-state, which exhibits a discontinuity at this point as → ∞.

The Strong-Interaction Limit
Now we come to the point where we explore the limit of g → ∞ by use of a scheme developed in [29]. Due to the long-range nature of the DDI, the larger is g, the larger is the average distance between the particles. Hence, and because of the fact that the Taylor expansion of U ( √ 2|x|) around |x| = ∞ gives √ 2|x| −3 , one can conclude that at very large g the distance between the particles is large enough so that the interaction among them does not depend on . Thus, as long as we are interested in the regime when g → ∞, we can focus only on Eq. (6).
To begin with, we expand the relative motion potential in (6), V (x) = x 2 /2 + g √ 2|x| −3 , around its local minimum x c = 2 1/10 (3g) 1/5 , retaining only the terms up to second order, V (x) ≈ V (x c ) + 5(x − x c ) 2 /2. The relative motion Schrödinger equation with such a potential is of course exactly solvable and approximations to the lowest even (+) and odd (−) wavefunctions can be constructed as follows where C (±) (g) is the normalization factor that tends to 5 and then changing the variables back to x 1 and x 2 , we obtain final forms of the approximations to the lowest total symmetric (+) and antisymmetric (−) wavefunctions as with Subsequently we translate the coordinates by Due to the fact that the above function is symmetric under permutations of coordinates, it always has a Schmidt decomposition in the form [36]ζ where v k |v l = δ kl . By changing the variables back in (14), namely, byx 1 Finally, after substituting the expansions of ζ(x 1 , x 2 ) and of ζ(x 2 , x 1 ) into (11), we get It is easy to infer that in the limit as g → ∞, where x c → ∞, the integral overlap between the orbitals L k (x) and R l (x) vanishes for any k, l, L k |R l = 0. Hence, bearing in mind that L k |L l = R k |R l = δ kl , one can infer that the family {L l (x), R l (x)} forms asymptotically an orthonormal set as g → ∞. In this limit the orbitals L l (x) and R l (x) are thus nothing else but the asymptotic natural orbitals with the occupancy λ l related to the corresponding coefficient k l by λ l = k 2 l (a two-fold degeneracy occurs). The bosonic and fermionic ground states have thus as g → ∞ the same set of occupancies. Up to this point, our derivations are quite similar to those in [29], where the system of two Coulombically interacting particles was treated by the HA. Here we go one step further and derive closed form analytical expressions for the asymptotic occupancies and their natural orbitals basing on a methodology of [20], wherein a novel derivation of the occupancies of the analytically solvable two-particle Moshinsky model was given.
Following [20], we start with Mehler's formula: where H(l; .) is the lth order Hermite polynomial. Now it should be clear that the orbitals v l and their coefficients k l appearing in Eq. (14) can be found in closed forms by matching Eq. (13) with (16). Indeed, in the limit as g → ∞ that we are interested in, only (C (±) (g) → 5 ] 2 , we compute the Rényi entropy [37], that is, where the factor 2 comes from the normalization condition (2 l=0 λ g→∞ l = 1). It is easy to check that in the limit as q → 1, Eq.
and then, by taking the limit as q → 1, we arrive at S g→∞ B 1.24939. As discussed in Sect. 2, the asymptotic fermionic and bosonic ground states share the same set of occupancies, so that the values of their vN entropies differ from each other only by one, S g→∞ F 0.24939. In the next subsection we will confirm the correctness of our analytical results by comparing them with the results obtained numerically.

Numerical Analysis
As far as we know, except for the cases discussed in the previous subsections, neither analytical solutions to Eq. (5) nor to Eq. (6) are known and we have to resort to numerical methods. Here we apply the Rayleigh-Ritz (RR) method, that uses as a variational trial function a linear combination of a finite set of some basis functions, Since the interaction potential in Eq. (5), U ( √ 2|x|), is finite at the origin x = 0, the basis of the harmonic oscillator (HO) eigenfunctions appears to be appropriate in this case. We found it to work well over the whole range of values of and g. On the other hand, when it comes to Eq. (6), the interaction potential in it, ∼ |x| −3 , Table 1 The lowest approximate energy E (N ) 0 (γ opt ) determined using the optimized RR method as discussed in the text diverges at x = 0, and the interaction integrals are not convergent in the HO basis. Fortunately, this divergence can be cured by using the pseudoharmonic oscillator basis [38] 0 < x < ∞, (u n (0) = 0), where 1 F 1 is the Kummer confluent hypergeometric function and γ is a nonlinear parameter (γ > 3/2), which, as we have verified, has a strong impact on the rate of the convergence of the RR method.
We determine the value of the parameter γ according to an optimization strategy based on the stationarity of the trace of the RR matrix, H RR [39], Having the RR wavefunctions φ i (x) obtained in that way (where φ i = 0 for x < 0), we can construct the even (+) and odd (−) approximate solutions of (6) by putting ψ (±) . For the demonstration of the convergence, Table 1 shows the numerical results for the ground-state energy of (6) together with the corresponding optimal values of γ opt , for different values of g and N . As can be seen, the convergence is fairly good over the full interacting regime, getting increasingly better with an increase in g. It is apparent from the results that the occurrence of the TG regime, which is manifested by the closeness of the relative motion energy value to 1.5, takes place at a very small value of g, that is, at about g = 0.0001. and for S g→∞ F obtained in the previous subsections are marked by horizontal lines. As may be seen, the vN entropy increases, attains a local maximum, and then decreases, and in the limit of large g saturates at a value that is insensitive to . As a matter of fact, the last point has been already explained at the beginning of the Sect. 3.2. As the figure indicates, our analytical results concerning the limit g → ∞ are in excellent agreement with the numerical ones, which confirms their correctness and thereby proves the applicability of the HA to the case of strongly correlated dipolar particles. The dependence of the entanglement on is stronger in the bosonic case than in the fermionic one, which can be attributed to the fact that the probability of finding the bosons at the same place is, in contrast to fermions, affected by in particular. This effect is less pronounced at large g, where the bosons become spatially separated simulating thus the Pauli exclusion principle. Both in the bosonic and fermionic case, the vN entropy is largest in the 1D limit and its value attained in this limit by a local maximum is maximal. We see from Fig. 1 the TG regime starting to occur at a very small value of g, which coincides with our earlier conclusion drawn from the numerical results for the relative motion energy (Table 1). When it comes to the points appearing in the behaviour of the bosonic vN entropy at which S B = 1, we have verified their Sn are essentially different from 2, so they must be considered as truly entangled.
We close our discussion with a comparison of the present results obtained for the systems with the DDI with the ones obtained in [33] for quasi 1D systems of two harmonically trapped particles with a Coulomb interaction. Here we refer only to the strong correlation regimes of both systems, comparing them regarding their entanglement features. In particular, in [33], the value of the vN entropy of strongly interacting charged bosons was found to be about 1.14. Comparing this value with that obtained in this paper, S us to the conclusion that strongly interacting dipolar bosons are more entangled than the strongly interacting charged ones.

Summary
We carried out, within the SMA, a comprehensive study of the entanglement properties of systems of two dipolar particles confined in a harmonic trap. Our results show the effect of the anisotropy parameter on the entanglement in the bosonic and fermionic ground states over the whole range of values of g. Moreover, within the framework of the HA, closed-form analytical expressions for the asymptotic natural orbitals and their occupancies have been obtained by use of Mehler's formula. The corresponding vN entropies of the asymptotic bosonic and fermionic ground states have been also derived analytically and it has been shown that their values are in excellent agreement with the results obtained numerically.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.