Spectral density of interacting pairs in low-dimensional finite lattices

The spectral density of bound pairs in ideal 1D, 2D and Bethe lattices is computed for weak and strong interactions. The computations are performed with Green's functions by an efficient recursion method in real space. For the range of interaction strengths within which bound states are predominantly single pairs, the spectral profiles guide to the energy bandwidths where the bound pairs can be maximized.


Introduction
For the Hubbard model [1] in lattices, formation of composite objects can be observed for both attractive and repulsive interactions [2]. A doublon is such a composite object formed from two fermions. In Auger spectroscopy, certain lines and lineshapes in the spectra can be explained within the formalism of formation of bound pairs [3][4][5]. In general, pairs formed from interacting particles show the properties different from constituent individual particles. For strongly interacting particles, their lifetimes are generally high [6] and depend on interaction strengths. These strongly correlated pairs are now realized and observed in optical lattices within broad range of tuneable interaction strengths [7]. Recently, the pairing phenomena have been also observed in optical lattices for spin excitations [8] which are relevant for excited states near the ground ferromagnetic state. The bound states that arise in lattices from onsite and nearest neighbour interactions are also related to solitons [9,10].
Two interacting particles show nonclassical correlations depending on the statistics of the particles [11]. The interaction changes the effective statistics between particles. However, a clear mathematical account on this is yet expected to emerge. The constraints of the two-particle density matrices are also not known [12]. One can study the correlations from quantum walk experiments and simulations which are fundamental building blocks of many-body systems. The dynamics of particles in lattices changes significantly depending on the binding. The quantum walk of correlated particles has promising applications in quantum technologies [13][14][15].
A general understanding on the mechanism of superconductivity is also still unresolved. The Hubbard model with onsite and nearest neighbour interactions is mostly used for explaining the phenomena. Pairing of fermions with attractive interaction is known to cause the mechanism [16,17]. In two-dimensional Hubbard model, the super-exchange in quasi two-dimensional oxides and nearest neighbour coulomb interaction between polarons in organics [18] is known to be responsible for the phenomena. It is suggestive that an unified approach to the mechanisms may have relations to the two-particle bands with weak and strong interactions within the range where two-particle bound pairs have a higher density than the bound states with a larger number of particles.
The effects of magnetic fields on single particles or excitations in two-dimensional lattices are well known from Hofstadter's approach [19]. The density of states for two particles for various interaction strengths had also been computed before [20]. In optical lattices, the external magnetic fields can be synthesized with periodic modulation of the lattice potentials [21][22][23]. This article specifically looks at the density of bound pairs for both zero and non-zero magnetic fields at different interaction strengths. In the context of condensed matter, the effective fields can also be realized through doping with isotropic magnetic impurities.
The solution for the two-particle states are known exactly for strong interaction cases which are described with continuum and bound states [24]. Within many particle environment, there can be several bound states and interplay of these bound states with the continuum states are highly complex. These many-particle bound states can appear at different energy windows with different densities. Green's functions provide a direct approach to study these states at different interaction strengths. Recursion methods [25] are known for efficient computation of these Green's functions. Several other methods have also been applied recently for such computations, e.g., exact diagonalization with the method of Lanczos [26], DMRG [27] and so on. In this article, the Green's functions are computed using Berciu's approach [28] which is best described in the original article. A brief description is given in the method section for the reader. With Berciu's method, the Green's functions can be computed efficiently for the cases of zero and non-zero external magnetic fields as well as for disordered systems [29]. The bound pairs can be found in general graph architectures where the quantum transport may prove relevant to quantum technologies. The Bethe lattice is one such graph for which the density of bound pairs at different interaction strengths is considered in this article. The Bethe lattice which is important for its mathematical form, also has presence in cores of several photosynthetic complexes.

Method
In this article the density of a bound pair of two interacting spinless particles is studied. The particles are taken as hardcore bosons in 1D and 2D lattices and softcore bosons in Bethe lattice. The case of hardcore bosons is similar to two spinless fermions. In 1D and 2D lattices, the particles get bound with nearest neighbour interactions. For the models with zero nearest neighbour interaction and non-zero onsite interaction, the same bound pairs appear for spin paired particles. The simplified model Hamiltonian is that of extended Hubbard model H H with the hopping limited to the nearest neihbour sites on the lattices.
The onsite interaction term is taken as infinite in 1D and 2D cases. The nearest neighbour interaction term is neglected for the Bethe lattice. Onsite energies are scaled to zero for regular lattice systems. For a Hamiltonian H , the Fourier transformed Green's function is defined as where = + i is a complex number with a very small positive real number. In real space, the two particle Green's function 2 ( , , ) = | ( )| is the propagator for two particles from sites , to sites , in a lattice, where | = † † |0 and the initial site indices , are omitted for brevity. The bound pair spectral density can be calculated for the bound state from the Green's function 2 ( , , ) in real space with initial distance | − | = 0 and 1 for softcore and hardcore bosons, respectively The total density of states (DOS) is obtained with summing over densities of all the states On the graphs with translational symmetry, a few states with increasing relative distance (| − |) prove sufficient for converging results.
Substituting equation (2.1) in equation (2.2) generates the recursion relations between the Green's functions for every site indices , of the two particles. The Green's functions with constant + are then grouped with a vector V ( = + ) which forms a 1D chain and generates the recursion equation with C ≠ 0 for = + = containing the initial state. The hopping matrices α , β connect Green's functions between nearest neighbour vectors with the hopping integrals. Solving this equation, any Green's function can be computed for a finite lattice. At the boundaries, equation (2.6) becomes where 0 and are the minimum and maximum indices for . The recursion modifies to this general form These A and B matrices can be computed recursively starting from equation (2.7) before one reaches = from both sides of the chain with Once V is found, all other V are given by equation (2.8). The procedure accounts the full self-energy term which can be obtained from renormalized perturbation expansion [30] on any ordered or disordered lattice. This method can also be used for the cases where sites have internal structures.
In presence of external magnetic field, a gauge field appears in the Hofstadter Hamiltonian with the hopping terms which is known as Peierls substitution [31].
For an ideal 2D lattice, the onsite energy terms can be scaled to zero with = ∞ in the phase for zero magnetic field. The phase term contributes with the hopping on one of the axes of 2D lattice for non-zero magnetic field. The subscripts 1, 2 of the site index denotes coordinates of two axes of site . These terms are simulated in optical lattices with periodic modulation of lattice potentials [22]. The essence of the method depends on the classical Floquet theory where the momentum terms gain phases.

