Quantum Simulation of 2D Quantum Chemistry in Optical Lattices

Benchmarking numerical methods in quantum chemistry is one of the key opportunities that quantum simulators can offer. Here, we propose an analog simulator for discrete 2D quantum chemistry models based on cold atoms in optical lattices. We first analyze how to simulate simple models, like the discrete versions of H and H$_2^+$, using a single fermionic atom. We then show that a single bosonic atom can mediate an effective Coulomb repulsion between two fermions, leading to the analog of molecular Hydrogen in two dimensions. We extend this approach to larger systems by introducing as many mediating atoms as fermions, and derive the effective repulsion law. In all cases, we analyze how the continuous limit is approached for increasing optical lattice sizes.

The field of theoretical quantum chemistry has experienced an extraordinary progress due, in part, to many advances in computational methods [1]. For instance, Density Functional Theory [2, 3] has enabled a better description and understanding of both static [4][5][6][7] and dynamic [8] properties of a large variety of molecules. The capability of such computational methods, whose main challenge is to address electronic correlations, are however sometimes hard to assess experimentally. One approach is to use another (classical) computational technique that is exact in some restricted conditions, but can deal with large systems where exact calculations were not possible. The most prominent example is DMRG [9] which, despite the fact that it operates in 1D lattice systems, offers an ideal platform to benchmark DFT methods [10][11][12][13]. In more general scenarios, the field of quantum computing [14][15][16][17][18][19] can play a key role to overcome numerical limitations in the long-term, offering an excellent setup to benchmark quantum chemistry computational methods. Recently, we have proposed the alternative approach of analog quantum simulation [20], based on the experimentally mature field of ultra-cold atoms [21][22][23], where fermionic atoms play the role of the electrons. While quantum computers and analog simulators would certainly help to push quantum chemistry, the exploration of their full potentiality requires the development of techniques that go beyond the state of the art.
In this Letter we propose and analyze a scheme for analog quantum chemistry simulation that can be im-plemented with present technology. Our approach uses ultracold atoms to address lattice models in two spatial dimensions (2D), where the electron-electron interaction takes different forms. While not exactly reproducing all aspects of the real quantum chemistry scenario, this simulator still retains the most relevant ingredients, enabling the observation of the most representative phenomena in quantum chemistry. Furthermore, it offers a suitable platform to benchmark computational methods in that field. In particular, it allows us to extend the benchmarking offered by DMRG beyond 1D [24].
For the sake of clarity, we will discuss several scenarios, with increasing experimental difficulty, for the simulation of quantum chemistry problems in 2D discrete lattices that could later be compared to contemporary theoretical lattice methods, such as DFT or DMRG. We start with simple one-electron systems, the analogous to the Hydrogen atom, and the H + 2 molecule. Then, we show how to simulate two electron problems, here exemplified by the H 2 molecule. Finally, we show how the system can be scaled-up to more electrons, although with a different dependence of the repulsion with the distance.
Model. In the following, we will consider a discrete version of quantum chemistry models in 2D. First, we start by considering a 2D square optical lattice of size N ×N . N f fermionic atoms, playing the role of electrons, can localize within the local minima of this optical lattice, and hop with nearest-neighbor tunneling rate t F . The Hamiltonian describing their dynamics is then given by: where f † i and f i , are the creation and annihilation operators for a fermionic atom in the i-th lattice site [28], each of them separated by a lattice spacing a, and where  (red), where the nuclear potential is imprinted (blue). For a single simulated electron, this pattern can lead to, e.g., atomic Hydrogen ((a), one nucleus) or H + 2 ((b), two nuclei). For more than one fermionic atom, two different schemes are proposed to mediate an effective repulsion between them. (c) A single atom (green) is used. It tunnels with constant ta through a lattice with the same spacing as the fermionic one. There is an on-site repulsion with strength U when the mediating atom occupies the same site as the fermion. (d) We use as many mediating atoms as electrons need to be simulated (2 in the case of the figure). The on-site repulsion with the fermions now appears in a different internal level b, whose tunneling is slower as compared to level a, using a state-dependent lattice. Both levels are coherently coupled with coupling constant g. the sum is taken over all nearest-neighbor pairs of lattice sites. Fermionic atoms are subject to an external potential that induces the attraction to N nuc nuclei that we consider placed in fixed positions {r n } n=1...Nnuc [29] (Born-Oppenheimer approximation [30]), where Z n is the atomic number of nucleus n, and V (r) is the attractive nuclear potential [31]. In 2D lattices, this potential can be obtained by combining the light shift induced by an external laser orthogonal to the lattice and a fully programmable intensity mask using, for example, a digital mirror device [32]. Depending on the model to be simulated, we will also consider the Hamiltonian H med describing a set of bosonic atoms that mediates fermion-fermion interactions according to some effective potential, V eff . We consider now the simplest situation of simulating atomic Hydrogen. By choosing a potential with a unique nucleus Z 1 = 1 centered in the lattice site r 1 = ( N/2 , N/2 + 1/2), the total Hamiltonian reads as, To begin with, we consider the attractive Coulomb potential on its standard form, V (r) = V 0 /r, for moderate finite lattice sizes, e.g. N = 40. In order to gain intuition, one can compare this discretized Hamiltonian to the continuum limit, where an analytical solution is also known in 2D [33]. As a consequence of the reduced dimensionality, electrons get closer to the nuclei than in the 3D case [34]. Each energy level corresponds to E * n = −Ry (n−1/2) 2 , for n = 1, 2, . . . In that limit, one can also identify, that are the equivalent Bohr radius (a 0 ), and Rydberg energy (Ry), for the 2D discrete model [35]. The first ultimately determines the size of the orbitals and thus how the continuum limit is recovered. In particular, it is needed that the orbitals fit in the lattice (to avoid finite size effects), and that this Bohr radius occupies several lattice sites (to avoid discretization errors), leading to the inequalities, In Fig 2(a) we show the lower part of the spectrum of the discretized Hamiltonian (3) for different values of t F /V 0 and N . First, we observe that we have quantized levels, and thus the discrete model qualitatively reproduces the continuous one. In fact, this can be observed with small lattices (N = 40). Quantitatively, we see that by increasing the ratio t F /V 0 and making the lattice larger, one approaches the continuum limit, as intuitively expected. The error for this approximation as a function of t F /V 0 is shown in Fig. 2 Let us now explore a system with a single fermion and two equal nuclei, Z 1,2 = 1, separated by d/a lattice sites, r 1,2 = ( N/2 ± d/(2a) , N/2 + 1/2), i.e. the analog of H + 2 . This internuclear separation measured in number of lattice sites can be directly expressed in terms of the Bohr radius as d/a 0 = (d/a) · (V 0 /t F ), and therefore compared to tabulated values [37]. In Fig. 3(a) we plot the energy of the ground state as a function of the distance. We obtain a molecular potential, as it is expected for H + 2 , already for the moderate size N = 40. Increasing t F /V 0 favors accuracy, up to the point where finite-size effects appear. At this point the difference in energies to the continuum (dashed line) deviates from the universal scaling ∆E ∝ (t f /V 0 ) −1 , which identifies the optimal configuration for our finite system and a given choice of d/a 0 . In Fig. 3   Two-fermions model. Let us now explore the situation with two fermionic atoms emulating two electrons, where the interelectronic repulsion between them needs to be mediated. For this, we use an additional bosonic atom trapped in an optical lattice potential with the same geometry as the fermions. First, we start with a simple scheme that only considers one of the bosonic internal states, which allows them to tunnel at a rate t a to nearest-neighboring sites. As they coexist in the same lattice sites, elastic scattering processes between the bosonic and fermionic atoms occupying the same position induce an on-site repulsion U , that translates into an effective repulsion between the fermions when the effect of the mediating atom is tracedout: To obtain this expression, we assume to be in the regime in which the bosonic atom dynamics is faster than the movement of the fermions. In this first scheme, and for separations d/a 0.06 e 2πta/U N , this effective repulsion corresponds to, V I (d) ≈ V I,0 /(d/a) , where V I,0 ≈ 6.4e −2πta/U t a (see [36]). This simple scheme then mediates an effective repulsion between the two fermionic atoms that scales as 1/r, matching the dependence of the distance of 3D molecular interactions, but now restricted to 2D [38]. We illustrate the dependence of this potential and its effect in the 2D H 2 molecule in Figs. 4(a-b), respectively. There, one can observe molecular potentials also for relatively small lattices and assess the error. The continuum limit is obtained in a similar regime than the H + 2 molecule case. Many-fermion models: By increasing the number of fermionic atoms in the lattice while maintaining a single mediating boson, one would see that not all interactions among pairs of fermions are equally weighted, precluding scalability. Intuitively, it is more favourable for the mediating atom to localize among the pair of fermions that are closer to each other, rather than in an equal superposition, so that not all interaction are equally considered. In Ref.
[20], this challenge was overcome by including a cavity that symmetrizes these interactions. This cavity interaction is not available in the present, much simplified experimental setup, where interactions are mediated by a hopping atom, instead of a spin-excitation. Another option to induce a pairwise effective repulsion between these fermionic atoms would be Rydberg excitations, that enable for long-range strong atomic interactions. In particular, one can induce dipole-dipole repulsive interactions that depend on their separation as 1/d 3 for distances smaller than the Rydberg blockade radius [39][40][41][42].
Here instead, we present a second scheme that induces pair-wise interactions by including as many mediating bosonic atoms as electrons need to be simulated. This proposal is scalable, at the price of modifying the scaling of the repulsive interaction (see Fig. 1(d)). For these N f mediating atoms, we are going to consider two of its longlived energy levels, that we call a and b, separated by an energy shift ∆. Level b experiences an on-site repulsion U when occupying the same site as a fermion, while the atoms in level a live on a shallow lattice that allows them to move with tunneling rate t a . Both levels are coupled through a Raman (or direct) transition of strength g. Besides, bosonic atoms in the b level suffer an additional hard-core boson interaction |W | |U | which prevents doubly occupied states. The bosonic Hamiltonian then reads as, In particular, we are interested in the regime in which both levels are weakly coupled g ∆, and when the atomic states trapped in the a lattice hop faster than in any of the other levels: t a t b t F (see Fig. 1 . This allows one to trace-out the effect of the mediating atoms and write an effective Hamiltonian for the fermions. By using as many bosonic atoms as fermions, the hard-core boson interactions leads to a bound state in which all fermionic sites are equally occupied, getting a configuration in which the repulsion among each pair of atoms is equally weighted, as required by Eq. (6). For this configuration, the pair-wise mediated interaction scales as, a δII , and δ II = U − 4t a + O g 2 /∆ , (see [36]).
While this system differs from the molecular Hamiltonian observed in nature, it already captures the key features of the interactions appearing in molecular chemistry: nuclear attraction and electronic repulsion. It is then expected to reveal some of the features of chemical systems, including their electronic correlations. In Fig. 4(c), we show the effective repulsive potential induced by the second scheme for different values of detuning δ II , that controls the characteristic length of the interaction. In Fig. 4(d), we illustrate the effect that this modified effective repulsion controlled by δ II has on two fermionic atoms hopping in the lattice, whose dependence on the distance is also mimicked by the tunable attractive nuclear interaction. This leads to a molecular potential of a "pseudomolecule" of hydrogen, where the bonding length and dissociation limit are observed.
Conclusions & Outlook. To sum up, we have shown how ultra-cold atoms moving in 2D optical lattices can be used to simulate simplified models for quantum chemistry in todays experimental setups. We have observed that early experiments with a single simulating atom can pursue the timely goal of simulating the simplest discretized atom and molecule in this platform. In richer scenarios, bosonic atoms can mediate an effective repulsion between the simulated electrons, making repulsive interactions more experimentally accessible with stateof-the-art setups possibilities for further research. First, they provide an experimental platform for which numerical methods used in quantum chemistry can be adapted and benchmarked. Lessons learnt from these simulators, could then be transferred back into improved algorithms for quantum chemistry. Second, one of the main challenges of these discretized 2D simulators is that their solutions approach the continuum result slower than in the 3D case. Fully characterizing this scaling may well lead to improved protocols that are less sensitive to the system size. Third, while this Letter provides strategies to engineer a pseudochemical Hamiltonian in ultra-cold atoms using bosonic atoms as a mediator, other platforms and strategies may also serve for this purpose. Identifying good candidates to simulate specific interactions in chemistry is a promising open field of research.
where me is the mass of the electron and Zn is the atomic number of nucleus n. The first term then describes the kinetic energy of the electrons, the second its nuclear attraction following the potential V (r) , and the third the electronic repulsion. ter. Section A discusses the scaling of the spectrum of the discretized 2D Hamiltonian as the lattice size increases. Section B derives the effective interaction mediated by a single boson with one long-lived state. Section C focuses on the effective interaction mediated by several mediating atoms with two long-lived internal states. Section D includes further details about the numerical calculations shown in Fig. 2 In Fig. 2 we observed that the discretized solutions of the Hamiltonian approached the analytical result following a scaling ∆E ∝ (t f /V 0 ) −1 . This differs from the three-dimensional case, in which accuracy improves as (t f /V 0 ) −2 [20]. To analyze this effect, it is useful to have some insights on how the discretization of the space affects the approach to the continuum solution. A back-ofthe-envelope dimensional analysis can be presented for the 2D case, where we consider the ground-state electronic wave-function, ψ 0 (r) = a −1 0 2/πe −r/a0 . For the two main sources of discretization error, the calculation of the energy terms is based on integrals that are discretized as a Riemann sum. The difference between this sum and the continuum limit is defined to first order by the second derivative of the integrand. For the Coulomb term, this reads as, In the 2D case, this sum does not converge in the continuum limit, and the leading order error corresponds to the diverging term, that is dictated by our choice of the cutoff for the position closest to the nuclei. Normalizing by the Rydberg energy, this error terms scales as (t f /V 0 ) −1 in 2D, and dominates the scaling of the 2D setup as the effective Bohr radius increases, as numerically observed. As an introductory step to gain intuition, in this section we derive how a mediating boson affects the motion of a single fermion by localizing around it. This is the key ingredient responsible for the effective repulsion appearing when more than one fermion are present, that we derive in the next sections. In the limit t F /t B 1, one can make an approximation similar to Born-Oppenheimer. For a single fermion occupying the position j 0 , one can then expand the Hamiltonian H 1B in the basis |j 0 F |φ j0 B , where |j 0 F = f † j0 |0 F and |φ j0 B is the ground state of F j| (H 1B ) |j F where, in the continuum limit, H 1B takes the form, k ω k a † k a k + U a † j a j f † j f j , being ω(k) = −2t b (cos k x + cos k y ) the dispersion relation for a free boson.
In the single fermion subspace, let us start choosing the fermion to be positioned in j 0 = 0 = (0, 0). The eigenstate writes as β † λ f † 0 |0 , where the bosonic operator β † λ = k φ λ (k)a † k . The Schrödinger equation writes as, where φ(j) = 1/N k e −ikj φ(k). In general, for the bound state, describes the single boson localized around the fermion and its bound state energy E B is determined by Its wavefunction writes as, where the normalization factor, We define a pair creation operator F † j , which generates the local Wannier mode F † j |0 = |j F |φ j B by acting on the vacuum state. In terms of F j and F † j , the Hamiltonian under this approximation becomes, by projecting on the bound state energy surface, where the effective hopping strength, of the bound boson-fermion pair is determined by the Franck-Condon coefficient φ j |φ j+1 , i.e., the overlap of the bosonic Wannier states.

Single boson localized around two fermions
By introducing a second fermion, the boson forms a bound-state whose energy depends on this interfermionic separation, inducing an effective repulsion between these two fermions. Aided by the intuition gained in the previous section, here we characterize the properties of this bosonic bound-state.
In the single-boson subspace, the eigenstate writes with parameters, The bound state solution of Eq. (B7) gives rise to the self-consistent equation which determines the relation C 1 = ±C 2 . Focusing on the bound state on the upper-band, that provides the repulsive interaction, and definingk x,y ≡ −π + k x,y , the bound state energy E up corresponds to, This equation encodes how the energy of the bound state depends on the interfermionic separation. Note that d is a 2D-vector with integer components. Equating (B2) and (B11), one gets, The solution to this equation admits a solution given by a recurrence relation on d [44]. Using instead the expansions derived in Sec. B 3 and B 4, one gets for d/a 1/ δ B /t a ), where δ B = E B − 4t a ≈ 2 5 e −4πta/U t a , and γ ≈ 0.577 . . . is the Euler-Mascheroni constant.
This simple model then provides an effective repulsion between the two fermions that scales as δ up (d)/t a ∝ V 0,I /d with V 0,I = 2 7/2 e −γ−2πta/U t a .
From the wavefunction (B3) and the expansion in Sec. B 4 one sees that the characteristic length of the bound states is L I /a ≈ (δ B /t a ) −1/2 . For the previous expansions in (B13) to be valid, one needs to satisfy the regime d/a L I /a. To prevent finite size effects, it is also necessary, that L I /a N . To illustrate this, in Fig. S1 we observe that this expansion for δ B /t a is valid for U/t a > 1, so that L I /a N = 100. In Fig. S2 we also confirm that for this size, the scaling 1/d is maintained for d/a 10, so that d/a L I /a. One can now see that this pairwise interaction does not maintain when more than two fermions are present. To reach this scalability, in Appendix C we will consider a second internal level of the mediating atom.

