Composite-boson formalism applied to strongly bound fermion pairs in a one-dimensional trap

We analyze a system of fermions in a one dimensional harmonic trap with attractive delta-interactions between different fermions species, as an approximate description of experiments involving atomic dimers. We solve the problem of two fermion pairs numerically using the so-called “coboson formalism” as an alternative to techniques which are based on the single-particle basis. This allows us to explore the strongly bound regime, approaching the limit of infinite attraction in which the composite particles behave as hard-core bosons. Our procedure is computationally inexpensive and illustrates how the coboson toolbox is useful for ultracold atom systems even in absence of condensation

The possibility to engineer atomic and molecular many-body systems by controlling and assembling simpler components has made enormous progress thanks to Feshbach resonances.In this way, molecular Bose-Einstein condensates have been formed starting from ultracold atomic gases [1,2].Similar setups have been used for the controlled observation of relevant phenomena in statistical physics such as Wigner crystals [3], and the BEC-BCS crossover [4,5].Within the field of ultracold Fermi gases, one-dimensional systems are known to exhibit very peculiar properties [6].In particular, strongly bound fermion pairs reach a limit in which they behave as hard-core bosons, which in turn can be treated through fermionization [7].
We consider a one-dimensional scenario, with fermions of two different kinds in a harmonic trap and an attractive contact interaction leading to fermion pairing.Much effort has been devoted to this kind of system [8][9][10][11][12], but the usual numerical treatment takes as a basis the harmonic oscillator eigenstates, making computations very costly for strong attraction.Alternative procedures which are more efficient for strong attraction have been proposed in [13,14].Here, as a different approach, we tackle the problem of two pairs with two fermions each in the context of coboson theory [15,16].This theoretical framework, originally developed for excitons in semiconductors [15][16][17], has by now been applied to a variety of systems, including Bose-Einstein condensates [18], superconductors [19,20] and Feshbach molecules [21].
A very useful simplification often encountered in this treatment is the so-called coboson ansatz, which is analogous to a condensate formed by composite bosons [15,20].Using tools from the coboson formalism, we show that the coboson ansatz does not provide a good approximation of the true ground state for the case of two pairs in the limit of strong interaction.This is to be expected in the light of previous results [22,23] and also because the limit of infinitely bound pairs corresponds to hard-core bosons which are known to form only a quasi-condensate in 1D traps [24][25][26].However, the coboson formalism also provides tools to describe the state beyond the coboson ansatz [19,20].We thus develop a representation of the problem in the coboson basis, i.e. in terms of the eigenstates of one pair of interacting fermions in the trap.This is specially convenient and is expected to work better for the regime of strong attraction.In this respect, our method is related with the perturbative approach in [27].
The basic steps of the procedure are the same as in [23] and are as follows: 1. We solve the problem of a pair of interacting fermions in the trap.The operators B † n that create each one-pair eigenstate, and the corresponding energies E n , will be the starting point of the treatment.We truncate the basis considering the states with the lowest energies, up to some quantum number n max .
2. From the single-pair basis operators B † n we form the two-coboson basis generated by the action on the vacuum of operators of the kind B † n B † m .
3. We calculate the form of the Hamiltonian in this truncated coboson basis.
4. Solving the corresponding generalized eigenvalue problem, we estimate the ground state for two pairs and analyze its properties.
This method allows us to interpolate from the interaction strengths for which the singleparticle basis is suitable [9][10][11][12], all the way to very strongly bound pairs approaching the limit of hard-core bosons.Using coboson-theory tools combined with Taylor expansions, we calculate several quantities of interest, including the energy and two-particle correlators.The work is presented as follows: in Sec. 2 we review how to write the problem in the coboson framework.Section 3 is devoted to analytical considerations for infinite attraction.In Sec. 4 we discuss our numerical results.A summary and conclusions are given in Sec. 5. Finally, several appendices with detailed calculations are included.
2 The procedure, step by step

