Continuum limit Tonks-Girardeau matrix elements. Part I. The ground state and the uniform density state

The Tonks-Girardeau model is a quantum mechanical model of N impenetrable bosons in 1+1 dimensions. A Vandermonde determinant provides the exact N -particle wave function of the ground state, or equivalently the matrix elements with respect to position eigenstates. We consider the large N limit of these matrix elements. We present a binning prescription which calculates the leading terms of the matrix elements in a time which is independent of N, and so is suitable for this limit. In this sense, it allows one to solve for the ground state of a strongly coupled continuum quantum field theory in the field eigenstate basis. As examples, we calculate the matrix elements with respect to states with uniform density and also states consisting of two regions with distinct densities.


Introduction
A quantum state is completely characterized by its matrix elements with any basis of the Hilbert space. These matrix elements generalize the notion of wave function from quantum mechanics to more general quantum systems. The Coordinate Bethe Ansatz (CBA) [1] computes these matrix elements for the coordinate basis, in which the degrees of freedom have concrete values at distinct spatial points. The CBA only exists for systems with finite numbers of degrees of freedom. However, these systems all have a number N , for example the number of lattice sites in a spin chain or the number of particles in the Tonks-Girardeau [2,3] or Lieb-Liniger [4] models, such that in the formal large N limit one obtains a quantum field theory.
Were it possible to take the large N limit of the CBA, one would, for the first time, have the explicit wave functions for the states of a strongly interacting quantum field theory. This would open exciting possibilities. For example, one could use the equivalence between the high spin XXX spin chain and the CP 1 σ model to understand just how the fractional instanton plasma of [5,6] is realized in Minkowski space. By studying the difference between the ground state and first excited state wave functions, one may achieve an understanding of the nonperturbative generation of a mass gap [7] in a full quantum field theory with the same concreteness that is achieved in our understanding of the mass gap in the double well quantum mechanics model.
So far, such a large N limit of general wave functions has been lacking because the Bethe Ansatz is complicated. It consists of N ! terms, each corresponding to an element of the symmetric group S N . Many approaches have been developed which simplify the calculation of the matrix elements so that it can be done in a time which is only polynomial in N . But this is still too long for the large N limit. So far, to our knowledge, the large N limit of only one matrix element has ever been computed [8].
In this note we will consider the simplest model which allows a Bethe Ansatz solution, the Tonks-Girardeau model. In this model we will introduce a novel binning approximation scheme for evaluating the matrix elements in the large N limit. Our strategy is essentially to decompose S N into cosets for each of which the sum can be computed explicitly. For questions of physical interest, like the roles of fractional instantons or unitons in a mass gap, one is interested in matrix elements of Hamiltonian eigenstates with fairly simple coordinate states, or at least with states whose complexities do not themselves depend on N . We will show that for some families of such simple states, our binning prescription allows us to systematically approximate the matrix elements in a time which is independent of N . Therefore we feel that our method can provide a suitable starting point for a large N limit of these states.
While we suspect that a generalization of our method to any Hamiltonian eigenstate is possible, we restrict our attention to the ground state. This choice is largely driven by the fact that the ground state matrix elements are all given by Vandermonde determinants. This allows us to numerically test our approximation scheme for values of N of order 10 4 . We find that the logarithms of the matrix elements can be expanded in powers of 1/N and our method allows the leading term in this expansion to be estimated. In some cases we attempt to improve the accuracy of this estimation and we find that we can increase the accuracy as desired at the cost of a longer computation. Critically, the length of computation never depends on N , but only on the complexity of the state and the desired precision.
We begin in Sec. 2 with a review of the Tonks-Girardeau model and its solution using the Bethe Ansatz. Next in Sec. 3 we consider matrix elements with states in which the gas has a constant density. First we give an explicit and exact formula for the matrix elements using Vandermonde determinants. Next we introduce our binning approach and use it to calculate a set of approximations of the matrix elements of ever increasing accuracy. Finally, in Sec. 4, we consider matrix elements with respect to states with two densities. We consider two cases, one in which one density vanishes and another in which the region with each density contains the same number of particles. As this case is more complicated than the single density case, we compute only the leading approximation. We find that it calculates the logarithm of the matrix elements with an error which is roughly independent of N , and in general less than 10%. for any j, l. Periodic boundary conditions are assumed

The Model and Its Ground State
for all j.
The periodicity condition (2.4) restricts the space of wave functions to the subspace of free particle states generated by e 2πikx/L with k ∈ Z. The impenetrability condition (2.3) further restricts the space of wave functions to the subspace of Slater determinants of N such plane waves In other words, if one defines the matrix Note that if two k i are equal then the wave function ψ(x) vanishes.
The ground state wave function can and will be taken to be real and positive. It corresponds to the case We will impose the positivity condition by simply taking the absolute value of ψ in each expression Adding a constant to the k i has the effect of multiplying the wave function by an x-dependent scalar which would change the state (2.7) and in general even change the energy. However, the expression (2.9) for the ground state wave function, unlike (2.7), is invariant under this operation because the overall phase is zeroed by the absolute value. Therefore, for simplicity, we will add (N − 1)/2 to k i and so set This still describes the ground state wave function so long as we use the expression (2.9). Note that the k i can only be interpreted as momenta if the wave function is defined by Eq. (2.7) and so we have not actually modified the momenta, or even the state.
With this choice of k i , the wave function can be brought into Vandermonde form. Defining the Vandermonde matrixM the wave function is a Vandermonde determinant