Calculation of the first integral in (B11)
Defining the energy and length units t a ≡ 1, a ≡ 1 in the coming sections, let us now calculate, One can write an analytical solution [44], In order to extract the scaling of (B11) for frequencies close to the band-gap, it is useful to explore the continuous version of this sum. This will introduce a divergence, that was prevented by the natural cutoff of the lattice. Now, we are interested in the calculation of, for D = [0, 2π] ⊗2 . In the limit kd 1, we can expand the dispersion relation for frequencies close to the upper bandedge, [(k x , k y ) = (π, π)]. Taking the translationk x,y ≡ −π + k x,y , we expand ω(k) ≈ 4 −k 2 , and extend the integration domain to infinite. Note that the numerator e ikd prevents the otherwise divergent integral, and the frequency shift introduces a sign factor, e iπd , that does not enter in the mediated potentials for the strategies presented in this Letter. W.l.o.g., we align vector r in the z-axis, and use spherical units, where K n [x] is the modified Bessel function of the second kind [45] and d ≡ |d|. For small arguments (0 < x 1), Appendix C: Mediating atoms with two long-lived states When more than two fermionic atoms are introduced, the effective repulsion mediated in the previous section by the single-boson bound-state is not purely described by the pair-wise separation between each pair of fermions. To gain this feature, let us introduce in this Section a modified scheme, where we consider two internal levels of as many mediating atoms as fermions there are in the system. We will denote the two levels as b and a. Atoms in b level experience an on-site repulsion when occupying the same site of a fermion, while atoms in state a live on a shallow lattice that allows them to hop with tunneling rate t a . Both levels are coupled through a Raman transition of strength g and are shifted by energy ∆. In order to equally account for repulsion among each pair of fermionic atoms, we include an on-site repulsion W among them when they occupy the same lattice site, obtaining the mediating Hamiltonian, Intuitively, mediating atoms localize around the fermionic positions, and double occupations are prevented by the hard-core boson interaction W U . This then creates a bound-state in which each mediating atom localizes in a different fermionic position. As compared to the previous scheme, hopping from one fermion to the others now becomes a fourth-order process in the coupling g between the two atomic metastable states, as the movement of two mediating atoms is needed.
In particular, we are interested in the regime in which both levels are weakly coupled g/∆ 1, and atoms in level a hop in a lattice much more shallow than the rest: t f t b t a . As it occurred in the previous case, this last inequality allows to trace-out the effect of the mediating atom, writing an effective Hamiltonian for the fermions, ij V (|i − j|)f † i f i f † j f j . Let us now derive this regime using perturbation theory for g/∆ 1 and N f fermions occupying fixed positions j 1 . . . j N f . For this, let us separate the bosonic Hamiltonian (C1), as In particular, we are interested in the energy correction of the bound-state |ψ B,II = , that depends on the interfermionic positions. For this, we need to expand the perturbed Hamiltonian. One can see that only even orders enter the calculation, and expanding to fourth order, one gets the equation, This latter term originates from the pairwise repulsion introduced by the fourth-order correction of two mediating atoms swapping the fermionic position they localize around. This then leads to an effective pairwise potential, These two independent sums can be calculated as in Sec. B 4. Note that the alternating sign derived in Sec. B 4 cancels after the double product e ikd e −iqd .
, which, to lowest order in the regime d √ δ II > 1, scales as, This then leads to a pairwise repulsion between the fermionic atoms that decays exponentially with their separation, following a decay length L II ≡ 2 √ δ II −1/2 . In where E (2) B,II (δ) approximates the second order correction in (C4) as, that can be expanded as in (B15).
Appendix D: Numerical methods