Single-pair solution
For definiteness we will assume that both fermion kinds, which we call a and b, have the same mass, and that the creation and annihilation operators corresponding to different fermion species commute (this last choice does not affect the final results).We also assume that the trapping potential is the same for both species.
The first step requires the solution of the single-pair problem, with a Hamiltonian given by: This problem can be solved by separation of the center-of-mass and relative variables.The center-of-mass solution is given by the harmonic oscillator eigenfunctions corresponding to mass 2m.The relative motion has been solved in the general case in Refs.[28,29] but for simplicity we focus only on strongly bound pairs, so that the relative motion has a wavefunction of the form of an exponential, and the energy associated with the relative motion can be approximated by: In this regime, the single-pair eigenfunctions are then approximately of the form: where ϕ n are the harmonic oscillator eigenfunctions for a particle of mass 2m.The corresponding energies are: From these solutions, we define the coboson creation operators B † n such that: where |ñ is the n-th single-pair eigenstate, and |v is the vacuum.In particular, the coboson operators B † n can be written in terms of field operators as: For consistency, neglecting states where the internal motion is excited implies also a truncation in the center-of-mass states.In particular, we keep only states where the index n is such that the excited internal states are well above the energy scales considered, i.e.: For convenience here we have defined a spatial scale x ω associated with the harmonic oscillator, Once more, this stresses the fact that this restricted basis is only appropriate for strong attraction, when the size of each bound pair is very small compared with the spatial scale of the trap.

Basis for two pairs
From the set of states corresponding to the lowest energies of the single-pair Hamiltonian, one can form states of the form: with n ≤ m (we note that the coboson creation operators commute) and |v the vacuum.Because of the fermionic character of the constituent particles, states generated in this form are neither normalized nor orthogonal [15].We truncate this two-pair basis with the condition n + m ≤ n max .An often useful approximation for the ground state of dilute systems of N pairs with short-range interactions is given by what we call the "coboson ansatz" [15].This corresponds to the state obtained from the repeated application on the vacuum of the operator B 0 that creates a single pair in its ground state: where χ N is a normalization constant.However, this can only provide a good approximation of the true ground state in systems which are expected to exhibit condensation at zero temperature This is not the case in the problem we analyze [22,23,25,26].In order to quantify the quality of the approximation, we study the fidelity F between the true ground state for two pairs, |GS , and the coboson ansatz: Even if the coboson ansatz is not a good approximation, one can still compute the ground state by means of the coboson formalism.In order to do this, we will work with the space generated by the coboson operators as in Eq. (10).First, we compute all overlaps between the relevant states from the expression: ) Here X α with α = a, b is an operator that exchanges the states of the two fermions of kind α, and it acts on a fictitious space where fermions of equal kind are treated as distinguishable.Since our goal is to find the ground state, instead of building an orthonormal basis, we keep the overlap matrix S to solve the corresponding generalized eigenvalue problem.
The matrix S can be calculated following different strategies.In the coboson literature [15], the overlaps are evaluated in terms of matrix elements of the change of basis between single-pair eigenstates and the separable single-fermion basis.However, this procedure can be numerically costly and lead to large errors when many coefficients are non-negligible and no analytical expression exists for the sums required.Thus, we resort to a different form of evaluation.Plugging the explicit form of the operators B † n given by Eq. ( 7) in all formulas, and using (anti)commutators, we can obtain an expression for the elements of the overlap matrix as: Since we are interested in the case of strong attraction, the factors of the form e −λ|y j | allow us to perform a Taylor expansion in 1/(x ω λ) for the harmonic oscillator functions ϕ n .This is possible given the truncation of our basis in Eq. ( 8), which implies that the spatial scale associated with the center of mass is much longer than the pair size λ −1 .In this form one can find approximate expressions for S from a lengthy but straightforward evaluation of spatial integrals.This procedure is explained in detail in Appendix A.

Construction of the Hamiltonian
We now need to compute the Hamiltonian in the coboson basis.The Hamiltonian can be split in two parts, corresponding to the non-interacting terms and the interactions.The first part, which contains the kinetic energy and trap potential, is the one-body term.It will not play a role in the following calculation (apart obviously from appearing through quantities that were calculated when solving the single-pair case).The interaction part is quarctic and can be written in terms of field operators as: The Hamiltonian matrix elements in the coboson basis can be obtained from the expression: which is just a rewriting of the formulas in [15].In a similar spirit as for the calculation of the overlap matrix S, instead of following the standard expressions in [15] we estimate the Hamitonian elements using a Taylor expansion of spatial integrals.
In particular, the last line of Eq. ( 16) can be written as: The details of the procedure involving the Taylor expansion of the Hamiltonian elements are also provided in Appendix A.

