The momentum distribution of two bosons in one dimension with infinite contact repulsion in harmonic trap gets analytical

For a harmonically trapped system consisting of two bosons in one spatial dimension with infinite contact repulsion (hard core bosons), we derive an expression for the one-body density matrix $\rho_B$ in terms of centre of mass and relative coordinates of the particles. The deviation from $\rho_F$, the density matrix for the two fermions case, can be clearly identified. Moreover, the obtained $\rho_B$ allows us to derive a closed form expression of the corresponding momentum distribution $n_{B}(p)$. We show how the result deviates from the noninteracting fermionic case, the deviation being associated to the short range character of the interaction. Mathematically, our analytical momentum distribution is expressed in terms of one and two variables confluent hypergeometric functions. Our formula satisfies the correct normalization and possesses the expected behavior at zero momentum. It also exhibits the high momentum $1/p^4 $ tail with the appropriate Tan's coefficient. Numerical results support our findings.


I. INTRODUCTION
Two-body models have the merit that they can be solved exactly and give direct access to the wave function and to the one-body density matrix. We consider here two bosonic atoms with mass m in one dimension with infinite repulsive contact interaction (two hard core bosons or the so-called Tonks-Girardeau regime). The system is also subjected to an external harmonic potential trap with frequency ω. A system of Tonks-Girardeau (TG) gas has been already realized in ultra-cold experiments [1]. On the theoretical side the TG model is exactly solvable through the Bose-Fermi mapping theorem, which relates this gas to a system of noninteracting spin polarized fermions [2, 3]. The two systems -bosonic and fermionic -exhibit some identical properties such as the particle density. However, their corresponding one-body density matrices are not the same resulting into a considerable difference between the momentum density of TG gas from that of an ideal Fermi gas [4]. While it is rather easy to calculate the momentum density of an ideal system of fermions, the task is difficult for a TG gas. The main purpose of the present work is to provide a definitive solution for the case of two particles, deriving a closed form expression of the momentum distribution for two harmonically trapped interacting bosons in the TG regime. The analytical approach allows us to clearly identify the deviation from the fermionic counterpart.
(2) built on the two lowest single particle normalized eigenfunctions ϕ 0 (x) and ϕ 1 (x) of the harmonic oscillator potential trap. Explicitly we have where = /mω denotes the harmonic oscillator characteristic length. After insertion in (1), we get A few words about this wave function are in order. It may also be considered in the context of one dimensional spinor quantum gases studies, a topic which has recently received a lot of attention [5][6][7]. Let us consider the particular case of two ultra-cold bosons in the Tonks-Girardeau regime occupying or prepared in two hyperfine states; this constitutes a pseudo-spin-1/2 system of different spins (which we denote by ↑ and ↓), the wave function ψ B (x 1 , x 2 ) = |ψ F (x 1 , x 2 )| in (4) being the spatial wave function of a spin-1/2 hardcore bosons corresponding to the full wave function of the system, Φ = |ψ F (x 1 , x 2 )| ⊗ χ S=1 , where χ S=1 = (|↑↓ + |↓↑ ) / √ 2 is the triplet (symmetric) spin state [8]. Another possible configuration [8] would be to use a singlet (antisymmetric) spin state χ S=0 = (|↑↓ − |↓↑ ) / √ 2. In this case the full wave function would simply read Ψ = ψ F (x 1 , x 2 ) ⊗ χ S=0 , and calculations are very easy. In the present investigation we consider the much more difficult triplet spin state case for which the calculations are by far not straightforward because of the presence of the absolute value in the wave function.
The plan of the work is the following. In sect. II we use the two-particle wave function in (4) to derive an expression of the one-body density matrix in terms of centre of mass and relative coordinates. This first result is then used in sect. III to derive a closed form for the momentum density which holds for arbitrary values of the momentum p. This second result is analyzed and tested in sect. IV along three lines: (i) we check that one recovers the correct value of the momentum density at p = 0; (ii) we check that our density satisfies the required normalization condition n B (p)dp = 2; (iii) we prove analytically, and we verify numerically, that our momentum distribution exhibits a well-known result in quantum gases with contact repulsive interactions [9], namely that the momentum distribution of the atoms decays asymptotically as 1/p 4 for large momentum p, and as a consequence we derive the exact expression of the so-called Tan's contact coefficient for the case of two interacting particles [10]. In the last section the results of the paper are summarized.

