Quantum Monte Carlo for correlated out-of-equilibrium nanoelectronic devices

We present a simple, general purpose, quantum Monte Carlo algorithm for out-of-equilibrium interacting nanoelectronic systems. It allows one to systematically compute the expansion of any physical observable (such as current or density) in powers of the electron-electron interaction coupling constant $U$. It is based on the out-of-equilibrium Keldysh Green's function formalism in real-time and corresponds to evaluating all the Feynman diagrams to a given order $U^n$ (up to $n=15$ in the present work). A key idea is to explicitly sum over the Keldysh indices in order to enforce the unitarity of the time evolution. The coefficients of the expansion can easily be obtained for long-time, stationary regimes, even at zero temperature. We then illustrate our approach with an application to the Anderson model, an archetype interacting mesoscopic system. We recover various results of the literature such as the spin susceptibility or the “Kondo ridge” in the current-voltage characteristics. In this case, we found the Monte Carlo free of the sign problem even at zero temperature, in the stationary regime and in absence of a particle-hole symmetry. The main limitation of the method is the lack of convergence of the expansion in $U$ for large $U$, i.e., a mathematical property of the model rather than a limitation of the Monte Carlo algorithm. Standard extrapolation methods of divergent series can be used to evaluate the series in the strong correlation regime.

Figure 1

Sketch of a typical mesoscopic system: a central interacting region (red) is connected to several (semi-infinite) noninteracting electrodes (blue) with finite temperatures $T_i$ and chemical potentials ${\mu }_i$.

Figure 2

$Q_n$ as a function of ${\epsilon }_d$ for $n=0$, 1, 2, and 3. The calculations are performed using a direct evaluation of the integrals in Eq. (14) using the Simpson rule for $t=20, \gamma =0.5$, and $T=0$.

Figure 3

Color plot of the integrand of $Q_2$ as a function of the two times $u_1$ and $u_2$ for model A with ${\mu }_L={\mu }_R=0, {\epsilon }_d=0, T=0$, and $t=10$. The four panels correspond to the four possible values of the two Keldysh indices $a_1$ and $a_2$. The explicit form of the integrand is $f(u_1,u_2,a_1,a_2)=-\Im m{(-1)}^{{\sum }_i{a}_i}det{\mathbf{M}}_2(u_1,u_2,a_1,a_2)$.

Figure 4

Same parameters as in Fig. 3 but the integrand has now been summed over Keldysh indices. The color plot represents $f(u_1,u_2)=i{\sum }_{a_1,a_2}{(-1)}^{{\sum }_i{a}_i}det{\mathbf{M}}_2(u_1,u_2,a_1,a_2)$ ($f$ is real). Note that the integrand is now real, positive, and concentrated around $u_1=u_2=t$.

Figure 5

Same parameters as in Fig. 3 but the integrand now uses the matrix ${\mathbf{P}}_2$ instead of ${\mathbf{M}}_2$, i.e., is associated to the partition function instead of an observable. The color plot represents $f(u_1,u_2,a_1,a_2)=-\Im m{(-1)}^{{\sum }_i{a}_i}det{\mathbf{P}}_2(u_1,u_2,a_1,a_2)$.

Figure 6

QMC results for model A at $\gamma =1/2, T=0$, and $t=10$ for ${\epsilon }_d=0$ (green squares) and ${\epsilon }_d=-0.5$ (blue circles), as a function of the order $n$. Red diamonds: model B with $\gamma =1/2, T=0$, and $t=10$ for ${\epsilon }_d=-0.5$ and $\alpha =0.5$. Top panel: $Q_n$, central panel: $c_n$ and bottom panel: average sign $s_n$.

Figure 7

QMC results for model A at $\gamma =1/2, T=0$, and $t=10$. (Top) Charge $Q(N,U)$ as a function of $U$, for different $N$ and for ${\epsilon }_d=0$: $N=5$ (black), 7 (red), 13 (blue), and 15 (green). (Bottom) $Q(N,U)$ as a function of $1/N$ for different $U$: $U=0.2$ (orange), $0.4$ (red), $0.6$ (black), and $0.8$ (blue).

Figure 8

Hartree-Fock series of model A as a function of $U$: exact curve $Q^{\mathrm{HF}}\left(U\right)$ (thick line) and partial sums $Q^{\mathrm{HF}}(N,U)$ for $N=10,20,30, {\epsilon }_d=0$, and $\gamma =1/2$. The series has a divergence at $U=1$. The inset shows the analytical structure of $Q^{\mathrm{HF}}\left(U\right)$ in the complex plane $(\mathrm{Re}U,\mathrm{Im}U)$: branch cut of $Q^{\mathrm{HF}}\left(U\right)$ (thick magenta line) and convergence radius of the Hartree-Fock series (dashed line).

Figure 9

Role of the noninteracting Hamiltonian around which the expansion is done; $Q\left(U\right)$ at ${\epsilon }_d=0$ with an extra potential $\delta {\epsilon }_d$ in the noninteracting Hamiltonian and a corresponding compensating term $\alpha =\delta {\epsilon }_d/\overline{U}$ in the perturbation. The green line corresponds to the original series with $\delta {\epsilon }_d=0$ in its regime of convergence. All curves for $\gamma =1/2, T=0$, and $t=10$.

Figure 10