Analytical considerations for infinite attraction
Before presenting the results of our numerical approach, we note that the case of infinite attraction can be solved exactly.In this limit, fermions of different species are so strongly bound that they behave as point-like hard-core bosons of mass 2m, and the problem can be solved by means of fermionization [7].The wavefunction of these hard-core bosons is thus given by: where the subindex "hc" stands for "hard-core".
In this form, we obtain for λ → ∞ an asymptotic interaction energy ∆E = hω.Here, the interaction energy is defined as ∆E = E 2 − 2E 1 , where E M is the ground-state energy of M = 1, 2 pairs.We also find an asymptotic fidelity of F ∞ = 2/π 0.64 between the coboson ansatz and the true ground state.This fidelity is lower than the one obtained in the same regime for two fermion pairs in translationally invariant models [22,23].
Following the same lines, one can find the joint density of composite particles at positions x and x for the limit of infinite attraction.This is of the form: One can also write down the conditional probability P(x |x) of finding a composite pointlike particle at position x provided that another one was found at position x: Furthermore, one can calculate the asymptotic values of the coefficients in the expansion of the ground state in the coboson basis according to: and one obtains for λ → ∞: This expression is valid for even and nonzero n + m, and here δ mn is the Kronecker delta.Due to symmetry reasons the coefficients c mn vanish for odd n + m, and for n = m = 0 we find: Since the coboson ansatz corresponds to the repeated application of the coboson operator B † 0 , and for λ → ∞ the wavefunctions associated with the different B † m become orthogonal, the asymptotic value of c 00 determines the asymptotic fidelity between the correct ground state and the coboson ansatz.The additional factor √ 2 in the fidelity comes from the definition of the coboson basis in Eq. ( 10), which does not include a prefactor 1/ √ 2 for m = n.
Before tackling the numerical treatment of the problem for strong but finite attraction, we note that also the limit of infinitesimal attraction can be treated analytically.For g = 0, the ground state of the system is separable, with the two lowest oscillator levels occupied for both kinds of fermions.Then, the energy ∆E approaches 2hω.It is very important to notice that in this separable limit, the coboson normalization factor χ 2 in Eq.( 11) vanishes, and thus the coboson ansatz is not defined for g = 0. Nevertheless, using perturbation theory together with analytical results for the Schmidt coefficients [30] one can calculate the limit value of the fidelity between the true ground state and the coboson ansatz, and find that as the attractive interaction strength approaches zero, F approaches a value of approximately 0.37.Indeed, for g ∼ 0 we obtain χ 2 ∼ 0.342 λ 2 and F ∼ λ 2 /8χ 2 with λ ∼ g/ √ 2πhωx ω .We note, however, that the weakly bound case is not within the scope of our present study, and it has been extensively analyzed before [9][10][11][12].