II. ONE-BODY DENSITY MATRIX OF A ONE-DIMENSIONAL SYSTEM OF TWO BOSONS WITH INFINITE CONTACT REPULSION IN HARMONIC TRAP
The one body density matrix for the system under study is given by [11] Using the explicit expression (4) and denoting z 1 = x 1 / , z 2 = x 2 / , ξ = x / , the density can be written as where u is the dimensionless coordinate of the centre of mass of the two particles and v is half the dimensionless relative coordinate of the two particles: Performing now the change of variable w = (ξ − u), eq. (6) takes the form To carry out the above integration, we remove the absolute value symbol and write the density matrix as where The three functions A(u, v), B(u, v) and C(u, v) are calculated in Appendix A in terms of the error function erf(x) = 2/ √ π x 0 e −t 2 dt. Substituting the result (A7) into (9), we obtain the one-body density matrix ρ B (x 1 , x 2 ) in the form For two noninteracting fermions, the one-body density matrix is given by where we have used the explicit expressions (3) of ϕ 0 (x), ϕ 1 (x). Remark that these densities are correctly normalized to the particle number, so that ρ B (x, x)dx = ρ F (x, x)dx = 2. By comparing the bosonic to the fermionic density, we may write and thus clearly identify the deviation Note that if x 1 = x 2 = x then u = x, v = 0, and we obtain ρ D (x, x) = 0 and therefore ρ B (x, x) = ρ F (x, x), as required by the Bose-Fermi mapping theorem [3]. Expression (13), together with the explicit deviation (14), constitute the first outstanding result of the present study. It should be pointed out that, although the Tonks-Girardeau gas and the ideal Fermi gas have identical local spatial density ρ B (x) = ρ F (x), they display different one-body density matrices yielding to very different momentum density profiles [12]. In the sequel, this property will be explicitly shown and analytically verified in the case of two particles.
Before proceeding further, it is interesting to note that, at this level, one can perform an expansion of the density matrix in (13) with (14) in terms of |v| = |x 1 − x 2 |/(2 ). This corresponds to a short-distance expansion and allows one to directly access to momentum density profile at momentum high-values [13]. This asymptotic behaviour is of particular interest in the study of short-range interacting systems [14].

III. CLOSED-FORM EXPRESSION OF THE MOMENTUM DISTRIBUTION
In what follows we will derive an analytical expression of the momentum distribution for the system of two bosons in the Tonks-Girardeau regime. In terms of the one-body density matrix, the momentum density denoted by n B (p) is given by [11] and is normalized to the total particle number, i.e., n B (p) dp = 2. In terms of centre of mass and relative dimensionless coordinates (u, v), we have Now, substituting (13) into (16), we can write which is nothing but the momentum density of two noninteracting fermions in harmonic trap, and which is the momentum density corresponding to the deviation ρ D . Next we are going to evaluate the integrals in (18)- (19). It will be convenient to use the momentum dependent dimensionless variable q = p/ .

A. Evaluation of nF (p)
The integration with respect to u in (18) yields By using the eq. (3.952), number (9) of ref. [15] ∞ 0 the integrals in (20) are easily computed, giving the final expression of the momentum density n F (p): It is easy to check that it is normalized as n F (p) dp = 2. The same n F (p) result can be derived also by using the expression, n F (p) = | ϕ 0 (p)| 2 + | ϕ 1 (p)| 2 , where ϕ 0,1 (p) denote the single particle wave functions in momentum representation, i.e., the Fourier transforms of the single particle wave functions ϕ 0,1 (x) [16].

B. Evaluation of nD(p)
Let us start with the second term of the first double integral I 1 in (19), and perform an integration by parts with respect to u, to get Combining now the two terms, and integrating over u (simple Gaussian-type integral [15]), the first double integral in (19) simplifies into For the second double integral I 2 in (19), we use the property erf(−x) = −erf(x), to write Although it is not at all straightforward because of the presence of the error function, we were able to proceed again by reducing the double integral (25) to a single one (see result (B10) in Appendix B): Collecting the results (24) and (26), the momentum density (19) may be written as in terms of three single integrals we proceed with their evaluation.

C. Closed-form expression of the momentum distribution nB (p)
Collecting n F (p) from (22) and the result in (40), we obtain the closed-form expression of the momentum distribution which is the principal result of our paper. In the next section, we provide strong tests that support its validity. We have computed n B (p) for any value of p using either of the two series representations (36) or (37). For the evaluation of the Ψ 1 function, from a numerical convergence point of view it turns out that the series in terms of confluent hypergeometric functions is more convenient because |z 1 | < 1 while |z 2 | can take any positive value. The numerical results for n B (p) has been further confirmed by using (27) and direct numerical integration of (28), (29) and (30). A plot of n B / as a function of q = p/ is shown in fig. 1, and is compared to both contributions n F and n D . We clearly observe a substantial shape difference between the fermionic and bosonic cases. The fermionic distribution presents a minimum at q = 0, a maximum at q = 1/ √ 2, and decreases rapidly to zero due to the e −q 2 term. The bosonic distribution, on the other hand, has almost double value at q = 0 (see ratio n B (0) /n F (0) ≈ 2.02 discussed in sect. IV A), and decreases thereafter but in much slower fashion; approximately beyond q > 3, we have n F 0 and n B is essentially due to the deviation n D .

