Eliminating the wave function singularity for ultracold atoms by similarity transformation

A hyperbolic singularity in the wave-function of $s$-wave interacting atoms is the root problem for any accurate numerical simulation. Here we apply the transcorrelated method, whereby the wave-function singularity is explicitly described by a two-body Jastrow factor, and then folded into the Hamiltonian via a similarity transformation. The resulting non-singular eigenfunctions are approximated by stochastic Fock-space diagonalisation with energy errors scaling with $1/M$ in the number $M$ of single-particle basis functions. The performance of the transcorrelated method is demonstrated on the example of strongly correlated fermions with unitary interactions. The current method provides the most accurate ground state energies so far for three and four fermions in a rectangular box with periodic boundary conditions.

Strongly interacting quantum gases are paradigmatic model systems for the study of strongly correlated manybody physics [1], and can be probed in experiments with exquisite control in the many-particle [2,3] and fewparticle [4,5] regimes. At ultracold temperatures their interactions are accurately described by only a few parameters [6]. Often the s-wave scattering length a s suffices as a sole parameter, or no parameters are needed when a s diverges and the universal regime of unitary interactions is reached [7]. Despite this appparent simplicity, it is nevertheless a great challenge to represent the complicated many-body wave functions in computational approaches [8]. Specifically, a 1/r divergence when two particles with distance r approach each other [9] introduces divergent short-range correlations into the wave function. While exact approaches are limited to no more than three particles [10][11][12][13], computational approaches for larger particle numbers rely on lattice discretisation with renormalised interactions [14][15][16], which is closely related to a renormalised Dirac-delta pseudopotential [17], or the more sophisticated effective Hamiltonian approaches [18,19]. These approaches do not accurately describe the short-range correlations and suffer from slow convergence upon increasing the numerical dimensions of the models (e.g. lattice sites).
In this Letter we apply the transcorrelated method [20] to remove the short-range correlations from the wave function by a similarity transformation of the many-body Hamiltonian. Previously, the transcorrelated method has been applied to Coulomb-interacting electrons [21][22][23][24] and to ultracold atoms in one dimension [25]. In these cases the wave function is non-singular but has a cusp, i.e. is continuous with a discontinuous first derivative [26]. Here, we extend the transcorrelated approach to the hyperbolic singularity and show that it is com-pletely removed. The similarity-transformed, transcorrelated Hamiltonian is free from singular zero-range interactions, which are replaced by non-hermitian two-body as well as three-body terms, and has non-singular eigenfunctions. The advantages of the method are demonstrated by ground-state calculations with stochastic projective diagonalization in Fock-space [27]. For a few fermions with unitary interactions in a box with periodic boundaries we find that the error of the energy is the smallest among the available methodologies. Moreover, this error decays with 1/M , where M is the number of the singleparticle plane wave basis functions. This is the fastest convergence rate so far for these systems.
Zero-range s-wave interactions, as they occur in ultracold dilute atomic gases, are fully described by the Bethe-Peierls boundary condition [9] Ψ(r 1 , r 2 , . . . ) ∼ 1 r ij − 1 a s + O(r ij ) for r ij → 0, (1) where r ij = |r i − r j | is the distance between the interacting particles i and j, and Ψ is the exact many-body wave function. We aim to deal with the divergent short-range part of the wave function with a Jastrow factor e τ by writing Ψ(r 1 , r 2 , . . . ) = e τ (r1,r2,... ) Φ(r 1 , r 2 , . . . ), which defines the transcorrelated wave function Φ. The function τ is chosen as a sum of pair correlation factors In order to satisfy the Bethe-Peierls boundary conditions (1), we choose the correlation factor to obey whereH = e −τ He τ is the transcorrelated Hamiltonian [20]. Here we choose a convenient correlation factor that satisfies Eq. (17) and solve Eq. (5) for the energy E and the non-singular, transcorrelated wave function Φ by diagonalization in a Fock basis. Since Eq. (5) was obtained by similarity transformation from the original Schrödinger equation, it has the same energy eigenvalue spectrum.
It is convenient to define the correlation factor in momentum space, and we choosẽ where k c is a momentum cutoff. The real-space correlation factor is obtained by Fourier transform u(r) = (2π) −2 ∞ 0 dkũ(k)k sin(kr)/r and the corresponding Jastrow factor exp(u) is shown in Fig. 1. More details are provided in the Supplementary Material (SM) [28], where it is shown that u(r) satisfies the asymptotic condition (17). The momentum cutoff k c serves to avoid a singularity in momentum space and damps out the real-space u(r) for large r. The idea behind this construction is that long-range correlations in the transcorrelated wave function Φ can be effectively dealt with by standard approaches, e.g. by expansion in a Fock basis, as we will show, while the Jastrow factor e τ very effectively removes the singular short-range correlations.
For definiteness we consider a system of ultracold atoms of mass m under the influence of an external trapping potential V trap where H 1 = i 2 2m ∇ 2 i + V trap is the single-particle part. The zero-range s-wave interactions between atoms are represented by the Fermi-Huang pseudopotential [29] where g = 4π 2 a s /m is the potential strength. The derivative term regularizes the otherwise pathological contact interaction and enforces the Bethe-Peierls boundary conditions of Eq. (1) [29,30]. This pseudopotential has been applied in exact [10,31] and perturbative [29] treatments, but it has a limitation in the Fock-state based approaches. As the Fock-state basis functions are smooth, the Fermi-Huang pseudopotential reduces to a simple Dirac-delta function, which is pathological in two and three dimensions [17,[32][33][34]. It is suitable for use with the transcorrelated method, however, as long as the Jastrow factor e τ is designed to fulfil Eq. (1).
The transcorrelated Hamiltonian is found by applying the similarity transformationH = e −τ He τ term by term on Eq. (7). Simple functions of the coordinates are unaffected by the transformation because the correlation factor is local in the coordinates, e.g. e −τ V trap e τ = V trap for the trapping potential. The kinetic energy and the Fermi-Huang pseudopotential contain coordinate derivatives and thus generate additional terms according to the Baker-Campbell-Hausdorff expansion [35], which terminates because only first and second derivatives are present. Specifically, e −τ V FH e τ = V FH + [V FH , τ ] and, as we show in the SM [28] for wave functions φ and χ that are bounded and have bounded first derivatives. It can further be shown (see SM) that the matrix elements of the similarity transformed Fermi-Huang pseudopotential χ|e −τ V FH e τ |φ vanish due to cancellation as long as the correlation factor u(r) is chosen to have the appropriate short-range asymptotics of Eq. (17). Thus, the singular pseudopotential is removed and the transcorrelated Schrödinger equation (5) can be solved with a non-singular wave function Φ. This insight presents the main result of this Letter.
The transcorrelated Hamiltonian still acquires terms that originate from the kinetic energy operator, and fi-nally reads The new terms that appear in the square brackets represent an effective interaction potential that is less singular than the Fermi-Huang pseudopotential. The leading singular term is (∇ i τ ) ∇ i ∼ j r ij /r 2 ij ∇ i , which has a 1/r divergence and is also non-hermitian. Similar to the Coulomb potential it leads to a cusp feature, where the transcorrelated wave function Φ is continuous with a discontinuous first derivate, as is shown in the supplementary material for two particles [28]. As a consequence, the Fourier-transformed wave function Φ in momentumspace decays with 1/k 4 instead of the 1/k 2 decay of the original wave function Ψ for large k. It is this rapid decay for large k that makes it feasible to expand the problem in a plane-wave basis without the need for a renormalized (or running) coupling constant.
Although the similarity transformation eliminates the singularity from the wave function, it introduces new challenges for numerical calculations. The non-hermitian term (∇ i τ ) ∇ i prevents in general the variational minimalization of the energy. However, the energy still can be well approximated by diagonalization in a finite basis set, or by stochastic projection to the ground state as we will demonstrate. As a consequence of the nonhermiticity, the approximate energies no longer provide an upper bound to the exact ground state energy.
The terms 1 2 ∇ 2 i τ and 1 2 (∇ i τ ) 2 scale like a long-range interaction with 1/r 2 ij and partly compensate each other but leave an uncompensated three-body interactiona harbinger of the Efimov effect [36,37], which permits three-body bound states for resonantly interacting bosons. Three-body interactions are common in the transcorrelated method and the technical challenges have already previously been handled in the Hubbard model [38] and for electronic structure calculations in atoms [24]. The term (∇ i τ ) 2 also produces an infinite summation when evaluated in a plane-wave basis, which, however, can be performed efficiently as described in in the SM [28].
We have performed numerical calculations with the transcorrelated method on the example of two to four fermions with unitary interactions (i.e. 1/a s = 0) in a cube with periodic boundary conditions. The transcorrelated HamiltonianH is expanded as a finite matrix in a Fock basis of antisymmetrized products of single-particle plane waves with a momentum cutoff. For two particles the ground state energy is calculated with (numerically) exact diagonalization. For three and four fermions the full matrix diagonalization was not possible due to the enormous size of Hilbert space. Hence we used a stochastic projection method known as Full Configuration Interaction Quantum Monte Carlo (FCIQMC) [27,39] to  Figure 2. The ground-state energy of two particles with unitary interactions (1/as = 0) vs. the inverse number of single-particle basis functions 1/M from the transcorrelated method (circles) and with renormalized Dirac delta (crosses) as discussed in the text. The reference energy Eex = −3.786005 2 /mL 2 (green solid line) is obtained from Ref. [40]. The inset shows the difference to the reference energy on a log-log plot. The obtained linear relation indicates a power-law scaling for the energy error ∝ 1/M 1/3 for the renormalization approach and ∝ 1/M for the transcorrelated approach.
obtain the ground state energies. A detailed description of our numerical methodology is provided in the SM [28].
Numerical results for the ground state energy for two particles (either distinguishable fermions or bosons) are presented in Fig. 2. In this and the following figures, the energy is shown as a function of the inverse of the basis set size M , where 3 √ M is the number of single-particle planewave basis functions per linear dimension of the cube. Hence, the zero on the x-axis represents the complete basis set limit. The transcorrelated energies are compared to standard lattice renormalization [14], where a running coupling constant g 0 is scaled with the number of lattice points M as g −1 0 = m/4π 2 a s − mKM 1/3 /4π 2 L (K = 2.442749607806335...) [41]. It is not only seen that the transcorrelated method gives smaller errors by orders of magnitudes for the same M , but also that scaling of the errors with M follows a faster power law decay. For the renormalization approach, we find a scaling of M −1/3 , which is consistent with the previous results from lattice calculations [14,15,42]. In the transcorrelated approach the error decays with M −1 . This is the same scaling as obtained for Coulomb-interacting systems, e.g. the homogeneous electron gas, which is consistent with the Coulomb-like nature of the transcorrelated Hamiltonian. The M −1 scaling is the fastest scaling we found in the literature, and is shared, e.g. with the improved lattice action used for Auxiliary Field Quantum Monte Carlo (AFQMC) calculations by Endres et al. [16] or the renormalized lattice Hamiltonian with "magic" dispersion relation discussed in Ref. [14]. Results from the transcorrelated method ("TrCorr") are compared with semi-analytical results from Ref. [43] ("Scattering Theory") and AFQMC from Endres et al. [16] ("Endres AFQMC"). The yellow band marks the standard error of the extrapolated AFQMC results. For comparison we show renormalised lattice calculations with FCIQMC with different single-particle dispersions: "Hubbard", "Quadratic" and "Quartic" as in Ref. [15] and the "Magic" dispersion from Eqs. (122, 124) in Ref. [14].
A comparison of the transcorrelated method with literature results for the three fermion ground state energy is shown in Fig. 3. While the M −1 scaling (linear on the scale of Fig. 3) can be observed for the "Endres AFQMC" and "Magic FCIQMC" values, it is remarkable that the transcorrelated results are much more accurate already for very modest basis set size and hardly distinguishable from the reference values on the scale of the figure. Moreover, we find that approximate calculations avoiding the numerically expensive three-body excitations achieve the same accuracy within our statistical errors. Details are provided in the SM [28]. Figure 3 also shows renormalized lattice calculations with different single-particle dispersion relations as discussed in Ref. [15] obtained with FCIQMC. Since they are expected to show slower scaling than M −1 , the energy dependence does not appear linear in Fig. 3.
The transcorrelated energies for three particles are shown again in Fig. 4 with a magnified energy scale and with different momentum cutoffs k c in the correlation factor of Eq. (6). It is seen that the asymptotic regime of M −1 scaling of the energy error is only reached for the larger basis set sizes. With the known asymptotic scaling properties we can determine the ground state energies in the infinite basis set limit by extrapolation. The results for three fermions are illustrated in Fig. 4, where it seen that the extrapolations with two different k c values are consistent with each other as well as with the literature results from scattering theory The reference value from Ref. [43] ("Scattering Theory") is shown with error band in green. The error band of the extrapolated result from "Endres AFQMC" [16] does not fit on the scale of the plot but can be seen in Fig. 3. The purple and the red bands indicate the standard error band obtained from linear fitting of the transcorrelated FCIQMC ("TrCorr") results (four largest M values fitted) [28]. EF = 2π 2 /mL 2 is the Fermi energy.
[43] and AFQMC [16], while they have much smaller error bars than previous results. As the final value for the ground state energy for three fermions we obtain E/E F = 0.373444±0.000028. Compared to the results of Endres et al. [16] of E/E F = 0.3735(+0.0014/ − 0.0007) the error is reduced by more than an order of magnitude. The results from transcorrelated and renormalized calculations for a four-fermion system are shown in Fig. 7 where they are also compared to literature results with AFQMC (for details see the SM [28]).
While we have presented numerical results for the homogeneous Fermi gas, the approach can be easily extended to include trapping potentials or external gauge fields. The transcorrelated approach is thus well suited for highly precise calculations on correlated few-atom systems in microtraps [4,5]. Moreover, we hope that the transcorrelated formalism brings new insight into the treatment of the singularity in the wave function and that it provides a useful theoretical tool in other perturbative and exact computational approaches.  Figure 5. The ground-state energy of four fermions extrapolated to the infinite basis set limit from different methods. AFQMC results from Bour et al. [42] with Hamiltonian lattice 1 ("Bour 1") and 2 ("Bour 2") and Euclidian lattice ("Bour 3") and from Endres et al. [16] with O(4) ("Endres 1") and O(5) ("Endres 2") scaling are compared with FCIQMC calculations using the transcorrelated method with kc = 2π/L ("TrCorr") and renormalized lattice calculations following Ref. [15] using the Hubbard ("Hubbard") and quadratic ("Quadratic") dispersion relations. Details of the extrapolation procedure and numerical values are presented in the SM [28].
Apārangi. We also acknowledge support by the NeSI high-performance computing facilities.

THE CORRELATION FACTOR IN REAL SPACE
In this section, we examine the real-space form of the correlation factor, The Fourier-transform of function (11) can be calculated analytically, The boundary condition can be reproduced by expanding function v(r) and w(r) in Taylor series around r = 0, where γ is the Euler-Mascheroni constant. Calculating e u(r) , we obtain back the hyperbolic singularity for the Jastrow-factor, e u(r) = e 1−γ+ 4 askc π 1 k c r − 1 k c a s + O(k c r) .

MATRIX ELEMENT OF THE TRANSCORRELATED FERMI-HUANG PSEUDOPOTENTIAL
We consider the matrix element of the transcorrelated Fermi-Huang pseudopotential and show that it vanishes, if evaluated with wave functions that are bounded and have a bounded first derivative almost everywhere.
In order to show that, let us consider the transcorrelated Fermi-Huang pseudopotential, where the commutator can be evaluated if we apply the substitution V FH = g i<j δ (r ij ) ∂ ∂rij r ij , The partial derivative with respect to the separation r ij = |r i − r j | is defined in the usual way, where both particles i and j move while the center-of-mass 1 2 (r i + r j ) is held constant, as are the orientation of the vector r ij = r i − r j , and all other particle coordinate vectors r k for k = i, k = j. Since τ depends on the separations of all particle pairs, the chain rule will generate many term, most of which, however, vanish.
In order to evaluate the derivative ∂τ ∂rij , let us substitute in the expansion of τ in pair correlation functions, into ∂τ ∂rij , The last term on the right-hand side is zero as r kl does not depend on r ij . For the second and the third terms on the right-hand side we can apply the chain rule, ∂r kj ∂r j,p ∂r j,p ∂r ij where index p goes through the three spatial directions and θ jil is the angle between r ji and r il . Using the short-range behavior of the correlation factor, the first derivative of u(r) can be evaluated for short interparticle separations, Substituting Eqs. (15)-(18) into Eq. (14), the explicit expression can be obtained for ∂τ /∂r ij , Using the expression above, the matrix element of the commutator expression in Eq. (12) is can be expressed for short interparticle separations, Assuming that the functions χ and φ are bounded, the last two terms in Eq. (20) are zero as δ (r ij ) r ij gives zero after performing the integral either for r i or r j . Although the second and the third term on the right-hand-side of Eq. (20) can be finite at the coalescence points r i = r l and r k = r j , they is zero everywhere else. Since the coalescence points form a ste of measure zero, these terms yield zero after integrating over the remaining variables. This leads to a matrix element of the Dirac delta function: Due to the bounded nature of the functions φ and χ, the matrix element of the Fermi-Huang pseudopotential also reduces to the matrix element of the Dirac-delta function, but with the opposite sign, where we have assumed that χ and ∂φ/∂r ij are bounded. Equation (22) shows that a matrix representation of the the (physically meaningful) Fermi-Huang pseudopotential with sufficiently smooth (and bounded) basis functions is equivalent to the bare Dirac-delta pseudopotential, which is pathological in the sense that the infinite basis set limit does not exist. After the transcorrelated similarity transformation, however, we obtain the two matrix elements (21) and (22), which cancel each other and thus eliminate the irregular behavior in the matrix representation.

SMOOTHNESS OF THE TRANSCORRELATED EIGENFUNCTION FOR TWO PARTICLES
In this section we investigate the transcorrelated eigenfunction for two bosons or two fermions with different spin flavors. We show that the singularity is reduced in the transcorrelated Hamiltonian due to the similarity transformation. Consequently, the transcorrelated eigenfunctions are not singular, there is only a cusp at the particle-particle coalescence point.
We consider the two-particle Hamiltonian without trapping potential (V trap = 0), Separating the center-of-mass from the relative motion coordinates, we obtain, where r = r ↑ − r ↓ , µ = m/2 and the center-of-mass is described by free-particle motion. Let us apply a similarity transformation on the Hamiltoninan (23), Using Eqs. (13) and (17) τ can be given explicitly at small interparticle separation, with which the derivatives of τ in Eq. (24) can be expressed, Substituting back Eqs. (26)- (28) into (24), an explicit from for the Hamiltionian can be obtained for short distances, In order to obtain the transcorrelated eigenfunction, let us substitute the Hamiltonian into the Schrödinger equation, where E = E + O r 0 . Due to the spherical symmetry, we can transform the differential equation into polar coordinates, and we can consider only the s-wave solutions, The differential equation can be solved for small interparticle separation, where c 1 , c 2 and E can be determined only if we know the solution in the whole space.
Differentiating the wave function we notice that it has a linear term, where the prefactor b, before the linear term is Considering the spherical symmetry, we obtain a function which goes linearly to c 1 +c 2 around the origin and forms a cusp. This function is not-singular and continuous, however, its first derivative is discontinuous. Therefore, the transcorrelated transformation reduces the singularity both in the Hamiltonian and the eigenfunction.

SECOND QUANTIZED FORM OF THE TRANSCORRELATED HAMILTONIAN
In this section, we give an explicit expression for the second quantized form of the transcorrelated Hamiltonian in a rectangular box with periodic boundary conditions.
As we have been discussed previously in the main body of the paper the transcorrelated Hamiltonian only depends on the transcorrelated kinetic energy operator, where |χ K and |χ L are determinants for fermions and permanents for bosons. Following the description in Ref. [23] the secondquantized form can be easily determined, where a † k,σ creates a one-particle plane wave state with momentum k and spin σ, L is the length of the unit cell, and Θ σσ = δ σσ for bosons and Θ σσ = 1 − δ σσ for fermions. The tensors T and Q can be expressed explicitly,

NUMERICAL EVALUATION OF THE INFINITE SUMMATION IN EQ. (31)
In this section we describe the algorithm, which we used to evaluate the infinite summation in Eq. (31).
First, let us realize that we can restrict the indices in the summation from below due to momentum cutoff in the correlation factor (11), As we can see in Eq. (32), the summation goes to infinity, which prohibits the exact evaluation. However, an accurate approximate value can be obtained if we partition the summation in Eq. (32) to a summation inside a sphere with radius k int and a summation outside this sphere, As w(k, k ) decays with k −6 at larger values of k, we can approximate the summation with an integral in R(k, k int ), Due to the additional conditions in the sum (32), further restrictions apply at the boundaries of the integral, when In order to avoid the complicated limits of the integration, we choose k int large enough in a way to never satisfy condition (38). After some algebra, it can be shown that it is enough to consider condition which is usually fulfilled as k and k c is kept small to limit the size of the Hilbert-space and to enhance the effect of the correlation factor. In order to evaluate the integral (37), we consider the Taylor-expanded form ofũ(k) in Eq. (11), In this paper we consider unitary interactions, where the integral (37) can be exactly evaluated, In this paper we applied Eqs. (33)- (36) and (40) to evaluate the matrix element W (k) for up to k = 16π/L with k int = 1600π/L. Considering the maximal value of k to 16π/L we found uncertainty in the values of W (k) only in the seventh and in the eighth significant digits. As the energy scales linearly with the error in the matrix elements, the error should appear in the energy in the same order. Moreover, the accuracy of the integral approximation is also checked numerically by comparing the energies from k int = 1200π/L and k int = 1600π/L calculations. We did not find any significant difference in examples of two, three and four fermions.
For two particles the convergence of the energy is demonstrated in Fig. 6. This uncertainty seems adequate for our numerical calculations, where the uncertainty of our final results was in the fourth and fifth significant digits.

DETAILS OF THE NUMERICAL CALCULATIONS
For the numerical calculation we used the NECI code [39], where transcorrelated Hamiltonians including threebody excitations had previously been implemented for the homogeneous electron gas [23], the Fermi-Hubbard model [38], atoms, molecules [24], and the Fermi gas in one dimension [25]. In the context of this project we have further extended the capabilities of the NECI code by including the transcorrelated Hamiltonian for the unitary Fermi gas in three dimensions.
For two-particles, non-hermitian exact (deterministic) diagonalization is applied in NECI using an external Lapack library [45]. For three and four fermions the Hilbert-space is too large for deterministic diagonalization. Hence the Full Configuration Interaction Quantum Monte Carlo (FCIQMC) algorithm [27,39] is applied to obtain the ground-state energy.
One of the elementary parameters of the FCIQMC algorithm is the number of the walkers [27,39]. It controls the resolution of the wave function and the memory usage of the algorithm. In this algorithm a minimal number of walkers is required to eliminate the sign-problem. The minimal number is determined by the annihilation plateau [46], which appears in the number of walkers during the imaginary time evolution.
For all calculations for three fermions and for the lattice-renormalized calculations for four fermions, we were able to apply a large enough walker number to get above the annihilation plateau. However, for the transcorrelated calculations with four fermions the annihilation plateau was too high for the available numerical resources. In these cases we applied the initiator approximation [47], which has proved to be efficient for electronic structure calculations [51][52][53]. While this approximation causes a systematic bias in the calculations, the bias disappears when increasing the number of walkers. For all results shown, the number of walkers was increased until the changes in energy were insignificant compared to the statistical error bars. Another systematic bias, the population control bias, is known to affect FCIQMC calculations with small walker number but is well below the statistical error for the parameters considered in our calculations. We thus expect the FCIQMC results presented in this work to be essentially free of any systematic bias. Complete basis limit and uncertainty The complete basis limit of the energy, E cb , can be easily determined by considering the observed asymptotic scaling of the error in the energy. According to our numerical calculations for two particles in Fig. 2 of the main text, it should be inversely proportional to the number of the single-particle plane waves, with a fitting parameter β, where E F is the Fermi energy, and M is the number of single-particle basis functions. E cb is obtained from linear extrapolation using the POLYFIT code [44]. The standard error band of the linearly fitted curve is described in detail in Ref. [48]. Considering the linear parametrization (41), the standard-error SE(1/M ) for every value of 1/M can be calculated from where n is the number of the data points, M −1 i is the value of M −1 in the i-th data point, avg(M −1 ) is an average of M −1 over all the data points, and σ 2 is the residual sum of squares (RSS) divided by n. σ 2 can be calculated either from the standard deviation of M −1 (σ M −1 ) or the standard deviation of E cb /E F ( σ E cb /E F ),  Table I. Numerical values of the data shown in Fig. 5. Ground-state energies for two spin-up and two spin-down particles. The renormalization methodology for the different dispersion relations is described in the main text and in detail in Refs. [14,15].
The standard error bands from Eqs.  Table I.

Approximation for the three-body interaction
Evaluating the transcorrelated Hamiltonian during any diagonalization procedure requires increased numerical effort compared to the renormalized lattice Hamiltonian. The largest part of the increased effort can be attributed to the three-body term and thus scales with N 3 , where N is the number of particles. The results presented in the main part of the paper were computed by fully including all three-body excitations. In this section we discuss an approximate procedure that only requires evaluating effective two-body matrix elements, and reduces the numerical effort significantly, while still producing highly accurate results.
Ground state energies computed with this approxi-mate Hamiltonian are compared to the full transcorrelated Hamiltonian are shown in Fig. 8 for three fermions and in Fig. 9 for four fermions. We find that the approximate results agree within the Monte Carlo (statistical) error bars with the full transcorrelated results.