Evidence of attraction between charge-carriers in a doped Mott insulator

Recent progress in optically trapped ultracold atomic gases is now making it possible to access microscopic observables in doped Mott insulators, which are the parent states of high-temperature superconductors. This makes it possible to address longstanding questions about the temperature scales at which attraction between charge carriers are present, and their mechanism. Controllable theoretical results for this problem are not available at low temperature due to the sign problem. In this work, we employ worm-algorithm Monte Carlo to obtain completely unbiased results for two charge carriers in a Mott insulator. Our method gives access to lower temperatures than what is currently possible in experiments, and provides evidence for attraction between dopants at a temperature scale that is now feasible in ultracold atomic systems. We also report on spin-correlations in the presence of charge carriers, which are directly comparable to experiments.

Ultracold atomic gases have provided an alternative path to exploring the physics of high-temperature superconductors by emulation of a doped Mott insulator. This system, with charge carriers propagating on an antiferromagnetic background, is widely regarded as the parent state of high-temperature superconductivity [1]. In this setting, the dopants form quasiparticles in the form of polarons, as a result of competition between kinetic energy and super-exchange processes. The carrier can reduce its kinetic energy by delocalizing, but this naturally distorts the spin background. Minimizing the superexchange energy leads to a state with strong antiferromagnetic correlations, which increases the kinetic energy. The polaron emerges as the best compromise.
It is generally believed that pairing in the high-temperature superconductors occurs as a result of an attractive interaction between charge carriers that is mediated by the spinbackground. The precise nature of the microscopic mechanism remains debated, however. A circumstance which complicates this question is that superconductivity does not emerge from a conventional metallic state, but rather, from a non-Fermi liquid [1].
The search for a pairing mechanism has lead to considerable interest in how doping alters spin-correlations in Mott insulators. For example, the resonating valence bond (RVB) theory claims that the spins form a superposition of singlets, providing a better compromise between kinetic and magnetic energy [2]. If such a mechanism is at play, then it should have implications for the spin-background, possibly even above the onset of superconductivity.
Another point of contention is whether fermions form pairs before the onset of superconductivity, such that these merely condense at the critical temperature [3,4]. The conjecture of preformed pairs seems to be invited by the scale of the superexchange parameter J, which in the high-temperature superconductors corresponds to approximately 1500K [5,6]. This raises questions about the temperature range where attraction between charge carriers can be observed. In light of these considerations, spin-and charge-correlations have become the focal point of ultracold-atoms experiments.
Despite intense research, both the pairing mechanism and the nature of the normal state remain open questions, chiefly, as a consequence of a lack of theoretical methods that can reliably address strongly correlated fermions. The Hubbard model-which is believed to capture the essential physics of cuprates-has been investigated with a number of approximative many-body techniques, including density-matrix renormalization-group theory (DMRG) [7], dynamic meanfield theory [8][9][10][11][12][13] and auxiliary field quantum Monte Carlo [14][15][16][17]. Comparing these methods reveals discrepancies, which lead to the conclusion that the phase diagrams of strongly correlated systems are highly sensitive to uncontrolled approximations [18,19]. To date, there are few unbiased numerical techniques for strongly correlated fermions. Numerical linked cluster expansion has provided the equation of state as well as spin-correlations to exceptional precision [20,21], but has thus far not given any insights into polaron physics or superconductivity. Diagrammatic Monte Carlo techniques [22] can resolve pairing in the fermi-liquid regime [23]. Recently, this method has also been extended to strongly correlated systems, though currently published results are limited to relatively high temperatures [24][25][26]. However, for small systems, exact diagonalization and the Lanczos technique do provide indications of attraction between charge carriers [27]. The lack of theoretical methods for correlated fermions has motivated an intense effort to develop experimental techniques with access to microscopic observables. With the advent of quantum gas microscopy, imaging of entangled manybody states is now possible at the level of single-site resolution [28]. Using ultracold atomic gases, Mott-insulating states can be created and cooled to the point where antiferromagnetic correlations become significant [29].
Recent experimental work has focused on spin-correlations in doped systems [30,31]. In particular, this has resulted in the first observation of the internal structure of a polaron on an antiferromagnetic background, revealing microscopic details about the cloud of spins surrounding the carrier [32]. The search for correlations between dopants and thus signs of effective interactions has also been initiated [33].
In this work, we use worm-algorithm Monte Carlo [34] (WAMC) to extract completely unbiased charge and spin correlations in a doped Mott insulator, in the presence of thermal fluctuations. Our findings indicate that attraction between carriers is present at an energy scale corresponding to approximately 700K, significantly above the critical temperature. This is a regime that could realistically be achieved in ultra-cold atomic gases in the near future. We also report on spin-correlations in the presence of two interacting dopants, finding that the delocalization of multiple carriers can explain the reversal of correlators seen in recent experiments [31,33]. WAMC is a technique that allows extremely efficient sampling of world-lines, from which we obtain diagonal elements of the density matrix at thermal equilibrium. Our observables are therefore not subject to any bias, since they are exactly represented by the distribution of world lines. This sets our method apart from previous works, which are either based on uncontrolled approximations, or are confined to small systems, where finite size effects are considerable. WAMC is mainly applied to bosons due to the fermionic sign. However, for Gutzwiller-projected theories, the sign problem is only extensive in the number of charge carriers as opposed to the system size. Using a highly efficient sampling protocol, it is possible to address a small number of dopants in a system that is sufficiently large to avoid finite-size effects. Previously, this method has been used to extract spectral properties of a single carrier [35]. More recently, it has also provided the real-space structure of a single polaron [36].
In this work, we represent charge carriers and spins as world-lines in space and imaginary time. We use a 20 × 20 lattice with periodic boundary conditions to prevent finite-size effects. We use separate worms for the spin and charge sectors [34]. The former can wind in imaginary time, which alters the total spin in the system so that this sector is in the grand canonical ensemble. The worm corresponding to charge cannot wind, keeping the number of carriers to two at all times (in the partition function sector). By generating contributions to the trace of the density matrix at thermal equilibrium, we obtain access to essentially arbitrary operator expectation values, including spin and charge correlators. To confirm the accuracy of this method, we provide benchmarks against exact diagonalization, see supplementary material. We describe the system using the t-J model [37] which effectively captures the low energy physics of the Hubbard model [38] in the case of large onsite repulsion. Herê c iσ ,ĉ † iσ denote annihilation and creation operators of a spin σ electron at the site i. The total particle number on the site i is given byn i = σn i,σ wheren i,σ =ĉ † iσĉ iσ . Two energy scales thus describe the model (1). The kinetic energy is due to hole propagation, and is proportional to t. The superexchange energy originates in the virtual creation and annihi-lation of doublon-hole pairs and is proportional to J = 4t 2 /U , where U is the onsite repulsion.
In the low mobility limit, where t/J → 0, the attraction between two carriers can be understood from a simple brokenbond picture: An isolated hole negates 4 super-exchange interactions. When two carriers share a link, the number of broken bonds is reduced from 8 to 7, saving some magnetic energy. For realistic model parameters, where t/J is not vanishing, this picture becomes overly crude, however. Binding the two holes together impairs delocalization, resulting in competition between kinetic and exchange-mediated attraction. Therefore, it is expected that the binding energy decreases with increasing t/J, to possibly vanish at some critical value (t/J) c [39]. This picture is corroborated by zero-temperature calculations on small clusters based on the Lanczos algorithm, where it is found that with increasing t/J, the binding energy decreases, while the typical separation grows. However, being limited to systems of up to 26 lattice site, these results cannot resolve bound pairs where the distance between carriers is large compared to the system size [27]. DMRG calculations on notably larger clusters-up to 10 × 7 sites-also indicate attraction at zero temperature [40].
At realistic values of t/J, cluster calculations give a relatively small amplitude for finding the holes on neighboring sites, stressing the inadequacy of the broken-bond model [27]. Meanwhile, the delocalization of a charge carrier gives rise to frustration due to competition between kinetic and magnetic processes. This has motivated the suggestion that frustrated bonds may mediate attraction between carriers [40]. Figure 1.
Illustration of spin-spin correlators. The spin-spin correlators C s |d|,s (r) defined in Eq. (3) illustrated in the case of |d| = 1, |d| = √ 2, and |d| = 2. Here, the bond distance r is defined as the distance from one of the carriers to the point between the spins for which the correlation is considered and s is the distance to the second carrier.
To obtain indications of attraction between charge carriers that are directly comparable to experiments, we calculate the hole-hole correlator, which takes the form Heren h r = 1 −n r , is the hole-number operator at site r, s is the distance between the two holes and N is the number of lattice sites. Since the spin background is expected to mediate attraction between carriers, we also calculate spin-correlations of the form Here r 0 and r 0 + s are the positions of the two carriers, and similarly r 0 + r ± d/2 are the positions of the two spins, located a distance d from one another. We will consider the cases |d| = 1 (nearest-neighbor), |d| = √ 2 (next-nearestneighbor), and |d| = 2 (next-next-nearest-neighbor). These correlators are illustrated in Fig. 1.
Hole-hole correlations. (a-c) show contour plots of the correlator C h (s) (Eq. 2) for the case t/J = 2 and βJ = 2.1, 2.4, 2.7 respectively. The radial components of these are given in (g). For this parameter set, we observe a peak in the probability distribution at small separation which is visible already at temperatures corresponding to Jβ = 2.1. In (a) the effect is small, with a peak value of C h ≈ (2.3 ± 0.27) × 10 −2 at s ≈ 4.12. Decreasing the temperature bolsters the effect, and in (c) we observe a maximum value of C h ≈ (6.0 ± 1.7) × 10 −2 , also at a separation of s ≈ 4.12. This scenario corresponds to U = 8t in the Hubbard model, which is the parameter set that was realized in the experiment [33]. Slightly decreasing t/J increases the height of the peak of C h . In (d-f), contour plots of C h (s) are given for t/J = 5/3 and βJ = 2.1, 2.4, 2.7 respectively. The corresponding radial components are again shown in (g). At the lowest temperature (f), the peak in C h is found when the holes are next-nearest neighbors.
In Fig. 2, we show examples of hole-hole correlations for the cases t/J = 2 and t/J = 5/3. These reveal an attraction between the carriers which is manifested in a peak in the correlator C h (Eq. 2) at small separation.
The scenario when t/J = 2 corresponds to U = 8t in the Hubbard model, i.e. where the onsite repulsion is equal to the band width. For this parameter set we observe a peak in C h which is present at Jβ = 2.1 (Fig. 2, a). For comparison, the energy scale of the super-exchange in cuprate superconductors has been estimated to J = 128 ± 5meV, corresponding to approximately 1500K. This suggests an onset of weak attraction at a temperature equivalent to approximately 700K, though it should be stressed that the ratio of the kinetic and magnetic energy scales in the cuprates is somewhat larger, with estimates at t/J ≈ 3.3 [5,6]. This parameter set was recently realized in an optically trapped ultracold atomic gas. Examination of charge-correlations at the temperature Jβ ≈ 1.53 did not reveal any signs of attraction between carriers, consistent with our results [33].
Zero-temperature estimates of the binding energy obtained from the Lanczos algorithm and DMRG-with significant uncertainty due to finite size effects-give results in the span 0.6 ≥ binding /t ≥ 0.15 when t/J = 2. In a temperature range where βt = 4.2 − 5.4, this energy scale is sufficient to impose significant correlations, and the fact that these only reach C h ≈ 6% suggests that thermal fluctuations renormalize the interactions between charge carriers. This can partially be attributed to the loss of the magnetic correlations that mediate attraction. However, it is also the case that spin fluctuations increase the delocalization of the carrier, which also works against pair formation. The latter effect is manifested in the fact that the kinetic energy of a dopant may actually increase as the temperature is lowered [36].
The fact that fluctuations renormalize interactions-even to the point where attraction vanishes-suggests that ground state results provide limited guidance regarding the temperature ranges where we can expect to observe attraction in ultracold atomic gases, where strong temperature effects are ubiquitous.
Decreasing the scale of the kinetic energy to t/J = 5/3 suppresses delocalization. The peak in C h sets in at Jβ ≈ 1.8 and is now larger at comparable temperatures. Eventually, the correlator attains its maximum when the carriers are nextnearest neighbors (Fig. 2, f).
By generating hole-hole correlations for a wide range of parameters and systematically testing whether C h possesses a maximum (within two standard deviations, see supplementary information for details), we obtain the phase diagram shown in Fig. 3. To put this data into context, we have included the model parameters for three recent experiments. The polaron studied in [32] is situated far from the region of attraction, owing to a relatively high temperature and also large onsite repulsion. The experiment on correlations between dopants reported in [33] is, however, closer to the onset of attraction. For this parameter set, our results indicate that a peak in C h appears at Jβ ≈ 2.1. Even lower temperatures are achievable, as exemplified by [29], were long-range spin-correlations were reported at T /t ≈ 0.25. The phase diagram confirms that the attraction is sensitive to the ratio t/J. This is also the case at zero temperature, and stems from competition between delocalization and binding of the carriers [39].  Discs denote parameters for which the correlator C h exhibits a peak at a carrier separation smax, which indicates attraction. The bar intersecting the disc gives the orientation of smax, while the color gives its magnitude according to the legend. The shaded region is a guide to the eye. For comparison, we also indicate the parameter sets of three recent experiments. In Koepsell et al. [32], the internal structure of a single polaron was mapped out using quantum gas microscopy. In Chiu et al. [33], correlations between dopants in an ultracold atomic gas were examined. In Mazurenko et al. [29], spin-correlations in the Hubbard model were examined, though this parameter set does strictly speaking not correspond to a Mott insulator.
In Fig. 4, we show spin-correlators (Eq. 3) for the case t/J = 2 and Jβ = 2.4, which is in the attractive regime. For a single carrier (a, c, e) we observe a cloud of distorted spincorrelations as a result of competition between delocalization and super-exchange, in line with [32,36,41]. This includes a suppression of nearest neighbor correlators (a) as well as the reversal of next-next-nearest neighbor correlations across the carrier (c).
In the proximity of the carrier, the next-nearest neighbor correlators (b) are strongly suppressed due to competing processes: When the dopant hops a single lattice spacing, spins that were previously next-nearest neighbors-and therefore would be correlated-are brought into direct contact, where the interaction is antiferromagnetic. This leads to a highly frustrated state with competition between a kinetic-magnetic interaction and super-exchange, see Fig. 5, (a-c).
When two carriers are placed as next-nearest neighbors, the frustration is alleviated. The primary interaction between the two spins which are situated next to both dopants is now kinetic-magnetic (Fig. 5, d-f), causing remarkably strong anticorrelation as seen in Fig. 4 (d). This result is consistent with magnetic properties, which have been observed in several recent experiments: When doping a Mott insulator, it is found that next-nearest neighbor spin-correlations cross over from positive to negative at a carrier density of approximately 20% [31,33]. Notably, this cross over is not predicted by the structure of a single polaron, and the generation of strong anticorrelations must be understood as a multiple charge-carrier phenomenon. The alleviation of frustration has been identified as a possible source of attraction between carriers in doped  Mott insulators [40].
In conclusion, we present high-precision data for spin-and charge-correlations in a Mott insulator with two dopants at finite temperature, obtained via worm-algorithm Monte Carlo. To the best of our knowledge, this is the first time that unbiased theoretical results are reported for this problem. Our findings indicate that the interactions are renormalized by thermal fluctuations, and that attraction sets in at Jβ ≈ 2.1 when the onsite repulsion is equal to the bandwidth (i.e. U = 8t). This is a temperature range that is now experimentally feasible, and we argue that real progress in understanding this longstanding problem is now possible using optically trapped ultracold atoms and high-precision quantum Monte Carlo simulations in tandem.
Our technique gives access to lower temperatures than what is possible in experiments, allowing us to provide a phase diagram outlining parameter regions where attraction is present. We make detailed predictions about spin-correlation in the In the presence of a single hole (a), the spins S1 and S2 interact indirectly via super-exchange that is mediated by S3. This gives rise to an effective ferromagnetic interaction. However, delocalization of the hole (a-c) brings S1 and S2 into direct contact, where interactions are antiferromagnetic. We refer to this as kinetic-magnetic interaction. The result is a high level of frustration and almost vanishing spin-correlations near the carrier. When two holes are present (d-f), only kinetic-magnetic interaction is present so that S1 and S2 become anti-correlated.
proximity of the two carriers that are directly comparable to quantum gas microscopy and thus serve as a natural bench mark for future experiments.