IV. TESTS OF THE MOMENTUM DISTRIBUTION EXPRESSION (41)
We consider now three tests on the analytical expression of the momentum distribution n B (p). First, we will verify that the correct value of n B (p = 0) is obtained, then we will check the normalization condition, and finally we will make sure that we recover the correct high-p asymptotic behavior.
A. Test of nB(p) at zero momentum p = 0 We wish to compare the value of momentum density n B (p) at p = 0 by using its definition in (15) with that obtained through our general expression (41). While the idea of the test is simple, the calculations are not straightforward.
Using the definition (15) together with ρ B (x 1 , x 2 ) given by the first eqality in (6), gives for p = 0 For any fixed value of ξ, we make the change of variables z 1 = r + ξ, z 2 = s + ξ, and get After the integration on ξ is performed, the above expression reduces to Using the identity (see eq. (3.562) number (4) of [15]) The remaining integral can be found in the literature [20], and after simplifications we get It is interesting to compare this value of the momentum density at p = 0 (peak momentum distribution) with its two-fermions system counterpart given by (22) We clearly observe a large deviation measured by the ratio n B (0) /n F (0) = 4 π arctan √ 2 + 4 √ 2 π − 1 ≈ 2.02. This ratio is featured in fig. 1.

B. Normalization condition of momentum distribution nB(p)
Next, we want to show that momentum density obeys the normalization condition n B (p) dp = 2. It is more convenient not to use the form of n B (p) in (41), but an expression obtained a few steps before. That is, we use the momentum density as given in (17), n B (p) = n F (p) + n D (p), with the expressions for n F (p) and n D (p) given by (20) and (27), to write Taking into account the following representation of the Dirac delta distribution the integral over the momentum p (or equivalently over q = p/ ) is quite straightforward where the final result is obtained using erf(0) = 0. We have thus checked that, for the two-particle system considered in this work, we have the correct normalization condition for n B (p).

C. High-p behavior of the momentum distribution and Tan's coefficient
Here, we will examine the asymptotic behavior of the momentum density n B (p) obtained in (41) as p → ∞ and we will prove that where C 2 is the exact Tan's contact coefficient for the two-atom system under study [10], recovering the 1/p 4 dependence originated from the short-range character of the interaction.
To derive this result we search the asymptotic behavior up to 1/p 4 . For this purpose we need the asymptotic behavior (x → +∞) [22] To determine the asymptotic behavior of Ψ 1 1 + j, 1/2, 3/2, 1/2; −1/2, − 2 p 2 / 2 as p → ∞, with either j = 0 or j = 1, we use the second expression in (36) to write and then the asymptotic expansion (56) to get To retain only terms up to 1/p 4 , we have to consider the summation indexes up to m + k ≤ 1 − j. It is easy to show that Collecting expressions (57), (60) and (61) into (41), after the cancellation of the terms of order 1/p 2 , one finds which is the result announced in (55): we have recovered the well known 1/p 4 asymptotic decay of n B (p) for high momentum values. The coefficient in front of the high-p tail of the momentum distribution, called Tan's coefficient [10], is given by (2/π) 3/2 (mω ) 3/2 for the system of two particles under study. Note, finally, that from a basic dimensional analysis, expressions in (62) or (55) have correct units of momentum density. The high-p tail behavior can be observed in fig. 2 where we show n B (p), n F (p) and n D (p) multiplied by p 4 /C 2 . One can appreciate how n B (p) tends to 1 for large q, and that the asymptotic limit is dominated by the deviation term n D , since the n F is already negligible for q > 4.

V. CONCLUSION
In this paper we have studied a one-dimensional quantum system consisting of two bosons harmonically trapped and with an infinite contact repulsion. This system has been the subject of numerous studies, in particular with respect to its asymptotic momentum distribution. In the present work we have been able to provide an exact expression of the momentum distribution valid for arbitrary value of the momentum, and thus we have extended the list of exactly solvable models of two interacting particles. Starting from the wave function of the two particles, we have derived an exact expression of the one-body density matrix expressed in terms of centre of mass and relative coordinates. From this result we have calculated the corresponding momentum distribution. To clearly identify the deviation of this system from the two noninteracting fermions, we wrote this momentum density as the sum of two terms: a first one corresponding to the two noninteracting fermions, and a second one related to the short range character of the interaction.
In order to validate our analytic expression of the momentum density, three robust tests have been performed satisfactorily. In particular, the third test concerns the asymptotic decay for high momentum values; we have been able to derive the exact expression of the Tan contact coefficient, which we recall is the coefficient in front of the high-p tail of the momentum distribution, for the case of two particles [10].

(A7)
Appendix B In this Appendix we consider the double integral I 2 in (25) which, for convenience, we have split into three double integrals defined through In order to compute the u integral in (B4), we first perform an integration by parts to get The second integral on the right-hand side of (B5) can be carried out, giving Hence, substituting (B6) into (B4), one finds out that K 3 is related to K 1 through To carry out the u integrals in (B2)-(B3), we start from the identity Using the above relation for a = |v| allows us to write the expressions of (B2) and (B3) as single integrals, so that Finally adding the three results (B9) and (B7), we can write (B1) as ∞ 0 e −3v 2 /2 v cos(2qv) dv .