{'fraction_non_alphanumeric': 0.031441745634099405, 'fraction_numerical': 0.011182514613513415, 'mean_word_length': 4.707418151678409, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Simulating strongly correlated multiparticle systems in a truncated Hilbert space

Representing a strongly interacting multi-particle wave function in a finite product basis leads to errors. Simple rescaling of the contact interaction can preserve the low-lying energy spectrum and long-wavelength structure of wave functions in one-dimensional systems and thus correct for the basis set truncation error. The analytic form of the rescaling is found for a two-particle system where the rescaling is exact. Detailed comparison between finite Hilbert space calculations and exact results for up to 5 particles show that rescaling can significantly improve the accuracy of numerical calculations in various external potentials. In addition to ground state energies, the low-lying excitation spectrum, density profile and correlation functions are studied. The results give a promising outlook for numerical simulations of trapped ultracold atoms.


I. INTRODUCTION
Numerically solving the many-particle Schrödinger equation in one spatial dimension with short-range interaction requires a truncated Hilbert space to make the calculations computationally tractable. This introduces inaccuracies especially in strongly correlated regimes, systems with many degeneracies, and near avoided crossings. Here, we demonstrate a simple method that significantly improves results by rescaling the interaction strength.
Understanding the complexities of many-body quantum physics remains a grand challenge, in particular with regard to excited states and dynamical problems. Ultracold atoms have recently become a testbed for both analytic and computational approaches to understanding quantum many-body physics. These systems can be prepared with few to millions of particles and the experimental control over interactions, system size, and quantum states is rapidly improving [1][2][3].
However, exact theoretical results are rare, especially for few-particle systems. Approximate approaches, like the mean-field Gross-Pitaevskii equation, accurately model a large number of weakly interacting atoms. Now, there is a considerable amount of research that focuses on the creation of entanglement in ultracold atomic systems. These require calculations at a single-particle level. Much success with fixed basis-set expansions has been achieved using Wannier functions of an optical lattice leading to the Bose-Hubbard model [4] and with harmonic oscillator eigenfunctions [5] or a plane-wave basis [6,7] for trapped systems. These approaches are limited by accessing only a finite Hilbert space. Methods like time-evolving block decimation [8] and time-dependent density matrix renormalization [9,10] exploit spatial entanglement properties to allow for time-dependent simulations of larger multiparticle problems.
A different approach is taken by the multiconfigurational time-dependent Hartree (MCTDH) method [11][12][13]. It variationally optimizes a set of single-particle functions, that define the accessible Hilbert space at any given moment in time. It has been successfully applied to ultracold gases [14][15][16]. Even so, the truncated basis leads to errors in strongly interacting systems. Here, we analyze the origin of errors that arise from basis-set truncation, suggest a scheme to compensate for them, and benchmark the performance of the scheme.
Trapped bosonic atoms confined to a single spatial dimension undergo a crossover from a Bose-Einstein condensate to the fermionized and strongly correlated Tonks-Girardeau (TG) gas [17] as repulsive short-range interactions are tuned from weak to very strong. This system is well suited to benchmark computational many-body methods as exact results are available from integrable models in the limits of infinite interactions (TG [17]) or vanishing external potential (Lieb-Liniger model [18]). Previous calculations of the ground state with exact diagonalization in a fixed product basis [5,19,20] and with the MCTDH method [14,21] have shown a peculiar feature: Specific finite values of the interaction strength in these calculations provide a reasonable approximation of the completely fermionized state, which in exact theory is only reached for infinite interaction. Our calculations show that larger values of the interaction produce worse results in the approximate calculation for the energy and lead to spuriously enhanced density oscillations. In this paper we provide a clear explanation of this effect based on exact results for two particles.
We find that the origin of the problem lies in the truncation of the Hilbert space. Previously, elaborate schemes have been developed to generate effective Hamiltonians that remove this problem, however these techniques are cumbersome and require careful application [22][23][24][25]. A simpler, yet powerful method is desirable to accurately describe different interacting regimes for a variety of applications [26][27][28][29]. Here we present a rescaling method that increases accuracy or decreases the required Hilbert space by mapping the interaction in the simulations to the physical interaction. Therefore numerical simulations can be sped up and used to approach previously unaccessible regimes and systems.
In this paper we initially present an exact rescaling procedure for two particles in a one-dimensional periodic box as described previously in Ref. [7]. A more general scheme is proposed that only requires the knowledge of the numerical value of the interaction strength at the TG energy. This approach is closely related to the known renormalization of the two-body T matrix but extends this method into the strongly interacting regime. The approach is tested for different system geometries, particle numbers and external potential shapes. The results are compared to analytical results for two atoms in a harmonic potential [30], the Lieb-Liniger description of atoms confined to a one-dimensional periodic system [18] and infinitely strong interacting atoms in one dimension described by a TG gas [17].

A. Finite basis-set expansion
Consider N interacting quantum particles in one dimension subject to an external potential V ext (x, t), which may depend on time: where is the one-body Hamiltonian and for simplicity we assume that all particles have the same mass m and interact by contact interactions of strength g. This model describes, e.g., quantum gases of ultracold atoms in an elongated and tightly confining trapping potential [31]. Although we specifically deal with identical bosons later on, at this point considerations are not restricted to a specific quantum statistics and the Hamiltonian (1) could describe single-or multicomponent Bose or Fermi gases or mixtures thereof. All these possibilities are of interest and related to actual or possible experimental scenarios. The true N -particle wave function can be expanded in the form (2) in a basis consisting of products of single-particle wave functions where we have introduced the multi-index J = (j 1 , . . . , j N ) for the set of N single-particle indices. The single-particle wave functions (or mode functions) are mutually orthonormal by φ * j (x)φ k (x)dx = δ jk . In practice we chose a finite set of M single-particle functions to define the finite basis-set expansion. The size of the finite basis of multiparticle Hilbert space is nominally M N , although for identical particles this number can be significantly reduced by accounting for bosonic or fermionic exchange symmetry in a well-known manner by symmetrizing the products to permanents or Slater determinants, respectively.
The standard procedure for approximating the full problem (1) is to simply truncate the Hilbert space to span a finite basis and represent the Hamiltonian as a finite matrix with elements We denote the solution of the corresponding matrix Schrödinger equation for the coefficient vector withĀ J and the eigenvalues of the truncated matrixĒ ν . These approximate the exact coefficient vector A J and eigenvalues E ν , which form part of the spectrum of H mb . Since the basis-set truncation can be understood as a variational procedure with a restricted variational space, the approximate ground-state energyĒ G ≥ E G is an upper bound for the true one. By increasing the number of modes M and thereby computational space, the approximation improves and the approximate energy converges to the correct one. However, for an interacting system this convergence is painfully slow, as seen in the inset in Fig. 1. The purpose of this paper is to improve the approximations for a given, finite, basis-set expansion. While previous works have explored elaborate procedures to approximate an effective finite Hamiltonian with the exact eigenvalues E ν [22][23][24], we show that it is possible to obtain significant improvements by simply rescaling the interaction parameter g.

B. Two particles on a ring
For two particles without external potential an analytic expression for the rescaled interaction constantg can be found that reproduces the correct ground-state energy and the Fourier components of the wave function up to a cutoff. The rescaling procedure is thus exact for this case.
Consider N = 2 bosons in a one-dimensional box with length L and periodic boundary conditions. This model might be realized with ultracold atoms in a tightly focused ring trap. The Hamiltonian (1) simplifies to where x i is the spatial coordinate of atom i and g the physical interaction strength. A finite basis-set expansion for this system is where we assume the number of momentum modes M in the expansion to be odd, for simplicity. First consider the case for the full Hilbert space with M = ∞. To find the ground state we substitute Eq. (6) and Eq. (5) into the Schrödinger equation (E − H 2p )ψ = 0. Projecting to zero-momentum solutions by multiplying with exp[−i2πk(x 1 − x 2 )/L] and integrating over the particle coordinates leads to Here, E 0 = 2π 2 2 /(mL 2 ) is the smallest nonzero singleparticle energy. The right-hand side of Eq. (7) is independent of k and therefore C −q,q ∝ (E − 2E 0 q 2 ) −1 and which relates the exact energy E to the physical interaction strength g. Employing the finite basis-set expansion limits the sum over q to only M terms and thus yields a different, approximate, valueĒ for the energy. Here, we choose a different path and introduce the rescaled interaction strengthg by where the last equation defines the constant g 0 . While Eq. (9) guarantees thatg yields the exact energy E it is Eq. (10) that provides the rescaled interaction constant g as a function of g:g Therefore, solving the Schrödinger equation with this rescaled valueg instead of g in the Hamiltonian (5) in the finite basis-set expansion gives the exact energy E for the physical interaction strength g. Also the finite number of expansion coefficients C k,−k is identical to the full Hilbert space expansion and, according to Eq. (7), is given by where The value of g 0 is found by expanding the sums in powers of 1/M from For a large number of modes and sufficiently low energy E [with the approximate condition M 2 ≫ E/(2E 0 )] we thus find to leading order that g 0 is independent of the energy and given by Figure 2 shows how the rescaling in approximation (15) significantly improves the ground-state energy. For M = 6 the energy after rescaling is only 3% off the exact energy whereas before rescaling the difference was ∼ 20%. The constant g 0 is readily interpreted from Eq. (11) as the value of the rescaled interactiong where the physical interaction g reaches infinity. Physically, this is the TG limit where the contact interactions are arbitrarily strong and the system is strongly correlated. Our results thus indicate that a finite basis-set expansion with a finite value of the rescaled interaction constantg = g 0 can produce exact results for infinite values of the physical interaction strength. Since the energies and wave functions in the TG limit can also be calculated by the Bose-Fermi mapping [17], this limit is an important reference point.

C. Empirical rescaling
The analytical results of the preceding section motivate us to extend the idea of rescaling to more than two particles and to problems with arbitrary external potential, which, generally, cannot be solved analytically. We propose to replace the Hamiltonian (1) by the rescaled versionH where we have only changed the value of the interaction strength from the physical value g. Within a finite basis-set expansion we obtain a matrix with elements H IJ , which is to be used for numerical simulation. The smallest eigenvalue is the approximate ground-state en-ergyĒ G (g). For the rescaled interactiong we continue to use Eq. (11) and determine the value of g 0 by requiring that the approximate ground-state energy at g 0 equals the exact energy of the TG limit: A finite solution for g 0 can always be found since the approximate value ofĒ G (g) overestimates the real value andĒ G (g) is a monotonously growing function ofg with E G (0) ≤ E TG . The value of E TG is equal to the groundstate energy of a system of noninteracting fermions due to the Bose-Fermi mapping theorem [17]. It is given by the sum of single-particle energies, which can be determined analytically or numerically to very high precision.
This rescaling procedure is very simple to implement, since only knowledge of the value of the single constant g 0 is needed. We assume here for simplicity that g 0 is independent of the energy of the wave function and can be used for any set of (low-lying) excited states as well as time-dependent processes involving such states. As an approximation, this is consistent with the analytic findings for two particles from the previous section.

D. Connection to T matrix renormalization
The rescaling approach introduced in Sec. II C is closely connected to the well-known renormalization of the scattering T matrix (for a readable account see Ref. [32]). In this approach, the coupling constant of an effective many-body theory is renormalized to produce the correct T matrix for the two-body scattering problem, which is known exactly. Indeed, Eq. (11) for the rescaled interaction strength together with Eq. (15) for the value of g 0 are identical to the result of Ref. [32] for the renormalized T matrix of a high-momentum cutoff introduced there due to discretization of space (identifying the computational volume of momentum space 2πM/L with the Brillouin zone of Ref. [32]). This result was rigorously derived for weak interactions and high cutoffs. In Sec. II B we have extended the approach to the strongly interacting regime for two particles where now the rescaled interaction constant in principle becomes energy dependent. The extension to nonperturbative, strongly interacting multiparticle systems in Sec. II C with Eq. (8) for the value of g 0 is heuristic in nature. The purpose of the numerical studies reported in the following section is to establish the usefulness and to quantify the accuracy of this approach.
Rescaling improves numerical calculations by either increased accuracy or decreased computational effort for a given desired accuracy (CPU time ∼ M 3 for our five particle simulations). Forg > g 0 we find unphysical results. Therefore, CPU times are further reduced as only results forg ≤ g 0 are used in the rescaling.
Before detailing our numerical results and benchmarking, we introduce the MCTDH method suitable for general trapping potentials.

E. Variationally optimized basis
The finite basis-set expansion (2) and (3) expands the multiparticle wave function in products of a given orthonormal set of single-particle functions. Optimizing this set variationally is the aim of the MCTDH method. This is achieved with a time-dependent variational ansatz with a fixed number of basis functions without further approximation [11], in contrast to the density matrix renormalization group approach, which further decomposes the problem into bipartite systems [9,10] before finding an optimized basis.
The single-particle wave functions φ j (x, t) are now allowed to be time dependent in addition to the time dependence of the coefficients A J (t). The variational equations of motion for the coefficients and the single-particle wave functions are discretized in time (and space) and solved numerically with the MCTDH program package [12,13].
Eigenenergies and eigenstates can be obtained by the relaxation method [33], i. e., by propagation in negative imaginary time. Here we use a modification, improved relaxation [34], which can also be used in block form [35] to simultaneously compute a set of eigenvectors. For a fixed basis of single-particle functions, the MCTDH method reduces to the exact diagonalization method. Therefore, the key advantage is that the MCTDH method always stays in a variationally optimized basis. This generally ensures higher accuracy or keeps the Hilbert space as small as possible for a given accuracy. The MCTDH method is therefore effective at describing different trapping potentials. However, it requires a numerically unattainable number of single-particle functions to describe strongly interacting systems.
In practice, the single-particle wave functions are represented on a large number of primitive basis functions, the choice of which depends on the spatial geometry (typically 128 and 162 grid points in the harmonic oscillator discrete variable representation for the harmonic and the double-well potentials respectively and 155 points in the exponential discrete variable representation as per Ref. [12]). The setup of the MCTDH method for ultracold atom calculations follows Ref. [14]. In particular, the contact potential is replaced by a narrow Gaussian for practical purposes. We have carefully studied the dependence of our calculations on the parameters of the primitive basis and the Gaussian interaction to make sure that the results presented here are converged and the influence of these further approximations is negligible.

III. BOSONS IN A RING
We first benchmark the rescaling method of Sec. II.C. for the Lieb-Liniger model [18] of bosons in a box with periodic boundary conditions, where exact eigenstate energies are easily available for any values of the interaction strength and particle number. Figure 1 shows the relative deviation from the exact ground-state energy obtained from MCTDH calculations before and after rescaling. Rescaling is seen to significantly reduce the error. The rescaled results reproduce the exact energy in the TG limit by construction but the improved accuracy for finite interaction strength is nontrivial. As expected, both the raw and rescaled data improve with increasing M . The maximal errors for, e.g., M = 13 can be reduced from ∼ 23% down to ∼ 0.6%, i. e. a factor of 40. This is a significant improvement.
In order to test how well the wave function is approximated in the TG limit we investigate the reduced pair FIG. 3. Density-density correlation function for five particles and M = 13 single-particle functions in a periodic box. This plot compares the reduced two-particle density ρ (g) 2 of Eq. (18) for different interaction strengths from MCTDH calculations without rescaling with the exact TG expression (the area below this curve is shaded). We find a good agreement with the exact solution close to g0 = 5.3958E0L. This is shown in the inset where the overlap F (g) (ρ (19) is plotted as a function of the interaction g. There, the maximum occurs close to g0. density ρ (g) 2 (0, x) = · · · dx 3 . . . dx N |Ψ (g) (0, x, x 3 , . . . , x N , t)| 2 .
(18) The results of MCTDH calculations are shown in Fig. 3 for different values of the interaction g along with the exact solution in the TG limit ρ T G 2 (0, x) = [(2 + cos x) sin 2 x 2 ]/3π 2 [17]. Agreement between the finite basis-set expansion at the value of g 0 = 5.3958E 0 L and the exact solution is reasonable. The discrepancies can well be explained by the absence of short-wavelength modes in the finite expansion. Most importantly, the agreement between the approximate and exact functions clearly worsens for interaction strengths larger than g 0 , where the structure in the finite basis-set calculation overshoots.
In order to quantify the quality of the approximation, we define the overlap with the exact solution at the TG regime by which is shown in the inset of Fig. 3. We find that the maximum of F (g) is reached near g 0 , which justifies the choice of g 0 through Eq. (17). In order to study the scaling with particle number, we have performed exact diagonalization calculations for N = 2 . . . 6 particles in a ring of length L = N/ñ, where the particle densityñ was fixed to the same value as in Fig. 3. Figure 4 shows the maximal relative deviation of the rescaled energy as a function of particle number. The results show that retaining M = 2N − 1 singleparticle functions gives a maximal relative error bounded by |E/E LL − 1| 1.5% for any value of g > 0. Although for clarity only odd numbers of single-particle modes M are shown in Fig 4, we have verified that these findings also hold for even M .
If we assume that the scaling M = 2N − 1 is sufficient to maintain accuracy uniformly over the whole range of interaction strengths also for larger particle numbers, we still find that the numerical effort of MCTDH or exact diagonalization studies grows exponentially with N . This limits studies of this sort to fewer than 10 particles with current-day computers.

IV. BOSONS IN A HARMONIC POTENTIAL
In order to test the rescaling method in the presence of an external potential we consider a harmonic trapping potential For this and the following section we use as natural energy and length scales the harmonic oscillator level spac- ing E 0,H = ω and L H = mω , respectively, where ω is the harmonic oscillator frequency. Figure 5 compares the single-particle density of the TG wave function with MCTDH simulations for five particles at different values of the interaction strength without rescaling. The agreement between the exact and the finite basis-set results at the value of the interaction strength g 0 = 15.516 . . . E 0,H L H is much better than for larger values of g. This result is analogous to the ring case and supports our proposition that g 0 from Eq. (17) provides a good approximation for the wave function.
Full analytical solutions are available for two particles in a harmonic potential. The exact energies for the ground and excited states are found from the implicit equation [30] where Γ is the gamma function. Figure 6 compares the ground-state and a few excited-state energies with rescaled MCTDH results. The highest few states and the ground state have been chosen as they are more likely to show the largest errors within the block-diagonalization scheme (see below). We emphasize that a single value of g 0 = 29.75 obtained from Eq. (17) was used for rescaling. The procedure not only maintains the correct order of excited states but also quantitatively describes the excited-state energies remarkably well.

V. BOSONS IN A DOUBLE WELL
The rescaling method is further tested for the description of excited states of the TG gas in a double-well po- tential. Exact solutions are not available for interacting multiparticle systems in the double well except for the TG gas, where their exact energies are simply found from sums of single-particle energies (which can be determined with great accuracy) due to the Bose-Fermi mapping of Ref. [17].
We consider the double-well potential with varying barrier height h at the center. Figure 7 shows the energies for the ground and excited states from rescaled MCTDH calculations and analytical results. All MCTDH results were obtained in the TG regime atg = g 0 . The value of g 0 has been determined for each value of h from a dedicated ground-state calculation. The results reported in Fig. 7 were calculated using the block-improved-relaxation method [35] to efficiently determine several states simultaneously. Since an identical set of single-particle functions is used for the finite basis-set expansion for all states, the result is not as well variationally optimized for each individual state as the single-state calculation that was used to determine g 0 . This explains the discrepancy of the numerical and exact ground-state energies seen in Fig. 7, which also provides a convenient order-of-magnitude estimate for the error of the block-improved-relaxation scheme. A shift to high energies for the whole set of numerical values compared to the exact ones of the same approximate magnitude is clearly observed in the data.
In addition, the second and third excited states in  in the occupation number basis of fermions in the harmonic oscillator (h = 0) potential. The two states involve different modes in the fermionic description and thus have different sensitivity to the finite set of available single-particle functions in the block-improved-relaxation method. This is a problem that the proposed rescaling method cannot fix. The problem vanishes when more single-particle functions are used, but this is costly. Single-state simulations indeed show that the energy difference approaches zero for an increasing number of single-particle functions. These problems do not affect the validity of rescaling which is shown to work adequately and to be independent of the shape of the external potential. The energy differences compared to the analytical results range around about 1%.
For h = 10E 0,H a near degeneracy between the third and fourth state emerges. In contrast to the previously discussed degeneracy at vanishing barrier, it is very well described by the simulation. For large enough barrier height the double well resembles two separated anharmonic wells. According to the Bose-Fermi mapping, we identify the degenerate states by the superpositions of (1, 1, 1, 0, . . .); (1, 0, 1, 0, . . .) and (1, 0, 1, 0, . . .); (1, 1, 1, 0, . . .), showing the fermionic configurations in the left-and right-hand wells, respectively. Since these two states involve similar fermionic modes, their degeneracy is well described by the truncated basisset expansion with a limited number of single-particle functions.
Overall, we find good description of the excited states by the rescaled simulation in the double well for a large variation of barrier height h from a harmonic trap (h = 0) to almost completely separated wells.

VI. CONCLUSION
We have proposed a simple rescaling procedure for many-body calculations with a contact interaction in one dimension in a truncated Hilbert space based on Eqs. (11) and (17). The method requires knowledge of only a single parameter. This is found from the interaction strength where the numerical ground-state energy matches the TG energy, which is easily found from singleparticle calculations. For two particles in a ring geometry the rescaling exactly reproduces the correct energies and expansion coefficients of the wave functions. Empirical evidence suggests that the scheme can be extended to multiple particles, external trapping potentials, and lowlying excited states. Testing the method against exact results for ring geometries and harmonic and double-well potentials shows significant improvements in accuracy.
Because rescaling works for excited-state calculations and for arbitrary potentials it is expected to significantly improve time-dependent simulations as well, as long as the state of the system is well described by the low-lying part of the energy spectrum.