Numerical study of the ground state for strong attraction
In the following we perform a numerical study of the ground state according to the procedure outlined in Sec. 2. A delicate point in the calculation is the choice of the number of basis states.A very small number leads to a poor description of the system, whereas for a very large number it becomes unjustified to leave out the excited states of the relative motion, and it can also lead to numerical problems if the overlap matrix becomes worse conditioned.As a compromise, we choose the maximum center-of-mass energy included in our description to grow linearly with λ.
In Fig. 1 (a) we show our results for the interaction energy ∆E = E 2 − 2E 1 using a Taylor expansion for the calculation of both the overlap and the Hamiltonian matrices.We also plot in Fig. 1 (b) the fidelity F between the ansatz and the true ground state as a function of λ when choosing the energy in the truncated basis to be given by n max = λx ω .We note that our results show reasonable agreement with the known behaviour for infinite attraction.
As can be seen in the comparison provided in Fig. 2, for λx ω = 30 the coefficients are already very close to the ones obtained from the hard-core boson limit given in Eqs.(22)(23).This also hints at a procedure to perform approximate computations more efficiently: instead of taking the full basis as in Fig. 2, one can use a truncation inspired by the asymptotic values of the coefficients in Eqs.(22)(23).One can also directly approximate the state by taking the coboson basis in Eq. ( 10) to be a function of λ but the coefficients in this basis to be given by the asymptotic values, which gives a fast and compact approximation for the ground state.Indeed, the ground state found numerically for λx ω = 30 has a fidelity of 0.993 with the state obtained taking the asymptotic values of the coefficients and truncating the basis in the same form.From the numerical solution of the problem one can characterize the ground state through several key properties.In particular, in Fig. 3 (a) we illustrate the spatial correlations between fermions of equal kind through the joint density distribution This plot exhibits a strong diagonal correlation corresponding to particles that form a bound pair, with additional much broader peaks corresponding to particles belonging to different pairs.The calculation of D ab is explained in Appendix C. Another quantity that reflects the spatial correlations present in the ground state is the conditional probability P aa (x|0) to find one fermion of kind a at position x given that another identical fermion was found at the origin.This function is plotted in Fig. 4 (a), for the numerical solution with λx ω = 30 and N = 2.For comparison we also show the conditional probability P aa (x|0) obtained from the hard-core limit of λ → ∞ and from the coboson ansatz of Eq. ( 11) evaluated for λx ω = 30.The corresponding formulas are given in Appendix B. The plots show qualitative agreement between the numerical results and the point-like hard-core boson limit, in sharp contrast with the coboson ansatz in its standard form.Indeed, the form of the conditional probability P aa (x|0) is similar to the probability distribution corresponding to the first excited state of the harmonic oscillator, the maxima of which are indicated with dotted vertical lines in the Figure .In a similar manner one can compare the predictions for the spatial correlations of fermions of different kinds.To this aim, we consider the behaviour of the conditional particle density D ab (x|x ) indicating the density of fermions of kind a at position x conditioned on having found a fermion of kind b at position x .We plot this quantity with x = 0 for the numerical solution corresponding to λx ω = 30 in Fig. 4 (b), where we also plot the predictions of the point-like hard-core boson limit and the coboson ansatz for λx ω = 30.The derivation of the corresponding formulas is shown in Appendix C. All three curves have a narrow peak around the origin, associated with the probability to find a fermion paired with the first one detected (in the limit λ → ∞ this peak is a delta function).The curves however differ strongly in the behaviour related with the probability to find the remaining particle of kind a.This second contribution to the conditional density has the same shape as P aa (x|0), and closely resembles the probability distribution for the first excited state of the harmonic oscillator of mass 2m, a behaviour which is not properly described by the standard coboson ansatz.
Figures 3 and 4 were concerned with density distributions in space, associated with diagonal terms of the system's density matrix in space representation.Figure 5 a) shows in contrast an off-diagonal feature, namely the off-diagonal correlation function [9]: where ρ ab is the reduced density matrix for two fermions of different kind.The quantity g 2 is an indicator of two-particle coherence, and indeed the coboson ansatz predicts a constant value g 2 (x) = 1.The numerical results (in black) show that this coherence decays within the typical scale set by the harmonic oscillator, but it stays high for all values with non-negligible particle densities.Nevertheless, the off-diagonal correlation we find is always smaller than the one corresponding to the hard-core limit, depicted in red for comparison.This is not due to a variation in the decay of the spatial coherence, as can be seen in Fig. 5 b).Rather, the difference between our numerical results and the limit λ → ∞ is given by a different density profile, since the particle density at the origin is lower for finite λ than in the limit of infinite attraction.
For the same numerically found ground state one can also characterize the properties in momentum space using similar techniques.In Fig. 6 we show the joint probability distribution for fermions of different kinds in momentum space.This plot displays a strong anti-diagonal peak which is the counterpart of the diagonal peak found for the joint probability distribution in position space, shown in Fig. 3 (b).The remaining features of the plot do not ressemble the state of two identical trapped fermions of mass 2m; this difference in the behaviour of position and momentum is typical of hard-core bosons [7,24,31].The calculation of the joint density in momentum space is similar to the one of D ab (x, x ), but involves a Fourier transform of the coboson basis.The details are explained in Appendix E.

