Bond Bipolarons: Sign-free Monte Carlo Approach

Introduction. Broad interest in polarons—stable quasiparticles emerging as a result of strong renormalization of properties of a bare particle due to interaction with this or that environment—comes from their ubiquitous presence across all fields of physics. There exist electron-phonon polarons Landau33 ; Frohlich50 ; Feynman55 ; Schultz1959 ; Holstein59 ; Alexandrov99 ; Holstein2000 ; Mona2010 , spin-polarons BR70 ; Nagaev74 ; Mott06 , Fermi-polarons Lobo06 ; Bulgac07 , protons in neutron rich matter protons , etc. Coupling to the environment may also lead to formation of bound particle states called bipolarons. Their studies are motivated by the possibility of identifying a possible mechanism for high-temperature superconductivity when a gas/liquid of bound pairs of fermions undergoes a superfluid (SF) transition at temperature

$$T_{c}\sim{n^{2/d}\over m_{\rm*}}\,,$$

(1)

where $n$ is the bipolaron density in a $d\geq 2$-dimensional system and $m_{\rm*}$ is the bipolaron effective mass. The estimate for $T_{c}$ is very robust against repulsive interactions between bosons: in $d=3$ it barely changes between the gas and liquid densities Ceperley ; Landau ; tc3da ; tc3daa ; tc3db , and in $d=2$ the dependence on interactions is logarithmically weak tc2da ; tc2daa ; tc2db ; tc3db . Since the SF transition temperature increases with density, the highest value of $T_{c}$ for this mechanism is obtained when the inter-particle distance is of the order of the pair size, or $nR_{*}^{d}\sim 1$, where $R_{*}$ is the bipolaron radius. By increasing the fermion density further one goes through the so-called BEC-to-BCS crossover corresponding to a radical change of the microscopic mechanism behind the SF transition—the (quasi) Bose-Einstein condensation (BEC) of spatially separated pairs gets gradually replaced by the Bardeen-Cooper-Schreiffer (BCS) pairing in momentum space—and $T_{c}$ starts decreasing exponentially Leggett ; BCSBEC . Thus, the maximum value of transition temperature can be expressed though the bipolaron parameters as

$$T_{c}\sim{1\over m_{\rm*}R_{*}^{2}}\,.$$

(2)

When the dominant mechanism of electron-phonon coupling is of the density-displacement type, as in the Holstein model Holstein59 , large values of $T_{c}$ cannot be reached Chakraverty ; PhysRevLett.84.3153 ; PhysRevB.69.245111 . The reason is exponentially large effective bipolaron masses originating from small phonon overlap integrals for realistic values of the adiabatic parameter

$$\gamma\,=\,{\omega_{\rm ph}\over W}\,\ll\,1,$$

where $\omega_{\rm ph}$ is the phonon frequency and $W$ is the particle bandwidth; $W\approx 4dt$ for tight-binding dispersion relation with hopping amplitude $t$ between the nearest-neighbor sites on a simple $d$-dimensional cubic lattice. The corresponding bipolarons are also very sensitive to the on-site Hubbard repulsion $U$. Much stronger electron-phonon interaction (EPI) is required for their formation when $U\sim W$, leading to an additional exponential increase of the effective mass PhysRevLett.84.3153 .

The prospects for obtaining high $T_{c}$ appear to be far better when strong EPI originates from hopping dependence on atomic displacements by one of the two mechanisms. In mechanism A, tunneling is enhanced (suppressed) when the distance between the sites is reduced (increased) PhysRevLett.25.919 ; PhysRevB.5.932 . In mechanism B, tunneling is modulated by displacements of atoms in the barrier region between the sites KK . Previous studies Mona2010 ; us21 ; Sous21 found that these bond polarons remain relatively light when entering the strong coupling regime. However, properties of the corresponding bipolarons remain unexplored: the only available study Sous2018 considered the $d=1$ case in the antiadiabatic limit $\omega_{\rm ph}=3t$ for mechanism A. For realistic predictions of high $T_{c}$ one needs to look at higher dimensional systems in the adiabatic limit $\gamma\ll 1$ and account for large repulsive potential between the electrons. This poses a significant computational challenge for unbiased methods based on exact diagonalization PhysRevLett.84.3153 ; Bonca2001 or controlled truncation of the phonon Hilbert space Carbone2021 because bipolaron states in $d>1$ are extended and the number of excited phonon modes quickly increases with $1/\gamma$.

