Off-diagonal correlations of lattice impenetrable bosons in one dimension

We consider off-diagonal correlation functions of impenetrable bosons on a lattice. By using the Jordan–Wigner transformation the one-body density matrix is represented as a (Toeplitz) determinant of a matrix of fermionic Green functions. Using the replica method we calculate exactly the full long-range asymptotic behaviour of the one-body density matrix. We discuss how subleading oscillating terms, originating from short-range correlations, give rise to interesting features in the momentum distribution.


3
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT effects. This result allows one to determine the main features of the momentum distribution of impenetrable bosons on a lattice.
The paper is organized as follows: in section 2 we use the Jordan-Wigner transformation to obtain the determinantal representation for the one-body density matrix of impenetrable bosons. We discuss how the results ofVaidya and Tracy are obtained by considering the leading asymptotic behaviour of the determinant. In sections 3 and 4 we present the replica method and derive our main result, the full asymptotic expansion of the one-body density matrix. Section 5 contains conclusions and prospects. The mathematical details are given in appendices.

One-body density matrix as a Toeplitz determinant
Let a m , a † m be creation and annihilation operators of bosons on an infinite lattice with sites labelled by the index m. In the limit of strong repulsive on-site interactions the Hilbert space is projected onto a subspace of states with at most one particle in each site. Using the Jordan-Wigner transformation The one-body density matrix is defined by the following ground-state expectation value: where we have used the relation and the fermionic commutation relations (2). For R = 0 the one-body density matrix is related to the mean number of particles per site (filling factor) ν as while for R > 1 we use the relation (4) for every factor in the product in equation (3) and represent the one-body density matrix as an expectation value of 2R operators: where we have defined A l = c l + c † l , B l = c l − c † l . We now use the Wick theorem and take into account that free fermionic expectation values are given by where the Green function of free fermions is Here n q = n −q is the ground state momentum distribution of free fermions in a lattice with filling factor ν, so that Then the one-body density matrix (6) is represented as a determinant of the R × R matrix The determinant (9) has identical elements along each of its diagonals and thus belongs to the class of Toeplitz determinants: where g(q) is called the generating function (symbol). It is easily shown from (7) that in our case the generating function is It has two jump discontinuities at Fermi points q = ±πν. The large-size asymptotic behaviour of Toeplitz determinants with such singular generating functions has been the subject of extensive studies in mathematics (for overview and references see [19]). It has been shown for a certain class of singular generating functions that the large size asymptotics of Toeplitz determinants is correctly described by the Fisher-Hartwig conjecture [19]. In the case of generating function with two jumps as in equation (11) the analysis of the asymptotic behaviour of the Toeplitz determinant (9) is subtle and the Fisher-Hartwig asymptotic formula remains a conjecture. Being applied, it yields the following long-distance behaviour of the one-body density matrix: where ρ ∞ = πe 1/2 2 −1/3 A −6 , and A = exp(1/12 − ζ (−1)) 1.2824271 is Glaisher's constant related to the Rieman zeta function. Numerical computation of the determinant (9) for R ∼ 100 and various values of the filling factor shows that equation (12) provides an excellent estimate for the dominant smooth behaviour of the one-body density matrix at large distances. Subtracting this dominant contribution from the exact numerical expression reveals interesting oscillating corrections. These corrections are sensitive to the specific value of the filling factor and reflect an interplay between short-distance interparticle correlations and lattice effects. The Fisher-Hartwig conjecture is unable to capture this behaviour. In the next section we present the calculations of the asymptotic long-distance behaviour of the one-body density matrix, based on the replica method which has been recently applied to study correlation properties of integrable models. We shall see that this method is able to provide us with the dominant behaviour (12) as well as with the oscillating terms.

Replica calculations of the one-body density matrix
Below we describe the replica method for finding the large-R asymptotics of the determinant (9). By standard manipulations the determinant in (9) is transformed to the following R-dimensional integral: where is the Vandermonde determinant. The notation of average in (13) is justified by the fact that the measure of integration coincides with the distribution of eigenvalues of random unitary matrices drawn from Dyson's Circular Unitary Ensemble [20]. Denoting exp(iq l ) = z l and exp(iπν) = v and using the representation (11) we arrive at an equation The replica evaluation of the average (13) begins with modifying (replicating) the function being averaged in (15) by taking the 2n-th power of absolute values. Then the original expression is recovered by taking n → 1/2. It is however assumed throughout the calculations that n is an integer (and 2n is even) so that many simplifications occur. Eventually, the limit n → 1/2 will be taken with care and we shall describe it in detail.
For an integer n, a straightforward algebra shows that we are dealing with an average of a rational function:

Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
The crucial fact is the existence of a duality transformation relating the R-dimensional integral (16) to a 2n − 2 = m-dimensional integral: where The proof of the duality transformation is presented in appendix A. Note that in the dual representation (17) large distance R appears only as a parameter, which makes the dual representation an excellent starting point for the asymptotic expansion of g 1 (R). The main contribution to the integral (17) comes from the end-points of integration. Let p be a number of variables in the vicinity of x − = 0 and p = m − p be a number of variables close to x + = 1. We shift the integration variables and approximate the integrated high-degree polynomial as , The integration measure in (17) factorizes as Then, summing over all possibilities to distribute m variables among saddle points and performing the remaining integration over ξ c , ξ d we obtain where the integrals are performed using the Selberg integration formula Combining them with the constant (18) we observe that each term in the sum (21) is proportional to the total combinatorial factor: is a familiar combination from previous studies of replica theories [14,15]. Now we are at a position to take the limit n → 1/2. It is important to notice that for even m = 2n − 2 there is always a central term p = p = m/2 in the sum (21) which provides a smooth non-oscillatory contribution. It is characterized by a maximal degree of replica symmetry breaking, since it originates from the saddle point for which the number of components x d ∼ x − is equal to the number of components x b ∼ x + . This term will be retained in the limit n → 1/2 (m → −1) and this provides an operational definition of the correct analytic continuation.
Making a change of summation index from p to m/2 + k = n − 1 + k we factorize the combinatorial coefficient (24) (26) This leads to an expression and the summation over k in equation (27) was extended to all integers. Now we take the limit n → 1/2 term by term in the series (27). The only complication appears with the overall numerical factors C n and A n . The analytical continuation of these factors is explained in appendix B where we show that A 2 1/2 = √ 2ρ ∞ and C 1/2 = 1/R. We also have (1/2 + k) (1/2 − k) = π(−1) k . We therefore find that The first several terms in this asymptotic expansion of the one-body density matrix read where the leading term coincides with the result (12). As is seen from equation (30), for small k the coefficients in the series (29) rapidly decrease with increasing k. However, they will diverge for large k, which shows that we are dealing here with asymptotic rather than convergent series.
The result then is which leads to the following asymptotic expression for the one-body density matrix: Our result reproduces correctly the continuous limit. In this case it is useful to introduce a length scale, the lattice constant a, so that x = Ra is the distance and πν/a = k F is the Fermi wavenumber proportional to the density. Going to the continuous limit corresponds to having the filling factor ν tending to zero while the product k F x = πνR is kept fixed. Dividing g 1 by the lattice constant to have a correct normalization, equation (30) becomes in agreement with the Vaidya and Tracy result [21]. The sign of the oscillating cosine term in equation (36) has been corrected as discussed in [17].