Exact diagonalization
Once the kinetic term is approximated as a nearestneighbor hopping term (1), the Hamiltonian can be conveniently written in a position basis and the ground-state obtained using exact diagonalization (ED).
In Fig. 2 we use this approach to calculate the energies λ n associated to the lowest part of the spectrum of Hamiltonian (3) for different choices of the ratio t f /V 0 . These energies are shifted to correct the shift induced by the nearest-neighbor approximation, ω(k) ≈ 1 − k 2 /2. The result is divided by the Rydberg energy; this is, where, in this case N f = 1. The same strategy is applied to calculate the fermionic potential in Fig. 3(a), where only one mediating atom is involved.
This approach is also used in Fig. 3(a) to calculate the ground-state energy of an Hydrogen cation for a given internuclear separation d/a 0 . In addition to the previous shift, nuclear repulsion V 0 /(d/a) needs to be included before expressing the result in Rydberg energies. Similarly to the atomic case, for a fixed interatomic distance, accuracy improves by increasing the effective Bohr radius a 0 /a = t F /V 0 , up to the point in which finite-size effects become relevant. The number of lattice sites separating the nuclear positions d/a is then adjusted accordingly, identifying the optimal separation value as the one giving the lowest ground-state energy (see Fig. 3(b)).
The same strategy is also applied to obtain the ground state energy of H 2 in Fig. 4(b). The main difference is that now N f = 2, and further simplifications can be made taking into consideration the fermionic statistics. As each fermion can occupy N 2 sites which, together with the fermionics statistic {f i , f † j } = δ i,j leads to a Hilbert space of space of size N 2 (N 2 − 1)/2. ED is also used in 4(b) to calculate the effective potential mediated by a single bosons, for fixed fermionic positions separated by d/a sites and centered in the lattice.
The exponential decaying potential explored in Fig.  4(d), requires a more careful analysis, as the natural rescalings to the Bohr radius and Rydberg energy does not apply now. In particular, three parameters can be independently tuned: the fermionic hopping t f , the interacting potential V 0 and the decay length L II /a. As compared to the previous case, one can remap V 0 → V 2 0 /t f and L II → L II t f /V 0 , so that the final result is still dimensionless when normalizing the energies by the previous definition of Bohr radius V 2 0 /t f . As an illustration, in this Figure 4(d), L II /a = 5 is chosen, and t F /V 0 is fixed as the ratio providing maximum accuracy for the atomic case (one fermion an one nuclei) hopping in a lattice of side N/2 , so that the dissociation limit is properly captured. Modifying the separation d/a between nuclear positions then allows to scan the different internuclear separations d/a 0 = (d/a) · (V 0 /t f ) for this fixed value of t F /V 0 .