In this Letter, we show that path-integral representation for particles offers a unique opportunity for conducting comprehensive studies of $n$-particle polaronic states when EPI is based on mechanism B. The simplest Hamiltonian can be formulated as $\hat{H}=\hat{H}_{e}+\hat{H}_{\rm ph}+\hat{H}_{\rm int}$ with

$$\hat{H}_{e}=-t\!\!\sum_{<ij>,\sigma}\!\!(\hat{c}_{j\sigma}^{\dagger}\hat{c}_{i\sigma}^{\;}+{\rm H.c.})+\!\sum_{ij\sigma\sigma^{\prime}}\!V_{\sigma\sigma^{\prime}}(i,j)\hat{n}_{i\sigma}\hat{n}_{j\sigma^{\prime}},$$

(3)

$$\hat{H}_{\rm ph}=\omega_{\rm ph}\sum_{<ij>}\,(\hat{b}_{<ij>}^{\dagger}\hat{b}_{<ij>}^{\;}+1/2)\,,$$

(4)

$$\hat{H}_{\rm int}=-g\!\sum_{<ij>,\sigma}\,(\hat{c}_{j\sigma}^{\dagger}\hat{c}_{i\sigma}^{\;}+{\rm H.c.})(\hat{b}_{<ij>}^{\dagger}+\hat{b}_{<ij>}^{\;}),$$

(5)

where $V_{\sigma\sigma^{\prime}}(i,j)$ is the potential of interelectron interaction and $g$ is the strength of EPI. We use standard notation for creation and annihilation operators for electrons (on site $i$ with spin $\sigma$) and optical phonons (on bonds $<ij>$). The model can be further generalized to deal with several dispersive phonon modes and study effects of phonon dynamics and polarization. Here we focus on describing the numerical method and present illustrative results for bipolaron states in a two-dimensional system. These results indicate that the search for high-$T_{c}$ regimes requires comprehensive exploration of the model parameter space because compact states can easily end up to be heavy while light effective masses come at the price of larger pair radius, see Eq. (2). One counter-intuitive effect is that the $s$-wave transition temperature may increase with the on-site repulsion $U$ because less compact bipolarons are significantly lighter.

Basic relations. For a generic few-body system, the mathematical object containing all information about ita ground-state properties and potentially allowing sign-free Monte Carlo simulations can be formulated as

$$O_{ba}(\tau)=\langle\,b\;|e^{-(\tau/2)\hat{H}}\,\hat{O}\,e^{-(\tau/2)\hat{H}}|\,a\,\rangle\,.$$

(6)

Here $\hat{O}$ is a certain observable, $\hat{H}$ is the system’s Hamiltonian, $\tau$ is an appropriately large imaginary-time interval, $|\,a\,\rangle$ and $|\,b\,\rangle$ are any two states having finite overlap with the ground state, $|{\rm g}\rangle$. In the asymptotic limit of $\tau\to\infty$ the answer is dominated by the projection onto the ground state when $O_{ba}(\tau)$ is given by the product of the ground-state expectation value of $\hat{O}$ and the universal (for all observables) propagator ${\cal G}_{ba}(\tau)$:

$$O_{ba}(\tau)\;\mathop{\longrightarrow}\limits_{\tau\to\infty}\;\langle\,{\rm g}\,|\,\hat{O}\,|\,{\rm g}\,\rangle\,\langle\,b\,|e^{-\tau\hat{H}}|\,a\,\rangle\equiv\bar{O}_{\rm g}\,{\cal G}_{ba}(\tau)\,.$$

(7)

If the states $|\,a\,\rangle$ and $|\,b\,\rangle$ are free of phonons, ${\cal G}_{ba}$ is the standard $n$-particle Green’s function. Since ${\cal G}_{ba}(\tau)\equiv I_{ba}(\tau)$, where $\hat{I}$ is the identity operator, the expectation value of $\hat{O}$ in the ground state can be represented as

$$\bar{O}_{\rm g}\,=\,\frac{O_{ba}(\tau)}{I_{ba}(\tau)}\,\equiv\,\frac{\sum_{ab}W_{ab}\,O_{ba}(\tau)}{\sum_{ab}W_{ab}\,I_{ba}(\tau)}\,.$$

(8)

