Excitonic correlations in the equilibrium and voltage-biased bilayer Hubbard model: multi-orbital two-particle self-consistent approach

TL;DR

This paper develops a nonequilibrium, multi-orbital two-particle self-consistent (TPSC) framework and implements it on the real-frequency axis using Keldysh Green's functions to study a bilayer Hubbard model with interlayer interactions $V$ and hopping $W_\perp$. By self-consistently renormalizing spin and charge vertices through local two-particle sum rules and computing a spectral self-energy $\Sigma^{spc}$, the approach avoids spurious finite-temperature transitions in two dimensions and captures pseudogap phenomena driven by collective fluctuations. The authors show that interlayer bias $\Delta\mu$ can enhance excitonic fluctuations and bound-state formation, while large biases induce charge imbalance that suppresses these modes; they map phase tendencies across $U$, $V$, $\Delta\epsilon$, and $W_\perp$, including equilibrium and nonequilibrium steady states. The work provides a computationally efficient, gauge-friendly framework for investigating correlated multi-orbital systems under nonequilibrium conditions, with potential extensions to first-principles inputs and realistic materials.

Abstract

We develop a nonequilibrium multi-orbital extension of the two-particle self-consistent theory and apply it to the bilayer Hubbard model as a minimal platform to investigate correlation effects in the presence of interlayer interactions and tunneling. The method determines vertex corrections in the spin and charge channels self-consistently at the two-particle level, thereby avoiding the spurious finite-temperature phase transitions that limit dynamical mean-field theory in two dimensions. We derive the spectral self-energy and implement the framework directly on the real-frequency axis within the Keldysh nonequilibrium Green's function formalism, enabling the treatment of both equilibrium and non-equilibrium steady states without relying on numerical analytic continuation. As an application, we demonstrate that a pseudogap can emerge in the bilayer Hubbard model when spin, charge, or excitonic fluctuations become sufficiently strong. Instabilities in different channels are also evaluated in an unbiased manner across the parameter space. Remarkably, we find that the excitonic susceptibility grows with increasing interlayer bias, before it gets suppressed at large biases by the charge imbalance between the layers. This work establishes a versatile and computationally efficient framework for investigating correlated multi-orbital systems under nonequilibrium conditions.

Figures (5)

• Figure 1: (a) Illustration of the bi-layer structure coupled to two leads. (b) and (c) depict the energy diagram of the equilibrium ($\mu_A - \mu_B = 0$) and the nonequilibrium ($\mu_A - \mu_B = \Delta\mu$) setups. $\mu$ and $\epsilon$ refer to the chemical potential and the on-site energy of the lead and layer, respectively. Panels (d) and (e) show the RPA results of the local and $\mathbf{k}$-resolved spectral functions in the critical region close to the charge (d) and excitonic (e) ordering, respectively. The parameters are provided in the text.
• Figure 2: (a) Spin vertex $\Lambda^{sp}$ vs on-site Coulomb interaction $U$. (b) Imaginary part of the exciton susceptibility $\chi^{sp}_{ABBA}(\mathbf{q} = 0, \omega)$ with red and blue lines corresponding to the EQ and NEQ models, respectively. (c), (d), (e): Spin, charge and exciton instabilities $a$ as a function of crystal field splitting $\Delta \epsilon$. The black solid line shows the global minimum across the 1BZ, while dashed lines show the results for high symmetry points. Parameters: $U=4.0$, $V=2.0$ and $\beta = 4.0$.
• Figure 3: Phase diagram of the bilayer model in the $V$-$\Delta\epsilon$ plane for (a) $U=4.0$ and (b) $U=2.0$. The color map represents the minimum value of the $a$-coefficient across different sectors. Parameters: $W_\perp=0.0$ and $\beta=4.0$.
• Figure 4: Spectra of the exciton susceptibility, $-\Im \chi^{sp}_{\text{amp/pha}}(\omega,\mathbf{q})/\pi$, showing the amplitude and phase modes (Eq. \ref{['eq: susceptibility2']}) for different interlayer hoppings, $W_\perp = 0.0,,0.5,,1.0$. The top and bottom rows correspond to $V=1.0$ and $V=2.0$, respectively. The parameters used are $U=4$, $\Delta\epsilon=1.5$ and $\Delta\mu = 0$.
• Figure 5: (a),(b) Maximum value of $-\Im\chi^{sp,r}(\omega)/\pi$ across all frequencies for the amplitude and phase modes. The insets show the corresponding frequency positions of these maxima. (c) The exciton susceptibility for the phase mode $-\Im\chi^{sp,r}(\omega, \mathbf{q} = (\pi,\pi))/\pi$ for $\Delta \mu = 0.0$ (black), $\Delta \mu = 2.1$ (red dashed) and $\Delta \mu = 2.7$ (green dotted). Inset: densities vs $\Delta \mu$ for each layer. (d) Time evolution of the susceptibility $\chi^{sp,r}(t-t', \mathbf{q}=(\pi,\pi))$, which is real.