A Bethe Ansatz Formulation
If the x j or k i are evenly spaced then the wave function can be brought into Vandermonde form. Whether or not this is the case, the determinant formulas above can all be expanded using the Leibnitz formula for the determinant. We will now provide the expansion of (2.12), the following proceeds identically for (2.7) and (2.9) Here g is a permutation of the set of natural numbers [1, N ], realized by the bijection P : The symbol (−1) g is equal to 1 for an even permutation and (−1) for an odd permutation.
Let us define where θ is the Heaviside step function. Then This is the form of the Bethe Ansatz [1], which also describes many more interesting systems. The goal of this paper is to introduce a new method to perform the large N limit of the sum (2.15). The fact that this sum admits a Vandermonde form (2.12) will allow numerical checks of our formulas.

A Constant Density Gas
We will begin with the simplest case corresponding to a constant density gas. In Subsec. 3.1 we will use the Vandermonde form (2.12) of the ground state wave function to evaluate the ground state wave function ψ 0 (x) at this point, or equivalently to find the matrix element between the ground state wave function |Ω and the state |x corresponding to a constant density gas (3.1).
Next we will calculate the same matrix element using our binning formalism. We will find that the binning formalism estimates the logarithm of the answer. Successive approximations will be made which calculate this logarithm more and more precisely. The complexity of this calculation is independent of N .

The Exact Solution
Let us define the N th root of unity χ = e 2πi/N .
Then the Vandermonde form of the ground state wave function (2.12) may be simplified as follows Therefore we need to calculate the term β N in the absolute value Here β N is an (N − 1)st order polynomial in the χ. Let γ N,j be the term with j powers of χ We will proceed by calculating all of the γ N,j . As a warm up, let us start with small values of j. Clearly What about j = 2?
Let us defined the reduced summands Now the set {ĩ k } is a two element partition of i − 3. It is conventional to count the number of nonzero elements in a partition, so then this is a partition of at most two elements. The constants a .

(3.11)
Putting this all together where we used the fact that χ N = 1.
is no more difficult. Again one may define the coefficients a We then find and so the matrix element is It may seem odd that the matrix element grows with N , but using our conventions 3.20) and the norm of the |x states is infinite. To obtain probabilities from these matrix elements, of course one would need to normalize the states. This would not be difficult, but our goal is just to calculate the matrix elements.

The Binning Procedure
The dramatic simplifications in the above calculation occur in part because in Eq. (2.14) m j and K j are linear while Ψ jk is quite simple. The simplicity of Ψ jk occurs for all matrix elements in the model, and that of K j occurs for all matrix elements of the ground state |Ω . However, to evaluate the ground state wave function with general particle positions x one needs to relax the linearity of m j or equivalently of x j .
We are interested in the limit of a large number N of particles. In this limit, it makes sense to consider a local density of particles. We do not know whether the wave function will be supported on configurations for which the local density is smooth. Nonetheless, the matrix elements with respect to smooth configurations are the most interesting for physical applications. Our goal is to devise a method to compute such matrix elements.
Our wave function is symmetric under permutations of the x j [3] and so without loss of generality we may order the x j so that x j+1 > x j . Then if the density is smooth, on sufficiently small scales the x j will be linearly distributed in j to any desired precision for j in a sufficiently small interval. Therefore we expect simplifications to occur on such small scales. Our strategy will be to decompose the wave function into factors which depend only on such small ranges, to evaluate these factors and then to reassemble them to obtain an approximate wave function.
More concretely, our procedure for evaluating the sum (2.15) is as follows: Step 1: Divide the ordered set {x j } into q intervals S i called bins.
Step 2: For each bin S i with n i elements, choose a disjoint subset F i ⊂ [1, N ] which also has n i elements. Recall that each g ∈ S N corresponds to a map P : [1, N ] −→ [1, N ]. Let H ⊂ G consist of all elements g ∈ S N such that the corresponding map satisfies (3.21) Step 3: Calculate the terms in Eq. (2.15) corresponding to H ⊂ S N . In other words, sum over only those permutations which fix P (S i ) = F i .
Step 4: Consider all possible sets of images {F i } and sum over the result of step 3 for each summand.
This procedure, if carried out exactly, will produce the exact wave function. However step 4 is computationally difficult and so we will only perform this sum approximately, yielding an approximate answer.