Summary and conclusions
We have tackled the problem of two identical composite particles, each made of two distinguishable fermions, inside a harmonic trap and with contact attractive interactions between fermions of different species.We explored the strongly bound regime using the coboson formalism to build a compact basis of states, greatly reducing the computational requirements associated with the usual description in terms of single-particle eigenstates.We have studied the approach of the interaction energy to the limit of infinite attraction, corresponding to point-like hard-core bosons, and we have confirmed that the coboson ansatz in its standard form does not provide an accurate description of the ground state for any of the interaction strengths within our analysis.We have also shown that the point-like hard-core boson limit provides a good approximation of the coefficients when writing the ground state in the coboson basis.Furthermore, we have used the numerical results to characterize spatial correlations present in the ground state, both between fermions of different and equal kinds, complementing previous work [9].
The composite-boson procedure presented can be generalized to higher numbers of particles and different forms of the trapping potential.Most importantly, we expect this approach to provide an additional tool to the ones usually applied for the description of experiments involving bosonic Feshbach molecules made of fermionic constituents in quasi one-dimensional settings.

A Calculation of Hamiltonian and overlap matrix in position basis
In the limit of very strong interaction, it makes sense to use that the wavefunctions for the center of mass vary over a scale which is much larger than the one for the relative motion.Thus, we start from Eq. ( 14) for the elements of the overlap matrix, use that all y j are of the order of λ −1 , and perform a Taylor expansion in these small displacements.The lowest orders give: Here, all functions are evaluated at position x, the primes mean that a first derivative must be taken, and I mn,jk is an integral of a product of four single-particle harmonic-oscillator eigenstates: These integrals are evaluated using known properties of the Hermite polynomials.In turn, the integrals with derivatives of the eigenfunctions can be written in terms of the elements I mn,jk using the relation: keeping in mind that the ϕ are defined as the eigenfunctions of the harmonic oscillator with mass 2m.
In the same way one can write an expression for the part of the Hamiltonian involving the commutator, Eq. (17).The dominant contributions give, after some manipulations: where again all functions are evaluated at position x and the double primes mean that a second derivative must be taken.This formula can be calculated using similar steps as before.Putting this together with the part from (E j + E k )S mn,jk we can find a consistent expansion for the Hamiltonian up to this order.
For our numerical calculations, we include the orders reported for S and H.One could improve this evaluation by considering higher orders of the Taylor expansion.However, we checked that for the parameter regimes studied the results obtained with these formulas are not significantly altered by excluding the higher order, as can be seen in Fig. 1.

B Spatial correlations for two fermions of equal kind
We first consider the joint density distribution for fermions of equal kind: of course, taking two fermions of kind b leads in our model to the same result.We note that this definition means that: We now show how we calculate this joint density for the numerically found ground state.
Starting from the expansion of the state in the coboson basis, Eq. ( 21), we find: + same with n ↔ m, j ↔ l, and {n, l} ↔ {m, j}.(32) The result for the standard coboson ansatz corresponds to setting all c jl to zero except for c 00 .In the formula above we have introduced auxiliary integrals given by: and The integrals J 2 account for fermion-exchange terms and are negligible unless x and x are close together within a distance of order 1/λ.From these formulas one can recover the vanishing of the conditional probability for x = x for arbitrary states.For the limit λ → ∞, the joint density D aa tends to the expression given in Eq. ( 19) which was calculated from the ground state of two point-like hard-core bosons.In the limit of very large but finite attraction, one can resort to a Taylor expansion for the calculation of the integrals, in the same spirit as the calculations in Appendix A. For J 1 we obtain: (3) j ϕ m + 4ϕ j ϕ (3)  m + ϕ (4) Here, all functions are evaluated at position x and we are using primes (double primes) over the functions to denote derivatives (second derivatives), whereas derivatives of higher order are indicated with superindices between parenthesis.We remind the reader that the ϕ n indicate the oscillator eigenfunctions for mass 2m.
The integral for J 2 can be expanded as: From the joint density D aa (x, x ) one can also calculate the conditional probability P aa (x|x ) of finding a particle of kind a at position x when another of the same kind was found at position x .This can be computed from: so that: For the limit λ → ∞, the conditional probability P aa tends to the expression given in Eq. ( 20) calculated from the ground state of two point-like hard-core bosons.
On the other hand, the standard coboson ansatz predicts for λ → ∞ a behaviour of the form: so that