This is standard for MC methods setup when the stochastic sampling process is designed to sample the propagators ${\cal G}_{ba}$ while matrix elements of physical properties are taken care of by Monte Carlo (MC) estimators. The extended version of the relation—with sums over any subsets of states $|\,a\,\rangle$ and $|\,b\,\rangle$ with arbitrary weights $W_{ab}$—adds flexibility and efficiency in designing updating schemes.

According to Eqs. (6) and (7), the imaginary-time dependence of ${\cal G}_{ba}$ also contains direct information about the ground-state energy, $E_{\rm g}$:

$${\cal G}_{ba}(\tau)\;\mathop{\longrightarrow}\limits_{\tau\to\infty}\;\langle\,b\,|\,{\rm g}\,\rangle\langle\,{\rm g}\,|\,a\,\rangle\,e^{-\tau E_{\rm g}}.$$

(9)

Moreover, by selecting the state $|\,b\,\rangle$ or state $|\,a\,\rangle$, or both, to belong to a particular—non-ground-state—symmetry sector, we can employ (9) to determine energy of the ground state in a given sector. This way ${\cal G}_{ba}(\tau)$ can be used to obtain the quasiparticle dispersion relation, and, in particular, its effective mass.

A direct procedure of extracting $m_{\rm*}$ from ${\cal G}_{ba}(\tau)$ is based on the coordinate representation when for each of the states $|\,b\,\rangle$ and $|\,a\,\rangle$ we introduce the notion of the “center-of-mass” positions, ${\bf R}_{b}$ and ${\bf R}_{a}$, respectively, and consider the relative-coordinate dependence of the propagator ${\cal G}_{ba}(\tau,\,{\bf R})$, where ${\bf R}={\bf R}_{b}-{\bf R}_{a}$. In the limit of long time and large distance, $(\tau,{\bf R}^{2})\to\infty$, this dependence takes on the characteristic form of a single free particle propagator with effective mass $m_{*}$:

$${\cal G}_{ba}(\tau,{\bf R})\,\to\,{A_{ba}e^{-E_{\rm g}\tau}\over\tau^{d/2}}\,e^{-{m_{*}{\bf R}^{2}\over 2\tau}}.$$

(10)

Apart from the coefficient $A_{ba}$, the r.h.s. of (10) is insensitive to the particular choice of states $|\,a\,\rangle$ and $|\,b\,\rangle$, allowing one to average over them as in Eq. (8). In the asymptotic limit, the difference between the lattice and continuous space disappears, and $m_{*}$ is directly related to the mean square fluctuations of the relative coordinate:

$$\overline{{\bf R}^{2}}(\tau)=\frac{\sum_{ab}W_{ab}\,{\cal G}_{ba}(\tau,{\bf R})\,{\bf R}^{2}}{\sum_{ab}W_{ab}\,{\cal G}_{ba}(\tau,{\bf R})}\;\mathop{\longrightarrow}\limits_{\tau\to\infty}\;\frac{d}{m_{*}}\tau\,.$$

(11)

Starting with zero at $\tau=0$, the $\overline{{\bf R}^{2}}(\tau)$ function ultimately saturates to a straight line at long $\tau$ leading to an accurate estimate of the effective mass from its slope.

Perturbative expansion. Our approach to MC sampling is based on the general scheme proposed in Ref. worm and rendered here. Let $\hat{H}_{0}$ and $\hat{V}$ be the diagonal and off-diagonal parts of the Hamiltonian with respect to some basis ${\cal B}=\{|\alpha\rangle\}$ (out choice is site Fock states for particles and bond Fock states for phonons): $\hat{H}_{0}|\alpha\rangle=E_{\alpha}|\alpha\rangle$, $\langle\alpha|\hat{V}|\alpha\rangle=0$. Next, we decompose $\hat{V}$ into elementary non-vanishing terms, $\hat{V}=\sum_{\alpha\beta}V_{\beta\alpha}|\beta\rangle\langle\alpha|$. In the chosen representation, there are three types of the elementary terms each corresponding to the particle hopping along a specific bond in a specific direction: (i) bare hopping (ii) hopping assisted by the phonon creation, and (iii) hopping assisted by the phonon annihilation. For a sign-free MC method—apart from, possibly, the negative signs (or phases) originating from components of the vectors $|\,a\,\rangle$ and $|\,b\,\rangle$ in the basis ${\cal B}$—all matrix elements $V_{\beta\alpha}$ need to be non-negative real numbers. This requirement is satisfied with our choice of ${\cal B}$ for model (3)-(5).