The Two Bin Approximation: The Master Formula
Let N be even. We will divide [1, N ] into two bins We will refer to F as the image, reflecting Eq. (3.23).
We have thus completed steps 1 and 2. The q = 2 intervals are S A and S B and the image P (S A ) is F , whereas P (S B ) is the complement of F in [1, N ]. Next for step 3. Recall that we are still considering the matrix element of the ground state and the constant density state. We need to evaluate To evaluate this, we need to understand H F .
Elements of H F are represented by maps P such that P (S A ) = F . For a given image F , choose one such map P F corresponding to any fixed group element g F ∈ H F . Any other element of H F can be obtained by acting on g F from the right by a permutation of the image F and of its complement [1, N ] \ F . These two permutation groups are each S N/2 and so H F is a right coset Similarly, any element of H F can be obtained by acting on g F from the left by a permutation S N/2 of S A and a permutation S N/2 of S B and so H F is also a left coset Using either representation, the key point is that the two copies of S N/2 commute with each other and so (3.24) can be factored into one expression for each copy represent g A g F and g B g F respectively.
Simplification 1: Up to a total phase, ψ A and ψ B are Vandermonde determinants. This is true not only for the ground state, but also for all excited states. Each ψ I is always a Slater determinant and it is equal to a phase times the Vandermonde determinant whenever the x j with j in the bin j ∈ S I are evenly spaced.
In more general models whose states can be put in Bethe Ansatz form, matrix elements are not given by Vandermonde determinants. However, simple cases in which the x are distributed linearly, such as the Néel state in the XXX spin chain, do have known matrix elements [9,10,8]. Therefore we hope that in general our binning method will reduce the calculation of piecewise linear x matrix elements to that of linear x matrix elements, which at least in some cases can be found via traditional methods.
To bring ψ A into Vandermonde form, we pull out an overall phase where P F A is the map P A : [1, N/2] → F obtained from the action of g F on [1, N/2]. Note that each χ j appears in the product N/2 − 1 times. And so we may pull out the midpoints of the exponents from each factor to obtain Recall that we are free to choose any g F such that P F A (S A ) = F . Let us fix g F by demanding that the P F A (j) are monotonically increasing with j. In this case, the arguments π N P F A (k) − P F A (j) of the sine functions on the right hand side of (3.30) are strictly positive. The factor ψ B can similarly be put in Vandermonde form This is identical to the expression (3.29) for ψ A except for an overall F -dependent sign. Therefore (3.32) Again, we will choose g F such that P F B is monotonically increasing, so the arguments of the sine functions will be positive.
Multiplying ψ A by ψ B , the F -dependence drops out of the χ term, as the j sum is now over the entire range [1, N ]. Therefore the F -dependence of the phase of ψ F lies only in the overall signs (−1) Now we will show that this phase is just and so it is also independent of the choice of image F . First, consider the case F = [1, N/2]. In this case g F is the identity permutation. Therefore and so (3.34) is satisfied.
Next consider a single permutation which exchanges two elements j, k ∈ [1, N ]. Let j < k. If both j and k are in F or in [1, N ] \ F then F is invariant and so the quantity (3.34) cannot change. For concreteness, let j ∈ F and k ∈ [1, N ] \ F . Then a permutation which exchanges j and k, acting on g from the right, will yield a new image Let j be the number of elements of F which are less than j and let k be the number of elements of F which are less than k. Since there are j − 1 and k − 1 elements which are lower than j and k in [1, N ], then there must be j − j − 1 elements of [1, N ] \ F which are lower than j and k − k − 1 which are lower than j.
Now we are ready to calculate the change in (−1) g F resulting from the exchange of j and k. Note that the resulting permutation g F is not just the transposition (jk). This is because once j and k have been exchanged, the elements of F and [1, N ] \ F are no longer in order. k now appears in F in the old position of j, which has j elements before it. However k has k elements beneath it in F . One of these elements is j, which is not present in F , and so k has k − 1 elements beneath it in F . The result is that in order to order F , k needs to be moved to the right k − j − 1 places, requiring k − j − 1 transpositions. Similarly, j now appears after k − k − 1 elements of [1, N ] \ F but it is only greater than j − j − 1 elements, and so to order [1, N ] \ F requires (k − j − k + j ) transpositions. Summing the transposition (jk) and the transpositions in F and its complement we find The change in the other sign factor is and so Therefore the phase factor is invariant under all interchanges of two elements. But one may obtain any image F via interchanges of pairs of elements, and so the phase factor is the same for all F , therefore it is always equal to its value when F = [1, N/2] which is given in Eq. (3.35). Therefore we have proven (3.34).
Substituting (3.34) into (3.33) we find our master formula for the contribution of one Recall that the sine terms are all positive and the phases are manifestly independent of F , and so we have arrived at