C Spatial correlations for fermions of different kinds
We now calculate spatial correlations between fermions of different kinds.In particular, we are interested in the joint particle density Here, ρ ab is the reduced density matrix of two different fermions in position basis and is given by [9]: In the following we proceed to the calculation of the joint density for the general numerical solution.Replacing the expansion of the state in the coboson basis leads to: + same with n ↔ m, j ↔ l, and {n, l} ↔ {m, j} .(43) Here, the J 1 are given in Eq. (33), and the J 3 contain interference terms given by: Again, the result for the coboson ansatz is found setting all coefficients c jl to zero except for c 00 .
Resorting to the Taylor expansion J (jl|n) 3 (x, x ) can be approximated by We note that just as in P aa , the terms with J 1 are the only ones that are non-negligible when x and x are at a distance much larger than 1/λ.Thus, P aa and D ab behave in the same way for e −λ|x−x | 1, corresponding to detection of particles in different bound pairs.In the opposite limit of x close to x , D ab has a peak of width 1/λ corresponding to detection of the particle forming a pair with the fermion detected at x .
We now calculate a quantity analogue to P aa (x|x ) but applying to fermions of different kinds.In particular, we wish to calculate the conditional probability P ab (x|x ) of finding a fermion of kind a at position x conditioned on having found a fermion of kind b at position x .This, however, is trickier because after the detection of one fermion of kind b there are two remaining identical fermions of kind a.
Thus, we choose to work with a conditional particle density D ab (x|x ) indicating the density of fermions of kind a at position x conditioned on having found a fermion of kind b at position x .Since two identical fermions can never be found at the same place, this quantity is related with the conditional probability P ab , but its interpretation is more straightforward and, in contrast with a probability, D ab must be normalized to 2. More precisely, the conditional particle density is given by: For infinite attraction, it holds that D ab (x|x ) = P aa (x|x ) + δ(x − x ).For the coboson ansatz, in the limit λ → ∞ one has D ab (x, x ) = 2ϕ 0 (x) 2 [ϕ 0 (x ) 2 +δ(x−x )] and accordingly

D Off-diagonal correlation parameter
Here we provide the expression for the off-diagonal correlation parameter g 2 (x) defined in Eq. ( 26).The diagonal matrix elements appearing in the denominator are particular Here, the asterisk denotes a complex conjugation and we have defined: and J (jl|n) 3 (k, k ) = dqdq ψ j (k, q ) ψ l (q, k ) ψ * n (q, q ) .(55) Replacing the form of the wavefunctions in momentum space one finds the integral expression: Taking into account the restriction on the values of n, l within our basis, one can perform a Taylor expansion in k /λ in the expression above.We stress that the values of k cannot be assumed to be much smaller than λ, since λ is indeed the typical scale for the relative momentum.In this way we obtain: ω πλ e iπ(n−l)/2 dk ϕ n (x 2 ω k )ϕ l (x 2 ω k ) which can be evaluated using properties of the Hermite polynomials.
(58) Performing a Taylor expansion here is justified for the q, q divided by λ in the denominator, but not for the same variables inside the wavefunction ϕ.This makes this calculation much more involved.The Taylor expansion of the denominator gives: e −iπ(j+l−n)/2 dqdq ϕ j (x 2 ω q )ϕ l (x 2 ω q)ϕ n [x 2 ω (q + q − k − k )] Here one is still left with a non-trivial integral in q, q .This can be solved using the decomposition formula A ij|n ϕ i (x) ϕ j (y) , where the coefficients A ij|n = ϕ n (x + y) ϕ i (x) ϕ j (y) dxdy can be evaluated using properties of the Hermite polynomials.For numerical evaluation this summation must be truncated.Performing this one up to i, j = 100 good approximations are obtained.Then we are left with terms similar to those found in the calculation for J 1 .The first term in Eq. ( 53) contains contributions of order x ω /λ which decay in a scale of order λ for the relative momentum (k − k )/2 and of order 1/x ω for the center-of-mass momentum k + k .The second term, involving J 1 , contains contributions of order 1/λ 2 which decay on a scale of order λ for k and k separately.We note that this contribution is broad and has a Lorentzian decay, whereas the decay of the contributions in the first term is Gaussian for the center of mass.Thus, they may be of the same order depending on the point where they are evaluated.In any case, the dominant feature is the anti-diagonal resulting from the first term.The remaining terms, containing J 3 , have a similar behaviour as the first (i.e. with a strong anti-diagonal) but are one order smaller in 1/(λx ω ), which justifies using an expansion for J 3 to a lower order than for J 1 .