(Top) Absolute value of the charges $Q_n$ for different values of $\alpha ,\delta {\epsilon }_d$, at ${\epsilon }_d=0, \gamma =1/2, T=0$, and $t=10$: $\alpha =\delta {\epsilon }_d=0$ (black), $\alpha =\delta {\epsilon }_d=0.5$ (blue), $\alpha =\delta {\epsilon }_d=0.6$ (green), $\alpha =\delta {\epsilon }_d=0.7$ (red). (Bottom) Corresponding dependence of the charge $Q(N,U=1)$ as a function of $1/N$.

Figure 11

Charge $Q\left(U\right)$ for model A with ${\epsilon }_d=0, \gamma =1/2, T=0$, and $t=10$. Red symbols correspond the use of $\alpha$ parameters (different symbols for the same value of $U$ correspond to different $\delta \epsilon$ or $\alpha$). Thin lines correspond to partial sum with the bare technique ($\alpha =0$): $Q(N=6,U)$ (cyan), $Q(N=7,U)$ (green), and $Q(N=8,U)$ (blue). Dashed line: benchmark calculation performed with hybridization QMC in imaginary time. Blue squares: extrapolation using the homographic technique, see Sec. 7d.

Figure 12

Homographic transformation of the model A series for $Q\left(U\right)$ (same data as Fig. 7) with $b=6$. (Top) Effective radius $R_W\left(a\right)$ (ful black line) and homographic transformation $W(a,U=6)$ (dashed red line) as a function of position of the singularity $a$. (Bottom) resummation result $Q(U=6)$ as a function of $a$. The dashed line indicates the exact result $Q\approx 0.22$. (Inset) ${\overline{Q}}_n$ as a function of $n$ for $a=-0.25$ (squares) and the initial series (circles).

Figure 13

Current-voltage characteristics $I\left(V_b\right)$ for model B with $\alpha =0.5, \gamma =0.5$, and $T=0$. The current is in units of energy time $e/h$. Black line and black circles: particle-hole symmetric point ${\epsilon }_d=0$ for $U=0$ and $U=1.2$, respectively. Red dot-dashed line and red squares ${\epsilon }_d=-0.5$ for $U=0$ and $U=1.2$. Dashed line: perfect transmission $I=V_b$. (Inset) Convergence of the results as a function of $1/N$ where $N$ is the total number of moments included for $V_b=0.3,{\epsilon }_d=-0.5,U=1.2$ (red squares) and $V_b=0.2,{\epsilon }_d=0,U=1.2$ (black circles).

Figure 14

First six moments of the currents $I_n{U}^n$ vs time $t$. The moments have been rescaled with $U=3$ for visibility. Model B with $\alpha =0.5, {\epsilon }_d=0, \gamma =0.5, T=0$, and $V_b=0.5$.

Figure 15

First six moments of the sign $s_n$ associated with the currents $I_n$ vs time $t$. Model B with $\alpha =0.5, {\epsilon }_d=0, \gamma =0.5, T=0$, and $V_b=0.5$. Assymptotic $s_6=0.025$.

Figure 16

First six moments of the weight $c_n$ associated with the currents $I_n$ vs time $t$. Model B with $\alpha =0.5, {\epsilon }_d=0, \gamma =0.5, T=0$, and $V_b=0.5$.

Figure 17

Zero temperature magnetization. $M_z=(Q_{\uparrow }-Q_{\downarrow })/2$ as a function of magnetic field $h$ for model A with ${\epsilon }_d=-0.25, \gamma =1/4, T=0$ and $t=50$. Lines are analytical results (from Bethe ansatz Ref. [37]), in particular red line is at $U=0$, black line at $U=0.7\mathrm{\Gamma }$ ($\mathrm{\Gamma }=4{\gamma }^2=1/4$), blue line at $U=1.5\mathrm{\Gamma }$, dashed lines with points are results from QMC. The actual calculations were made with model B with $\alpha =0.5$ ($h\ge 0.025$) and $\alpha =0.8$ ($h<0.025$), while enforcing ${\epsilon }_d=\alpha U-0.25$.

Figure 18

Extrapolation of the series. $Q(U=1,\epsilon )$ as a function of $\epsilon$ for model A with ${\epsilon }_d=0, T=0$, and $\gamma =0.5$ for various maximum order $N=4,5,6,7$. The circles show a simple extrapolation from a linear regression in the region where the various curves overlap. We do recover $Q(U=1)\approx 0.32$ well beyond the apparent convergence radius of the series.

Figure 19

(Top) energy of the (unique) bound-state as a function of ${\epsilon }_d$. Bottom panel: absolute value of the diagonal part of the lesser Green's function $g_{00}^{<}\left(t\right)$ as a function of $t$, for ${\xi }_d=-0.5$ and $T=0$. For $\gamma =0.5$, there is no bound state.

Figure 20

$Q_2$ as a function of $1/t$, for $T=0$. Top panel is for $\gamma =1$ (bound state), bottom panel for $\gamma =0.5$ (no bound state). In each panel, the continuous (blue) line is for ${\epsilon }_d=0$, the dashed (red) line for ${\epsilon }_d=-0.5$

Figure 21

$Q_2$ as a function of ${\epsilon }_d$, for $T=0$ and different values of $t$. Top panel is for $\gamma =0.5$ (no bound state for $\left|\epsilon \right|<1.5$), bottom panel for $\gamma =1$ (bound state whatever ${\epsilon }_d$). dash-dotted line for $t=40$.