An All Even Image
As we may decompose the symmetric group into cosets we may also decompose the matrix elements into terms of the form (3.40) This is step 4 of our method.
As the contributions ψ F to ψ 0 are all of the same phase, our approach will be to find the largest and then to expand about it. The sine functions are maximized when the distances between the elements of F and of its complement are largest. Therefore a reasonable guess at the largest ψ F is the case when F consists of all even or all odd elements of [1, N ]. Let us introduce some notation. Define the vector So F is a sequence of N 0s and 1s with N/2 of each. Let us name two common subsequences So the all even subset F corresponds to F =↓ · · · ↓ with N/2 ↓'s. Similarly the odd subset F corresponds to N/2 ↑'s. For example, if F = {0, 1, 1, 0} then one could write the shorthand F =↓↑.
In Eq. (3.40) one sees that ψ F only depends upon the differences between the P F A (j). These differences are identical for the even and odd subsets, so they yield the same ψ F . We will now calculate ψ F for the even subset The term on the right looks quite similar to the first expression in (3.4).
In that case we saw Substituting the square of this identity into (3.46) one finds as we fixed the wave function to be real and positive. As all ψ F have the same phase, this phase is inessential. The key point is that ψ ↓···↓ contributes (N/2) N/2 of the total, which is N N/2 . Stated differently, we may expand The exact result is However just the even image ψ ↓···↓ alone contributes Therefore ψ ↓···↓ alone exactly yields the leading zero and nonzero terms, although the subleading term is incorrect. This is no surprise, as the other terms in the sum (3.50), corresponding to other images, have not been included.

Summing over Choices of Image: Arrows
In the rest of this section, we will include more images F in an attempt to improve our estimate of γ. Let us start by replacing the ith ↓ with an ↑. Now How does the expression (3.40) for ψ ↓···↑···↓ differ from ψ ↓···↓ ? In the latter case, the differences in the P F A are all even. Now, when j = i, P A (j) decreases by one unit and so the argument of the first sine decreases by π/N , while that of the second increases. When k = i, the first increases and the second decreases. As a result .

(3.55)
So far our calculation has been exact. Now let us consider an expansion around large N . Choose a natural number k << N . Consider just those terms of the product such that |i − j| ≤ k. If the large N limit is taken before the large k limit, then this product converges in k and we may approximate each sin(x) by x to obtain where the first product excludes j = i. The first approximation, k = 1, is (3.57) The approximation k = 1 is easy to understand. It means that one only considers adjacent arrows. Each pair corresponds to the following term More generally, one may use the k = 1 approximation (3.58) to find ψ F for any F that corresponds to a sequence of arrows, in other words any F that includes precisely one member of each pair {2j − 1, 2j}. The k = 1 approximation is reasonable because k = 1 yields the contribution to (3.56) which is furthest from unity, and so the most important contribution to the absolute value of the matrix element. This is the case because the corresponding |i − j| = 1 term is the most important in (3.55). But there is one exception to this observation. If there is a |i − j| = N/2 − 1 term, it will have the same value. Such a term only appears for the arrow at the beginning and at the end of F. Therefore the k = 1 approximation means that one considers only adjacent pairs of arrows in the N/2-vector F, where adjacency is understood modulo N/2 so that the first and last arrows are also adjacent.
For an arbitrary string of arrows F, let j be the number of occurances of ↓↑, with the above cyclicity condition understood. Then j will also be the number of occurances of ↑↓ and equivalently the number of strings of consecutive ↑'s and the number of strings of consecutive ↓'s. For example ↓↓↑↑↑↑↑↓ corresponds to j = 1. There are thus 2j subsequences ↑↓ or ↓↑, each leading to a factor of 3/4 in the k = 1 approximation according to Eq. (3.58). Therefore How many F 's are there with each j? One needs to place the 2j places where the direction of the arrow changes amongst N/2 arrows. At each of these places one inserts ↑↓ or ↓↑, which is of length 2. Thus there are a total of N/2 slots where one needs to place 2j transitions. The number of such choices is N/2 2j . Note that one can choose between ↑↓ and ↓↑ only once, leading to a single factor of 2 which we will ignore. For example, if the first transition is ↑↓ then the second will necessarily be ↓↑ and so on. The corresponding contribution to ψ 0 is ψ 1 where