5 3considerations for infinite attraction 6 4 7 5and conclusions 11 A 13 B 13 C
Analytical Numerical study of the ground state for strong attraction Summary Calculation of Hamiltonian and overlap matrix in position basis Spatial correlations for two fermions of equal kind Spatial correlations for fermions of different kinds 15 D Off-diagonal correlation parameter 16 E Correlations in momentum space 17 1 Introduction

Figure 1 :
Figure 1: a) Energy for two pairs, excluding the trivial contribution equal to twice the single-pair energy, as a function of λ. b) Fidelity between the coboson ansatz and the numerically found ground state as a function of λ.The results are obtained from the lowest non-trivial order of the Taylor expansion (green stars) and the next non-zero higher-order corrections (black circles) as reported in Appendix A. The horizontal dashed lines indicate the asymptotic values for λ → ∞.

Figure 2 :
Figure 2: Coefficients in the coboson decomposition from the numerical resolution of the problem based on a Taylor expansion for λx ω = 30, in black circles.For comparison we show the values according to the asymptotic expression in Eqs.(22-23) as red four-pointed stars, which overlap with the numerical resuls within the size of the symbols.The basis was truncated with n max = λx ω .The coefficients in the plot are normalized taking c t Sc = 1.The vertical dashed lightgray lines delimitate sections of the basis containing states B † j B † k |v with a fixed value of j + k.

Figure 3
Figure 3: a) Joint density distribution D aa (x, x ) in units of x −2 ω , for two fermions of kind a at positions x and x simultaneously.b) Joint density distribution D ab (x, x ), in units of x −2 ω , for finding a fermion of kind a at position x and one of kind b at position x simultaneously.Both densities were obtained from the numerical solution for λx ω = 30.Details of the calculations are given in Appendices B and C.

Figure 4 :
Figure4: a) Conditional probability P aa (x|0) to find a fermion of kind a at position x when another fermion was already found at the origin.b) Conditional density D ab (x|0) indicating the density of fermions of kind a at position x conditioned on having found a fermion of kind b at the origin.In both plots the solid black curve is the numerical result with λx ω = 30 and n max = λx ω .The dashed red curve is the analytical result for the probability obtained for the point-like hard-core boson limit, and the blue dash-dotted line is the probability predicted by the coboson ansatz in Eq. (11) for N = 2 and λx ω = 30.Details of the calculations are given in Appendices B and C. The vertical light-gray lines indicate the positions ± x ω / √ 2, which are the locations of the maxima of the conditional probability for λ → ∞.

Figure 5 :
Figure 5: a)Off-diagonal correlation parameter g 2 , b) off-diagonal matrix elements and c) diagonal matrix elements of the reduced density matrix ρ ab .Black solid lines correspond to numerical results for λx ω = 30, and red dashed ones for pointlike hard-core bosons.Details are provided in the main text and in Appendix D.

Figure 6 :
Figure 6: Joint density distribution D ab (k, k ), in units of x 2 ω , for finding a fermion of kind a with momentum k and one of kind b with momentum k simultaneously, obtained from the numerical solution for λx ω = 30.Details are provided in the main text and in Appendix E.