Diagrammatic Monte Carlo simulation of nonequilibrium systems

We generalize the recently developed diagrammatic Monte Carlo techniques for quantum impurity models from an imaginary time to a Keldysh formalism suitable for real-time and nonequilibrium calculations. Both weak-coupling and strong-coupling based methods are introduced, analyzed, and applied to the study of transport and relaxation dynamics in interacting quantum dots.

Figure 1

(Color online) Example of a Monte Carlo configuration corresponding to perturbation order $n=5$ and $n_{+}=3$, $n_{-}=2$.

Figure 2

Illustration of the band cutoffs considered in the simulations. The soft cutoff [(a) panel] is an exponentially decaying $\Gamma \left(\omega \right)$ which is centered at the chemical potential. The band with hard cutoff [(b) panel] is symmetric around $\omega =0$ and does not shift with the chemical potential.

Figure 3

(Color online) Segment representation of a Monte Carlo configuration corresponding to perturbation order four for spin up (upper contour) and two for spin down (lower contour). Dot creation operators are shown as full circles and annihilation operators as open circles. The segments represent the time intervals in which an electron of the corresponding spin resides on the dot.

Figure 4

(Color online) Segment configurations obtained from the expansion of the current [Eq. (62)] in powers of the dot-lead hybridization. There is a fixed operator $d_{\sigma }$ (red/gray open circle) at time $t$ and the hybridization functions connecting to this operator have only a left $\left(L\right)$ component.

Figure 5

Distribution of perturbation orders and average sign obtained using the hybridization expansion algorithm for a noninteracting dot with soft cutoff, $T=0$, $V=0$, ${ϵ}_d/\Gamma =-0.5$, and a single species of fermion. (a) The distribution of perturbation orders for different lengths of the contour ($t\Gamma =1.25,1.50,\dots ,3$ from left to right) and cutoff ${\omega }_c/\Gamma =40$. The average perturbation order grows $\sim t$. (b) The average sign as a function of time for indicated values of the cutoff.

Figure 6

Distribution of perturbation orders for different values of the interaction $U$, $T=0$, $V=0$, and ${ϵ}_d+U/2=0$. (a) The distribution of perturbation orders for the weak-coupling algorithm ($t\Gamma =1$, infinite bandwidth), where the average perturbation order grows $\sim U$. (b) The distribution of perturbation orders for the hybridization expansion algorithm ($t\Gamma =1.5$, ${\omega }_c/\Gamma =10$). Here there is almost no dependence on interaction strength.

Figure 7

(a) Perturbation order distribution for $t\Gamma =1$, $U/\Gamma =2$, $V/\Gamma =0$, and several positive values of $K$. The average signs in these simulations are 0.24 $\left(K=0.1\right)$, 0.04 $\left(K=1\right)$, and 0.0002 $\left(K=10\right)$. (b) Complex $K$ of norm 0.01. The best choice appears to be the negative $K$.

Figure 8

(Color online) Weak-coupling results for the double occupancy and density computed at $T=0$ and $V=0$ with parameter $K=-0.01$ and for infinite bandwidth. (a) The double occupancy for $U/\Gamma =2,3,4$ at half filling $\left({ϵ}_d+U/2=0\right)$. (b) The density per spin for ${ϵ}_d+U/2=\Gamma$ with the other parameters as the same. The full dots show the results from the real-time Hirsch-Fye method for $t\Gamma =1$ and indicated values of the number of time slices $L$.

Figure 9

Voltage dependence of the double occupancy for an infinite flat band, $T=0$, and half filling $\left({ϵ}_d+U/2=0\right)$. (a) The time evolution of the double occupancy for $U/\Gamma =2$ and indicated values of $V$. At finite voltage bias, the system reaches steady state much more rapidly than for $V=0$. (b) The steady-state value of the double occupancy as a function of voltage bias for $U/\Gamma =2$ and 3, $T=0$, and half filling.

Figure 10

Hybridization expansion results for the double occupancy and density. (a) The double occupancy for indicated values of $U$, band cutoff ${\omega }_c/\Gamma =10$, $\beta \Gamma =\nu \Gamma =10$, $V=0$, and ${ϵ}_d+U/2=0$. (b) The density per spin for the same parameters. The initial state is an empty dot decoupled from the leads and at $t=0$ the dot-lead hybridization is turned on.