Summing over Choices of Image: Beyond Arrows
In general, not every pair {2j − 1, 2j} will have precisely one element in F . It may have 0 or 2 elements in F . Let us include these cases by adding to our notation for pairs in ( What happens to the calculation in Subsec. 3.5 if we include a single adjacent pair 20 or 02 starting at the ith position? Recall that in that subsection we adopted the approximation k = 1 in which we only considered sine terms in (3.40) corresponding to differences between neighboring symbols, understood modulo N/2 so that the first and last symbols are neighbors. Let us refer to the sine terms that we consider as interactions. In that case, there were N/2 symbols ↑ or ↓, each of which had an interaction with the neighbor on each side leading to a total of N/2 interactions in S A and also in S B , for a total of N interactions. Each interaction came with a factor of 1/N , but these factors cancel in ratios ψ F 's because there were the same number of interactions, N . Now it is no longer obvious that the number of interactions is invariant when a 20 replaces two arrows. The two arrows were involved in 3 interactions: one with their neighbor on the left, one with their neighbor on the right and one between the two arrows. When the two arrows are removed, so are these three interactions. Now the 20 has two interactions with the arrow on its left, and one interaction between the two arrows in the 2. So the number of interactions counted is the same when two arrows are replaced by a 20. This would not be the case were the 2 and the 0 separated. Therefore, in the case of 20 and 02 insertions, we can continue to use our k = 1 approximation in which only interactions between neighboring symbols are considered, but one must recall that there is now also an interaction inside of the 2.
So what are the interactions of a 20 with its neighbors? The internal interaction gives the difference between the two elements of the 2. These are elements 2i − 1 and 2i so their difference is one, giving a factor of sin π N for both the S A and the S B terms in (3.40). The interactions with the neighbors depend on the value of the neighbor. If the left neighbor is ↓ then 2i − 2 ∈ F while 2i − 3 ∈ [1, N ] \ F . The former is separated by 1 unit from 2i − 1 and by 2 units from 2i, both of which are in F , and so one obtains a factor of sin π N sin 2π N from the S A terms. If the left neighbor is ↑, one similarly obtains a factor of sin 2π N sin 3π N .
As neither 2i − 1 nor 2i is in [1, N ] \ F , the left neighbor does not contribute to the S B terms. The right hand neighbor similarly has two possibilities which lead to these two factors. In all, we have three possible products of the six neighborly interactions involving the 20 (3.67) These can be compared with the sin 6 2π N that would be obtained for ψ ↑···↑ or the average (7/4) 2 sin 6 2π N obtained for ψ 1 in Subsec. 3.5.
Again let us take the large N limit by approximating sin(x) ∼ x. Then we obtain the ratios Here we have introduced the notationψ 1 to denote the sum over 8 terms ψ ↓···↓xxx↓···↓ where x =↓ or ↑. As a rough approximation, we will simply average these to obtain where ψ ···20··· is the sum over all ψ with F consisting entirely of arrows except for a single 20 insertion at a fixed position. The denominator is ψ 1 and notψ 1 as we have now summed over all combinations of arrows even far from the 20 insertion.
What about multiple insertions of 20? This is a length 4 string. In principle it may be inserted at N/2 different places in the string of N/2 arrows, 0's and 2's. However, many of these positions overlap. To simplify the calculation, let us make the approximation that there are N/4 places to insert it, and simply add a degeneracy factor of 2 reflecting the fact that it could be displaced by 1/2 of a location, ignoring the possible overlaps. Let us add another degeneracy factor of 2 for the fact that 20 may be replaced by 02. And let us ignore the fact that neighboring 20's and 02 s will have different interactions than those treated above, where we assumed that the neighbors were arrows.
With these approximations, the number of ways to insert j strings 20 is N/4 j with a weight of (4/49) j and a degeneracy factor of 4 j . Altogether this yields the sum In Table 1 we compare the k = 1 estimates above for the sums over subsets of the images F with the exact sums over those subsets calculated using (2.15) at N = 6, N = 8 and N = 10. The quantity ψ ↓···↓ was calculated exactly above, for any N , and so it is of no surprise that (3.49), reported in the first row of the first column, is equal to the exact results reported in the first column of the later rows. ψ 1 on the other hand is systematically over estimated by (3.63), reported in the second column of the first row, with respect to the exact results in the second column of the later rows. That is as a result of our k = 1 approximation. The k = 2 correction, for example, multiplies our (3/4) factor by (15/16) yielding 45/64 and so it reduces the expected 7/8 to 109/64. Each k correction in fact is negative. However it is important to remember that they can only be applied to those insertions which are at least a distance k from their neighbors. This correction of course also applies to ψ 2 , reported in the last column of the first row, as it was calculated from ψ 2 /ψ 1 in (3.70). With the k = 2 correction, ψ 2 will be less than ψ 0 , as it must be since it is a partial sum over positive contributions to ψ 0 . In our examples we also see that of course ψ 2 is less than ψ 0 . With more work, the estimates could become more precise. However in the table we see that ψ 2/N 2 already provides an estimate of ψ 2/N 0 with an error of less than 5%.

Summing over Choices of Image: A Systematic Approach
The above system of ever increasing precision by summing over more possible images F can be made systematic. P is a map from [1, N ] to [1, N ]. We always binned the domain into 2 bins, S A and S B . In the first step, in Subsec. 3.4, we choose only the value of F which gave the largest ψ F , corresponding to all even or odd elements of [1, N ]. This gave us ψ ↑···↑ = (N/2) N/2 . Next in Subsec. 3.5 we also binned the range [1, N ], into N/2 bins {2j − 1, 2j} and we only considered those images F which contained precisely one value from each bin. If f Ij (F ) is the number of elements of S I that are mapped into {2j − 1, 2j} then this corresponds to the case with all values equal to unity f Ij (F ) = 1. (3.72) We are trying to sum over all g ∈ S N . We have partitioned S N into cosets H F ⊂ S N each labeled by F . Now we have further partitioned F into bins labeled by the 2 × N/2 matrix f Ij with entries 0, 1 and 2. We considered only the sine terms, which we called interactions, in (3.40) which connect neighboring bins j and j . Therefore in Subsec. 3.5 we have again considered the contributions to ψ 0 arising from the dominant bin, which now is ( Actually one may make a similar interpretation of Subsec. 3.4 , where S is the first row of a 2 × N matrix f Ij which gives the number of elements of S I equal to j, which is necessarily 0 or 1 as j is a single element. In this sense the difference between Subsecs. 3.5 and 3.4 is that the first partitioned the image [1, N ] into N/2 bins while the second partitioned it into N bins, so that each element was its own bin.

Two Densities
We expect that the sum over permutations in the Bethe Ansatz (2.15) is simpler when the m i or equivalently the x i are distributed linearly. The purpose of binning is that while the x i are generally not linearly distributed for a matrix element of interest, by considering sufficiently small bins, the x i may be linear in each bin up to any desired precision.
In Sec. 3 the x i were already distributed linearly, as the density was taken to be constant. The purpose of that section was merely to show how the binning approximation can be implemented systematically. However it did not lead to any simplification, quite on the contrary.
In this section we will therefore consider matrix elements of the Tonks-Girardeau ground state with states with two regions with distinct densities. In Subsec. 4.1 one of these densities will be zero. In Subsec. 4.2, each region will have the same number of particles. In these cases, the Vandermonde determinant form for the ground state is of complexity N 2 , and so in principle requires at least time N 2 to evaluate. We will use the binning approach in a way which is manifestly independent of N , and so will arrive at results which may be a suitable starting place for a large N limit.
Of course in this case it would be possible to bin the Vandermonde product directly, and so bypass the complicated partitioning of the permutation group S N . However beyond the ground state in this model and also in other models treatable with the Bethe Ansatz, there is no Vandermonde product form, and so that method would have no hope of generalizing.

When One Density Vanishes
In this case we will consider particle positions of the form where α is an arbitrary real number. Then (2.11) yields η j = e 2πix j /L = e 2πiαj/N = χ αj/N (4.2) and the Vandermonde determinant of the matrix η k j provides arbitrary matrix elements This expression for the wave function is exact, but to calculate it requires a time polynomial in N and so it is not suitable for a large N limit. To simplify it, we will need to choose a regime for α. The case α = 1 was the subject of Sec. 3.
1 The Case: 0 < α << 1/2 In this case the argument of the sine function is always much less than π/2 and so we may use the approximation sin(x) = x to obtain Note that the leading term in the exponent, which is of order N 2 ln(N ) in the penultimate line of (4.4), cancels in the final expression so long as α is held fixed in the large N limit.
Even so, the norm of the state is of order e N ln(N ) and the volume of the N -dimensional coordinate space is of order L N , so matrix elements of order e −N 2 will have a measure zero contribution to the wave function even if they are all integrated over. Thus states with α << 1/2 do not contribute to the wave function in the large N limit.
In Fig. 1 we compare our approximation (4.4) to the exact Vandermonde formula and find good agreement at large N and α < 0.5. which were proved in Subsec. 3.1. These allow us to use the fact that α is near 1 to cleanly separate the N 2 terms in the exponent.
At large N , without loss of generality we can take N/α to be an integer. Then where A and B are defined to be the two factors in the penultimate term. These can be evaluated separately where the two factors in the first line were evaluated using the identities (4.6) and (4.5) respectively, with M = N/α.
This is similar to ψ 0 in Subsec. 1 but with so that 0 < α << 1/2. With these substitutions, Eq. (4.4) yields (4.11) Substituting A and B into Eq. (4.7) we find the final result As expected, when α = 1 at finite N , this expression equals (3.19). The limit of interest of course is α −→ 1, N −→ ∞. In this case, it differs from (3.19) when the first term differs from unity, which is the case in which N 2 (1 − α) 2 ln(1 − α) does not tend to 0.
To see this, let us define = 1 − α. (4.13) The leading order terms in (4.12) as −→ 0 are (4.14) The leading term is just (3.19). The corrections will tend to 0 in the large N limit if N → 0, in other words if shrinks faster than 1/N . If N tends to a nonzero constant, then both constants are of the same order. Otherwise, N → ∞, in which case the N 2 term dominates, corresponding to the first term in (4.12). Which case is relevant?
Recall that the normalization of our energy eigenstates is of order N ! and so its logarithm is of order N ln(N ). Also, if x is discretized then the norms of the position states will have logarithms of order N . Therefore we expect contributions from of order e N neighboring states, and terms in the logarithm of the matrix element of order N and N ln(N ) may contribute to finite quantities. For higher order terms however, no amount of normalization and integration can avoid the fact that one finds zero in the large N limit. This argument suggests that the dominant contributions to deformations will come from those that correct ln (ψ 0 (x)) with corrections of order N or N ln(N ), and so correspond to In this case it is the first term in (4.12) which dominates. This is fortunate, as we will soon see that the second is the most difficult to estimate in general.
In Fig. 2 we compare our approximation (4.12) to the exact Vandermonde formula and again find good agreement at large N and 0.5 < α < 1. which lead to Now the families of matrix elements which contribute to the wave function are neighborhoods of the cases in which α − 1 shrinks at least as quickly as N −1/2 .
In Fig. 3 we compare our approximation (4.18) to the exact Vandermonde formula and again find reasonable agreement at large N and 1 < α < 1.3. The exact formula produces a series of sharp dips, corresponding to the values of α at which αj = N for some j < N . Physically these are points where two particles coincide, which have zero wave function due to the impenetrability condition (2.3). The Vandermonde formula has zeroes at These zeroes are not captured by our approximation. At large N the zeroes become thin, but their density increases as N . We do not know whether they persist in the large N limit, or whether they lead to a physically relevant drop in the wave function at α > 1 in this limit. However we will see below that in the two bin case, this regime does lead us to systematically overestimate matrix elements.

Equal Number of Particles at Each Density
Finally we are ready for a nontrivial application of the binning formalism. We will consider an equal number of particles at density α and 2 − α in two contiguous regions (4.20) Recall that our Bethe Ansatz (2.15) describes the wave function in terms of a map P : [1, N ] → [1, N ] where the domain corresponds to our x's and the range to the momenta, which are integers. We will again consider the two bins S A and S B defined in Eq. (3.22). Our binning approach to the constant density case in Sec. 3 suggests that we can systematically improve our estimates of the leading order term in ln(ψ 0 ) by considering more and more images of S A . Our analysis of the constant density case in Subsec. 4.1 suggests that now this leading order will consist of the N 2 terms in ln(ψ 0 ).
We now face new complications with respect to Sec. 3. First, as the distance between some pairs of x is greater than L/N , the argument of some of the sine functions in (3.40) will be negative. Second, there is no symmetry condition which determines which set F should lead to the dominant contribution to the matrix elements.
As the x j in each bin have constant separation, the results of Subsec. 3.3 largely carry over to this case. The symmetric group S N can again be divided into cosets labeled by F , defined such that P (S A ) = F for all g in the coset. Then sum of (2.15) over the coset yields a formula similar to Eq. (3.24) The average value of x in S A still differs from that of S B by L/2. Therefore again, up to an overall phase, ψ F can be written as a product ψ A ψ B .
As the x j are linear in each bin, these can again be written as Vandermonde determinant and so as a product of sines, now with factors of α and 2 − α respectively. Recombining them we obtain the analogue of Eq. (3.40) (4.22) Again the phases are F -independent except for the sines themselves, which are now negative if their arguments are greater than π. Therefore the ψ F will all contribute to ψ 0 with the same phase up to a sign. In principle we could calculate the distribution of these signs to refine our estimate of ψ 0 , however we will leave this to future work.
The zeros and sign changes of the sine arise identically to those in Eq. (4.19) in the case of a constant density. In that case they had a very physical origin, they were the zeroes of the wave function ψ 0 when two particles overlap. Now the particle positions are given by (4.20) and no particles overlap, therefore the total wave function ψ 0 does not vanish.
Instead the zeroes in (4.22) are zeroes in ψ F , corresponding to the sum over a single image F . These zeroes come from the fact that the wave function would have vanished had the same x spacing α or 2 − α persisted over all of [1, N ] instead of just the bin S A or S B . More concretely, they arise from the fact that the wave function Slater determinant contains all behaviors from e −i(N −1)x/2L to e i(N −1)x/2L and so if two adjacent x's differ by L/(N − 1), for example, corresponding to α = N/(N − 1), then a minor in the determinant which includes e −i(N −1)x/2L for one x j and e i(N −1)x/2L for its neighbor will vanish. In this sense, the zeroes of (4.22) are zeroes in a minor and not in the determinant so they are spurious and it would be desirable if (2.15) could be resummed in such a way that they are not present from the beginning.
Our goal will now be to make the crudest possible estimate of ψ 0 , as a ψ F which is in a sense dominant. What do we mean by dominant? Different definitions give different estimates. One choice of image is just F =↓ · · · ↓ as before. In Fig. 4 this choice is compared with the exact ψ 0 calculated using (2.11).

Optimizing the Image
How does our estimate of ψ 0 depend upon our choice of F . Recall that F is a subset of [1, N ] consisting of N/2 elements. We can parameterize this interval by and in the large N limit we can associate a local density Clearly F determines ψ F . And F determines ρ(F ). But to what extent does ρ(F ) determine ψ F ? We claim that ρ(F ) fixes the dominant N 2 term in ln(ψ 0 ) and also perhaps the N ln(N ) term. We cannot prove this claim, but we will provide an example.
Consider two choices of image F . One is F 1 = F ↓···↓ which is the subset of even elements. The second, F 2 = F ↓↑↓···↑ alternates between ↑ and ↓, corresponding to F = {011001100 · · · 01}. For simplicity we will set α = 1, so that α does not dominate over the effect of this distinction. We do not expect that this choice will affect our conclusions. In both cases the separations between elements of F and [1, N ]\F are equal, and so ψ A = ψ B .
In the first case and in the second In this case the change in local density only affects the O(N ) term in ln(ψ).
Thus the leading behavior of ln(ψ) appears to be entirely determined by ρ(x). How do we optimize ρ(x)? Let us take ψ A and ψ B to be real and positive. We will now allow α to be an arbitrary real number. In this case the arguments of the sine terms in (4.25) and (4.26) do not necessarily lie in the interval [0, π]. However the absolute value on the left hand side of these equations means that each sine term should be |sin|. Then where c A and c B are defined by the preceding expressions. In the continuum limit the sums are replaced by integrals over x and y and the set F is replaced by its density ρ(x) c A ∼ where is a small number that cuts off the divergence as x → y, reflecting the fact that j = k in the original discrete expression (4.27). The expressions c A and c B are those for the energy of a system of density ρ with an interparticle potential given by ln(sin).
In this analogy, the potential V (z) for a particle at x is the functional derivative with respect to ρ(z) of c A + c B . However F has a fixed number of particles, so our extremization The optimal values of ρ found by solving (4.30) discretized into 25 parts J. The red, green, blue, black and brown curves correspond to α equal to 0.2, 0.4, 0.6, 0.8 and 1 respectively. Note that when α = 1 the optimal density is constant as expected.
condition must be that c A + c B is invariant under moving a particle from z 1 to z 2 Therefore one finds that extremization requires that the potential must be z-independent, or in other words that the force F (z) vanishes 0 = F (z) = ∂ ∂z V (z) = −π 1 0 dx [ρ(x) (αcot (πα(x − z)) + (2 − α)cot (π(2 − α)(x − z))) +(α − 2)cot (π(2 − α)(x − z))]. (4.30) We solve this equation for ρ by discretizing, and reducing it to a linear matrix equation. The equation is solved by inverting the operator which multiplies ρ and multiplying it by the term on the last line. The dimension of this matrix is independent of N , and so N has disappeared from the problem. We discretized into 25 parts, so that ρ(x) is defined by at 25 positions. Then we solved the equation to obtain our preferred value of ρ(x), displayed in Fig. 5. At various values of N , we then created a subset F ∈ [1, N ] with local density ρ.
The results are reported in red in Fig. 4. As a result of this clumsy optimization procedure, our results are quite noisy. In general at high α it systematically underestimates ψ 0 . However at low values of α it is closer to ψ 0 than ψ ↓···↓ .

Conclusions
State of the art computations of coordinate matrix elements in integrable models [9,10] are only able to compute matrix elements with respect to coordinate states that are in some sense homogenous or linear, such as Néel states. Our goal is to generalize these computations to states which are piecewise linear. The strategy that we adopt is the decomposition of the lattice sites or particle positions into piecewise linear bins, which induces a decomposition of sum in the Coordinate Bethe Ansatz.
More precisely, the CBA computes matrix elements as a sum of N ! terms, corresponding to maps P : [1, N ] −→ [1, N ] (5.1) or equivalently to elements of the symmetric group S N . Summing over all N ! terms in general is difficult. But we found that if the domain [1, N ] is divided into bins S i and one fixes the images P (S i ), the corresponding subset of terms can be summed exactly. Each such subset is a coset H ⊂ S N and S N is the disjoint union of these cosets. We provided this sum in Eq. (3.40). Thus the task of summing over all N ! group elements is reduced to the task of summing over the cosets, of which there are only of order 2 N .
More critically, in the case treated in Sec. 3 in which the density was constant, the sums over all cosets had the same phase. Thus one needed only to sum over the dominant cosets to obtain a good approximation. More generally, in Sec. 4 the sums over cosets had the same phase up to a sign. Our approximations to the matrix elements were generally too high because we did not consider the fluctuation in signs.
However for densities near the mean density, which are those that will contribute to wave functions in the large N limit, the fluctuations are caused only by the terms near the boundaries of the image [1, N ]. We thus feel that it should be possible to develop an approximation scheme, along the lines of that developed in Sec. 3 for transitions between ↑ and ↓, which allows us to estimate the correction caused by these sign flips. This will be done in future work.
Another critical generalization will be to other Hamiltonian eigenstates. Our motivation, understanding mass gaps in quantum field theory, suggests that we only need to investigate low lying states. For these, only a small number of momenta are changed and so we expect that we only need to modify our approximation scheme to keep track of which bin contains these excitations or, more simply, one may separately sum over P −1 of the excited momenta. The more challenging step comes next, we will try to generalize these results to the Lieb-Liniger model. As the Tonks-Girardeau model is the strong coupling limit of the Lieb-Liniger model, we will be able to test our results. If this can be achieved, we expect that spin chains such as the Heisenberg XXX model should present no new conceptual difficulties, at least at s = 1/2.