Worm-algorithm Monte Carlo
Using a continuous-time worm-algorithm Monte Carlo (WAMC), we simulate a pair of holes in the t-J model [1] on a quadratic lattice with 20×20 sites and periodic boundary conditions. The WAMC method effectively samples world-line configurations in real-space and imaginary-time by combining the partition function sector and Green's function sector: In order to go from one partition function world-line configuration to another, one needs to pass through the Green's function sector [2]. We use two different worms to represent spins and holes. The former is allowed to wind in imaginary-time such that the total spin can vary. Contrary, the hole worm is not, and the number of holes in the partition function sector is therefore kept fixed. All data presented are extracted purely from the sampled partition function configurations. In particular, the kinetic energy presented later in this supplementary information is obtained by counting the kinks involving the hole worm [3].

Sign problem
When simulating two holes in the t-J model using WAMC, there are two types of processes contributing to the fermionic sign problem: exchange of indistinguishable spins and exchange of holes. The former enters first at order t 2 J 3 [4] and the latter at order t 4 . When such a process is present in a world-line configuration, the configuration weight p might turn negative. Hence the name "sign problem".
In order to account for negative configuration weights in WAMC, we factor the weight p into its sign s = ±1 and magnitude |p|, i.e., p = s|p|. Then, when computing the acceptance ration, p is exchanged in favor of |p|, and every sampled quantity x is assigned the sign s of the world-line configuration. By performing this substitution, the acquired expectation value x is not that of the original fermionic system. However, the latter one is obtained through x = sx / s [5]. This way of circumventing negative configuration weights gives rise to an exponential decrease in signal-to-noise ratio (SNR) and is the hallmark of the sign problem. The decrease of SNR can, in part, be explained by the exponential decay of the denominator s with -in our case -increased values of tβ and Jβ. This decay of s is illustrated in Fig. 1.