Imaginary time evolution
For the calculation of the effective potential mediated by the two metastable levels of atoms in 4(c), we use Imaginary Time Evolution (ITE). This is a useful strategy to numerically obtain the ground state of a gapped Hamiltonian with purely positive eigenvalues, and consists on iteratively evolving an initially random state as, e −H·t . After each iteration the resulting state is normalized, and the contribution of the excited states is mostly reduced.
In more detail, one of the advantages of this method is that rather than writing the entire evolution operator O(N 4 ), one can choose to work in a diagonal basis, so that only O(N 2 ) terms are needed to describe the state at each point in time. From the computational perspective, this is specially useful when facing the multielectronic case. In principle, to calculate the interaction among N b -bosons one would need a state with (N 2 ) N b entries. Using ED, one would need to write the Hamiltonian, of size (N 2 ) N b × (N 2 ) N b . In contrast, evolving the state in imaginary time evolution only needs to store the diagonal terms [with size (N 2 ) N b ], once the state is expressed in a basis that commutes with the terms of the Hamiltonian. For our particular case, this corresponds to the position representation for the on-site interactions, and momentum representation for the kinetic term. The Hamiltonian H nuc is already diagonal in position basis, and one can define a momentum basis, where H K reads as being ω k,f = −2t F (cos(k x ) + cos(k y )) the dispersion relation. This induces a periodic boundary condition in the lattice, which does not affect the calculation as long as finite-size effects are prevented. To confirm that is the case, for each choice of parameters we check that the same result is obtained for the single-boson case using ED, evidencing that boundary conditions are not affecting the result.
To calculate the ITE of Hamiltonian (8), a constant energy shift is added to H during the calculation to make all the spectrum positive, which is later subtracted at the end of the calculation. To evaluate the operation, ψ(t) = e −Ht ψ(0) we use a Suzuki-Trotter [46] expansion of the first kind, dividing the evolution in n steps as e −Ht ≈ n−1 k=1 e −H∆t + O (∆t), and t k = k · ∆t/t. For each of these steps, we calculate e −H∆t ψ(t k ) ≈IFFT e −H K ∆t FFT e −HR∆t e −Hg∆t ψ(t k−1 where (I)FFT indicates the (Inverse) Fast Fourier Transformation, and normalize the resulting state. Here, H R denote the terms that are diagonal in the position basis, and H K the ones in momentum basis. H g denotes the coupling term, whose exponential can be directly calculated noting that, e −g(a † j b j +H.c.)∆t = cosh(g∆t)(a † j a j + b † j b j ) − sinh(g∆t)(a † j b j + H.c.). We iterate this procedure until the overlap between ψ(t k−1 ) and ψ(t k ) is smaller than 10 −5 . We initialize the algorithm with a random state for the smallest value of t F /V 0 , and use this converged solution as the initial state for the next configuration of t F /V 0 .
In our second scheme, N f atoms with two long-lived states are used to mediate the interaction among N f fermions. For a given fermionic configuration, we desire to numerically calculate the bound state, and compare it to the analytical expansion previously introduced in Eq. (C5). For this calculation, we use the ITE method (Sec. D 2), where now, each of the N f mediating atoms can occupy any of the 2 levels at any of the N 2 lattice sites, which a priori accounts for states of size (2N 2 ) N f . To reduce this space, we assume that |g/(U − ∆)| 1, so that level b is only populated in the sites where they interact with the fermions, and |W/U | 1, so that two mediating atoms in level b do not coexist in the same lattice site. For a configuration of 2 (3) fermions in sites r, s(, t), and given the indistinguishability of the mediating atoms, we can further reduce the Hilbert space to states written in the basis collected in Tables I and II. Within this basis, in Fig. 4 we calculate how the energy of the bosonic ground-state energy E(d) depends on their separation d between two fermionic atoms fixed in lattice sites [N/2−d/2, N/2] and [N/2+d/2, N/2], following the same strategy used in the previous case for N f = 2.
In the case N f = 3, we observe that the biggest demand on computational memory corresponds to describing processes in which the three mediating atoms simultaneously populate the a-level. Such processes scale as [g/(U − δ)] 6 in perturbation theory, and are subleading when compared to the second-order terms. Therefore, truncating 0 ≤ N a ≤ 2 would allow to push the calculation at a marginal error (see Table II).
To confirm this intuition, in Fig. S3(a) we use ITE to calculate the bosonic bound state for 3 fermions describing a triangular isosceles configuration. For moderate sizes (N = 16), we compare the numerical result given by this truncated space to the one obtained for the total basis. As desired, we observe that (1) the truncation to the space with up to 2 excitations in state a does not modify the solution, and (2) the scaling is in agreement with the calculation for a pairwise repulsion given by (C6).