9
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT

Conclusions and prospects
The result (35) provides an exact analytical expression for the long-distance behaviour of the one-body density matrix for impenetrable bosons on an infinite lattice. It is valid for all values of the filling factor ν, in particular it reproduces correctly the continuous limit ν → 0. The structure of equation (35) is in accordance with the hydrodynamic expansion conjectured by Haldane [22]. Equations (25), and (26) (at n = 1/2) provide the exact values for the non-universal coefficients of the leading oscillating terms. The oscillatory terms in (35) reflect the physics on short distance scales, of the order of the mean interparticle separation. They are analogous to the Friedel oscillations in the physics of fermions and provide yet another manifestation of fermionization. The interplay between these short-distance correlations and lattice effects is interesting, for instance, in the case of half filling (ν = 1/2), where the particles tend to occupy every second site. This results in the oscillations of the one-body density matrix with the period of the lattice.
The momentum distribution cannot be reconstructed from the asymptotic expression (35) for the one-body density matrix, as it requires the knowledge of the short-distance behaviour. However, we can predict some of the important features. For sufficiently small momenta their distribution diverges as an inverse square root of the momentum. The finite size would result in the finite peak value of the momentum distribution proportional to the square root of the total number of particles. For momenta of the order of 2k F = 2πν/a the derivative of the momentum distribution has a jump (the distribution itself has a cusp) resulting from the underlying Fermi surface of the fermionized bosons. Weaker singularities exist for higher multiples of 2k F . It is again interesting to see that for the case of half filling the jump in the derivative occurs just at the border of the Brillouin zone and therefore cannot be observed: the momentum distribution has no singularities.
These features can be probed in current experiments with cold atoms. The most promising method is expected to be a combination of Bragg spectroscopy with the time-of-flight technique. In this case a Bragg pulse is applied after switching off the trap, when the density of the cloud is already sufficiently small so that the number of scattered atoms (proportional to the momentum distribution) can be detected with precision [23]. This technique can allow the observation of the singularities such as cusps in the momentum distribution originating from the oscillating subleading terms in the coordinate correlation functions.
There are several directions in which the present work can be extended. An immediate question concerns the finite value of the interactions or, in other words the case of the soft-core bosons. This case is similar to the case of the Heisenberg spin chain, for which there are recent results based on the Bethe ansatz solution [24]. It is not clear at the moment how to obtain the explicit behaviour of the correlation functions away from the free fermionic point (equivalent to the hard-core bosons) where the Toeplitz determinant expression exists for the one-body density matrix. In the case of impenetrable bosons, our work can be extended to the domain of timedependent correlation functions which contains information about elementary excitations and are related to the response of the system to external perturbations.
with the following 2-fermionic Green function: We now take the limitsū 2n → 0, v 1 → ∞ using G(ū 2n v) = G(0) = 1 and G(ūv 1 ) (ūv 1 ) R+2n−1 . The factor v R+2n−1 1 common to the elements in the first column cancels with the corresponding term in the prefactor of the second line of equation (A.1). Taking out a factor u R+2n−1 c from each row we arrive at the following representation: . . , u m ). We take the singular limit u c → u and v c → v by a standard procedure, described, for example in [25]. Expanding in non-zero elements of the last row and first column leads to expressing equation (A.6) as a determinant of (2n − 2) × (2n − 2) matrix of partial derivatives: Now we use the following integral representation in which we recognize the expression for the hypergeometric function: , (A.9)

12
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT where F k (s) = [(s − 1)∂ s ] k F(−R − 2n + 3, 2; 4; s). In this function k derivatives ∂ s can be commuted through the factors s − 1 to the right, producing additional terms which has at most k − 1 derivatives. These terms do not contribute to the determinant, since they already appear in the k − 1-th column. Next, we use the known properties of hypergeometric functions to show that (1 − s) R+2n−1 F(R + 2n + 1, 2; 4 + k; z).
Substituting this expression into (A.9) and taking the limit s c → s = 1 −ūv we get In this expression we used the properties of hypergeometric functions [26] along with its integral representation. Substituting this result into (A.10) yields