Using standard interaction representation for the evolution operator in the imaginary time domain, $e^{-\tau\hat{H}}=e^{-\tau\hat{H}_{0}}\hat{\sigma}(\tau)$, and expanding $\hat{\sigma}(\tau)$ into Taylor series we arrive at worm :

			$\displaystyle\sigma_{\beta\alpha}(\tau)=\delta_{\alpha\beta}\,+\int_{0}^{\tau}d\tau_{1}\>V_{\beta\alpha}\,e^{\tau_{1}E_{\beta\alpha}}$		(12)
		$\displaystyle+$	$\displaystyle\sum_{\gamma_{1}}\!\int_{0}^{\tau}\!\!\!d\tau_{2}\!\!\int_{0}^{\tau_{2}}\!\!\!d\tau_{1}V_{\beta\gamma_{1}}e^{\tau_{2}E_{\beta\gamma_{1}}}V_{\gamma_{1}\alpha}e^{\tau_{1}E_{\gamma_{1}\alpha}}+\ldots,\qquad$

where $E_{\beta\alpha}=E_{\beta}-E_{\alpha}$. The MC scheme is based on the statistical interpretation of the r.h.s. of (12) viewed as an average over an ensemble of graphs representing strings of the hopping transitions, or “kinks”, $V_{\gamma_{i+1}\gamma_{i}}$. A string is characterized by the number and types of kinks, as well as by their space-time positions. Graphs are sampled according to their non-negative weights given by the values of the corresponding integrands in Eq. (12).

Sampling protocol and estimators. Our MC scheme includes $\tau$ into the configuration space of graphs and allows us, on the one hand, to sample the $\tau$-dependence of ${\cal G}_{ba}(\tau)$ and $\overline{{\bf R}^{2}}(\tau)$ and then use Eqs. (9) and (13) to estimate the ground-state energy and the effective mass. On the other hand, this setup dramatically simplifies the minimalistic set of updates that deal exclusively with the last (in the $\tau$-domain) kink in the string. To be specific, we stochastically (i) add and remove the last bare-hopping kink using a pair of complementary updates worm , (ii) switch between the three types of the hopping terms, and (iii) sample the length of the last time interval separating the last kink from the state $\langle\,b\,|$. This scheme is ergodic and produces states $\langle\,b\,|$ that admit any allowed configuration of excited phonon modes. Clearly, one can add additional updates dealing with other than last kinks for better decorrelation of the sting configuration, see worm .

The proposed setup provides access to all the relevant system’s properties, including the same- and different-time correlation functions. They are measured when the value of $\tau$ is large enough and the asymptotic limit (7) is reached; its maximal value $\tau_{\rm max}$ is an important parameter controlling the accuracy of the projection onto the ground state. Monte Carlo estimators for observables are based on Eqs. (6) and (8) and their straightforward generalization for the different-time correlators. In accordance with the general theory (see, e.g., Refs. worm ) estimators do not involve additional computational costs (i.e. modifications of the configuration space and the sampling protocol) when the corresponding operators are either (i) diagonal in the basis ${\cal B}$ or (ii) expressed in terms of elementary $\hat{V}$-terms. The estimator for observable of type (i) is simply the eigenvalue of this observable for the state $|\gamma\rangle\in{\cal B}$ created by the string of kinks at a certain moment of imaginary time, $\tau_{*}$, close to $\tau/2$. The estimator for observable of type (ii) is minus the total number of the corresponding $\hat{V}$-kinks within a certain time interval divided by the duration of this interval. It can be arbitrarily long provided its boundaries are appropriately far from the string end points $\tau=0$ and $\tau\sim\tau_{\rm max}$.

Diagrammatics for the phonon sector. In the considered model, phonons do not interact with each other. This allows one to treat them diagrammatically, while keeping the path-integral representation for the particle sector only. The benefit of the diagrammatic representation for phonons (compared to their path-integral treatment described above) is the flexibility of dealing with dispersive modes. Here the phonon dispersion is accounted for by simply modifying the zero-temperature (to comply with the diagrammatic rule that all averages be specified in terms of the vacuum state) phonon propagators, $D_{s}({\bf r}_{1},\tau_{1};{\bf r}_{2},\tau_{2})\equiv D_{s}({\bf r}_{1}-{\bf r}_{2},\tau_{1},\tau_{2})$, with $s=1,\dots,d$ enumerating directions of lattice bonds. For sign-free formulation, $D_{s}$ should be non-negative. One can still access (i) the particle configurations, (ii) the distribution of phonon numbers, and the correlations between (i) and (ii). The path-integral treatment of dispersive phonons is also possible, and comes at a price of dealing with yet another family of kinks originating from phonon hopping. The benefit of this approach is detailed information about spatial configurations of excited phonon modes.