Figure 11

(Color online) Current computed as a function of voltage bias in the infinite-bandwidth model using fourth-order perturbation theory in the interaction $U$ (from Ref. 26).

Figure 12

(Color online) Main panel: dashed-dotted line (red): current computed in mean-field theory for $U=12\Gamma$ as a function of voltage bias $V$ assuming negligible pseudothermal broadening. Dotted line (red): current computed in mean-field theory assuming a pseudothermal broadening equal to 20% of the voltage bias. Solid line (blue): current computed using mean-field theory with gap fixed at the $V=0$ value and negligible pseudothermal broadening; dashed line (blue): current computed using mean-field theory with gap fixed at the $V=0$ value and pseudothermal broadening equal to 20% of the voltage bias. Dashed double-dotted line (black): current computed for the noninteracting model (negligible pseudothermal broadening). Inset: dot magnetization as a function of voltage bias $V/\Gamma$ computed in mean-field theory for negligible pseudothermal broadening (dotted-dashed line, red on line) and moderate pseudothermal broadening (dotted line, red). All computations were performed for a hard cutoff with ${\omega }_c=10\Gamma$ and $\nu \Gamma =10$.

Figure 13

Weak-coupling expansion results for the current at temperature $T=0$. (a) The time evolution for $U/\Gamma =2$ and indicated values of the voltage bias. The value at time $t=0$ corresponds to the noninteracting current $I_0=-4\text{Im}A\left(0,0\right)$. As the interaction is turned on, the magnitude of the current decreases and eventually converges to the steady-state value for the interacting dot. (b) The time dependence of the current for $V/\Gamma =4$, and $U/\Gamma =2$ and 3.

Figure 14

Weak-coupling expansion results for the voltage dependence of the steady-state current for $T=0$ and indicated values of the interaction $U$. (a) The total current and (b) the reduction in the steady-state current produced by the interaction, $I\left(U\right)-I\left(0\right)$. The interaction correction peaks at a value of $V$ which is comparable to the interaction, or somewhat larger than twice the “level broadening” $\Gamma$. The thick black line in (b) shows the prediction for $U/\Gamma =2$ from fourth-order perturbation theory (Ref. 26).

Figure 15

Hybridization expansion results for the current through an interacting dot with $U/\Gamma =8$, ${ϵ}_d+U/2=0$, and a hard band cutoff ${\omega }_c/\Gamma =10$ $\left(\beta \Gamma =\nu \Gamma =10\right)$. (a) The left and right currents for $V/\Gamma =0$ and 5. $I_L$ is the current from the dot to the left lead and $I_R$ is the current from the right lead to the dot. (b) The average current $I=\left(I_L+I_R\right)/2$ and $dn/dt=I_L-I_R$ for the same parameters.

Figure 16

Hybridization expansion results for the current through an interacting dot with ${ϵ}_d+U/2=0$, and a hard bandcutoff ${\omega }_c/\Gamma =10$ $\left(\beta \Gamma =\nu \Gamma =10\right)$. The steady-state dot occupancy for this ${ϵ}_d$ is one (half filling). (a) The average current $I=\left(I_L+I_R\right)/2$ for $U/\Gamma =8$ and indicated values of $V$. (b) The current at time $t\Gamma =1$ (solid lines) and $t\Gamma =1.25$ (dashed lines) as a function of voltage for indicated values of the interaction. The current for the noninteracting dot has been measured at $t\Gamma =2$.

Figure 17

(a) Comparison to the mean-field result. The circles show the (exact) noninteracting current for ${\omega }_c/\Gamma =10$, $\beta \Gamma =\nu \Gamma =10$, which is in good agreement with the hybridization expansion result measured at $t\Gamma =2$, especially for $V/\Gamma \gtrsim 3$. The triangles show the current obtained with the “fixed gap” calculation for $U/\Gamma =8$ and 12, and an effective temperature $T=0.05V$ (open symbols) and $T=0.2V$ (full symbols). Stars show the Monte Carlo results for $U/\Gamma =8$ and 12 measured at $t\Gamma =1.25$. (b) Comparison of the Monte Carlo results measured at $t\Gamma =1$ (circles), 1.25 (stars), and 2 (diamonds) to the current (for infinite bandwidth) deduced from the fourth-order perturbation calculation of Ref. 26.