Results
The spectral density of bound pairs can be obtained from Green's functions computed using a recursion method or Chebyshev polynomials or exact diagonalization. The density of states had been computed before by Rausch et al. [32], Halpap et al. [33] and others for one dimensional regular lattices. The density of states for two-dimensional ideal lattices had been computed before by Barelli et al. [20]. I use the method of recursion that has been formulated recently [28]   The bound pair spectra in 1D have been recently studied by Rausch et. al for strong interactions. I obtained the results for weak interactions using the recursion which gives converging results at the band edges. The bound pairs within continuum can have different transport properties because of the scattering with the continuum states which I do not discuss in this article. The density of bound pairs at = 0 is maximum near = . The density of pairs at = 0 for noninteracting case is much smaller than for interacting ≠ 0 cases. The bound state splits from the continuum at = /2 where is the bandwidth of the two-particle continuum spectra. The bound state shows a sharp peak near = for weak interaction strengths ( < /2). For strong interaction cases, the density reflects the cosine dispersion curve of bound states [24] near = . The total density of states is similar to the density of states that is obtained analytically for the noninteracting single particles in 2D lattices [30] except the satellite peaks which are clearly observed for the strong interaction cases.
The bound pair density shows a similar structure of one bound state in 2D lattices in presence of interaction. The density profiles are very similar to that of 1D lattice for = 0. For weak interaction strengths, the pair density reaches the maximum within a range of the bandwidth rather than a sharp peak and decays to the band edge afterwords. For strong interactions, the pair density shows a sharp maximum at > /2. In presence of external gauge fields, splitting in the spectra can be observed from figure 2. The details of the spectra then depend on the flux per plaquette . At half flux per plaquette for the non-interacting cases, the pair density is reminiscent of Kondo profile. At strong interactions with half flux per plaquette, the density is regular dome-shaped while at non-half fluxes, the density is highly irregular. At weak interaction, the maximum density is shifted to higher . These density profiles will be reflected symmetrically about = 0 with changes of sign in the interaction term → − .
The bound pair state shows non-trivial features on Bethe lattices. At zero onsite interaction = 0, the spectra show discontinuity and multiple sharp peaks. At strong interactions, these peaks merge within a single envelope. The pair density in figure 3 is obtained for two bosons, and the interactions between particles at different nodes are neglected. The initial states are prepared at the root layer = 0 with two bosons on the node ( = 0, = 1). For two particles, the recursion can be performed with = 1 + 2  adding layer indices for the two particles on the Bethe lattice. For single particles, the recursion vectors map directly to the branches of each layer divided as and as depicted in figure 3. The Green's functions can be then efficiently computed for any initial superposition of states on these vectors, e.g., the rectangular box on the graph will be the vector V 1 in the recursion. The discontinuities of the spectra for weak interactions within small energy width suggest tunable control of bound pair transport on these graphs.

Conclusion
In this article the doublon spectra of regular 1D, 2D and Bethe lattices have been computed from real space Green's functions. In two dimensions, the effects of external magnetic fields on the spectra have been computed for weak and strong nearest neighbour interactions between two hardcore bosons. For weak interactions, the gauge fields enhance the density at few bandwidths. In Bethe lattices, the spectra show discontinuities for weak interactions within the band. This work points to energy windows where maximum bound pairs can be found in 1D, 2D and Bethe lattice systems with the change of interaction strengths between particles. The work might be helpful for understanding the many-particle states where bound pairs can lead to the difference in responses, e.g., conductivity and spectroscopy of the lattice systems. The scattering and dissipation of the bound pairs within continuum should also be analyzed. This work computes the spectral profiles which provide insights into the mechanism of interacting particles when two particle bound pair densities are much higher than the bound states with a larger number of particles. Many-particle bound states with a larger number of particles at different energies should be studied further which can also be important for the systems with strong interactions.