The mathematical structure of the perturbative expansion with diagrammatic treatment of phonons is similar to the series (12) provided the Greek subscripts are used exclusively for the particle states. Accordingly, $E_{\alpha}$ now refer to energies of the particle subsystem alone, while each kink in (12) representing the phonon-assisted hopping event is associated with one of the two end-points of the vacuum phonon propagator $D_{s}({\bf r}_{1},\tau_{1};{\bf r}_{2},\tau_{2})$. The other end belongs either to states $|\,a\,\rangle$ or $\langle\,b\,|$ or another phonon-assisted hopping event along the bond with the same direction $s$.

Illustrative results. For estimates of $T_{c}$, apart from the bi-polaron energy and effective mass, we determine $R_{*}^{2}$ from the probability distribution $P(R_{12})$ of finding particles at distance $R_{12}$ from their center of mass position

$$R_{*}^{2}\,=\,\sum_{R_{12}}R_{12}^{2}\,P(R_{12})\,.$$

(13)

Figure 1 shows the dependence of the ground-state energy and the effective mass on the strength of the on-site repulsion for coupling constant $g=t/\sqrt{2}$ and the adiabatic parameter $\gamma=1/16$. The ground-state energies are found to be smaller than the energies of two polarons for any realistic value of $U$, indicating that we are dealing with bound states (bipolarons). At $U=0$, the bipolaron effective mass is more than three times larger than its asymptotic limit $2m_{p}$, where $m_{p}$ is the mass of a single polaron. With increasing $U$, the effective mass drops significantly, and at $U\gtrsim 10$ it approaches the limiting value of $2m_{p}$.

In contrast to naive expectations that a significant drop in the effective mass should result in a substantial increase of $T_{c}$, the trend revealed in Fig. 2 is quite different. The estimate for the critical temperature remains relatively flat up to $U\sim 2$ and then starts to decrease, dropping by almost a factor of 3 at $U\sim 10$. This is because the decrease of $m_{*}$ is accompanied by a rapid increase of $R_{*}$, so that the product $m_{*}R_{*}^{2}$ does not show a pronounced maximum at any finite $U$.

A pronounced maximum of $T_{c}$ at a finite $U$ does exist at a larger value of the coupling constant, as illustrated by the data for bipolaron states at $g=t$ and $\gamma=1/16$ in Figs. 3 and 4. Here the behavior of both the ground-state energy and the effective mass is qualitatively similar to what we had at $g=t/\sqrt{2}$ (see Fig. 1). But the increase of $R_{*}$ with $U$ is not as dramatic and does not overcompensate the decrease of $m_{*}$ even at $U\sim 20$. The maximum of $T_{c}$ at $U>24$ with the maximal value of $T_{c}$ significantly larger than at $U=0$ does exist because for $U\to\infty$ the energy of the delocalized bipolaron state is extremely close to the bound state threshold.

Conclusions. Model (3)–(5) of bond (bi)polarons (in any spatial dimensions) allows a controlled and efficient numeric solution by Monte Carlo methods formulated in the path-integral representation for the electron(s) and either path-integral or real-space-diagrammatic representation for phonons. This is also true for the generalizations of (3)–(5) that include (i) dispersive bond phonons and (ii) additional density-displacement couplings as in the Holstein model. For illustrative purposes, we presented simulations of bipolaron states for model (3)–(5) in two dimensions.

In the context of the bipolaron mechanism of high-temperature superconductivity, both the effective mass and the size of bipolaron are equally important [see Eq. (2)]. Our results show that the positive effect of having a smaller effective mass can be readily overcompensated by the negative effect of having a larger bipolaron size; see Fig. 2.

There is a range of parameters where the net effect of the strong repulsion between the electrons is a substantial increase of the critical temperature for the superconducting transition; see Fig. 4. One may find this result rather counterintuitive given that normally the on-site repulsion suppresses Cooper pairing in the $s$-channel.

Acknowledgments. NP and BS acknowledge support by the National Science Foundation under Grant No. DMR-2032077.