Determining hole-hole attraction
The starting point for determining hole-hole attraction is the hole-hole correlator, which is project onto the separation distance s = |s|. The resulting radial correlation is then partitioned into two parts: the signal part in which the holes are in the close vicinity of one another s < s bg , and the background part where they are further separated s ≥ s bg . The condition for classifying a parameter point as attractive is that the peak value of the signal partition exceeds the background partition's peak value -including two standard deviations of uncertainty (noise). This is illustrated in Fig. 2, and the parameter region displaying attraction is shown in Fig. 3.
The value of s bg is chosen such that the attractive part of the correlation is not included in the background partition. We choose s bg = 6, which exceeds the separation where a peak is observed.
Due to the presence of the fermionic sign problem, the noise grows exponentially with an increased value of tβ and Jβ. Data at these parameter values, therefore, suffer from a poor SNR. Data at small values of tβ and Jβ also have a low SNR, but here the reason being a weak signal. To filter out noisy data, the SNR is defined as (C max − C min ) 2 / n 2 , i.e., the ratio of squared signal amplitude to mean square noise. A data point is deemed too noisy and ignored if SNR < 200.
To reduce noise, the outer edges of the background correlation, where s x , s y = 10 (we performed simulations on systems with 20 × 20 lattice sites), are omitted. These omitted correlation values merely carry half, or even a quarter, of the statistical weight compared to other background values. The SNR of the remaining correlator is, therefore, improved.
To further improve SNR, we ultimately performed largescale simulations, amounting to several million core-hours. The number of sampled world-line configurations for different values of tβ and Jβ is shown in Fig. 4.  Determining hole-hole attraction. Illustrations of how hole-hole attraction is inferred using radially projected holehole correlations obtained using WAMC. If the difference ∆C * = max(C(s < sbg)) − max(C(s ≥ sbg)) − 2 × noise > 0 (a, c), the data is said to indicate attraction, otherwise (b) attraction is deemed to be absent. The partition divider sbg = 6 is chosen such that it exceeds separation distances where a peak is observed.  Figure 3. Attractive parameter region. Parameter region scanned using WAMC. Red dots indicate observed hole-hole attraction, blue dots signify no observed attraction, and a gray cross mark points deemed to have inadequate SNR. The reason for large parameter values suffering from a poor SNR can be explained by the sign problem (c.f. Fig. 1). On the other hand, the low SNR for small parameter values is caused by a weak signal.
Additional hole-hole correlation data In Fig. 5, 6, we present additional data of the hole-hole correlator defined in Eq. (1).    Additional spin-spin correlation data In Fig. 7-9, we present additional data of the spin-spin correlator, defined by Further examples are given as animations [6].

Exact Diagonalization
In order to verify the accuracy of WAMC, it is been benchmarked against exact diagonalization (ED) [7] on a system of 4 × 4 lattice sites with periodic boundaries, containing a Hole-hole correlations.
(a-f) show contour plots of the correlator C h (s) for the case t/J = 5/6 and βJ = 0.90, 1.50, 2.10, 2.70, 3.30, 3.90, respectively. The radial components of these are given in (g). Attraction is present for all of these parameter values. single carrier. In Fig. 10 we present average kinetic energy results, which indicate a perfect agreement between the unbiased WAMC and exact ED method.     Figure 10. Kinetic energy. A single carrier's average kinetic energy, as a function of temperature, in a system of 4 × 4 lattice sites with periodic boundary conditions. The ratio of hopping strength to superexchange strength is t/J = 10/3. While error bars are given for the WAMC data, they are smaller than the marker size. The data indicate a perfect agreement between WAMC and ED.