Real-time dynamics in quantum impurity models with diagrammatic Monte Carlo

We extend the recently developed real-time diagrammatic Monte Carlo method, in its hybridization expansion formulation, to the full Kadanoff-Baym-Keldysh contour. This allows us to study real-time dynamics in correlated impurity models starting from an arbitrary, even interacting, initial density matrix. As a proof of concept, we apply the algorithm to study the nonequilibrium dynamics after a local quantum quench in the Anderson impurity model. Being a completely general approach to real-time dynamics in quantum impurity models, it can be used as a solver for nonequilibrium dynamical mean field theory.

Figure 1

Kadanoff-Baym-Keldysh contour, which starts at time 0 on the first branch, runs up to time $t=t_{\ast }$ then back from $t_{\ast }$ to $t=0$ and finally along the imaginary axis from $t=0$ to $t=-i\beta$.

Figure 2

(Color online) An example of configuration $\mathcal{C}$ sampled by diagMC algorithm.

Figure 3

(Color online) Probability distribution of different perturbative orders $k$ sampled during the simulation. Data refers to a local quench in the Anderson Impurity Model, starting from $U_{-}=0$ to $U_{+}=10\Gamma$ at particle-hole symmetry, for $T=0.1\Gamma$ and $W=10\Gamma$. Left panel shows the statistics of kinks along the whole contour, while the top and bottom right panels display the histograms for the Matsubara and Keldysh sector, respectively.

Figure 4

(Color online) Scaling of the average number of kinks with maximum time $t_{\star }$. Left panel shows data at fixed $U_{+}=0$, tuning the strength of the Coulomb repulsion $U_{-}=0,5,10$. See the perfect linear scaling with the same slope $\alpha =d\overline{k}/dt$. Right panel displays data at fixed $U_{-}=10$, $U_{+}=0$ tuning the strength of the conduction bandwidth $W$. We see how the slope $\alpha$ increases with bandwidth making increasingly difficult to access large time scales in the regime $W⪢\Gamma$.

Figure 5

(Color online) Average sign as a function of time $t$ for different initial preparations. We clearly see an exponential decay on a very short-time scale. Left panel shows data obtained fixing the final value of the interaction $U_{+}=0$ and tuning the initial value $U=0,5,10$. We see a slight increase in the average sign. Right panel shows the dependence of $\overline{\eta }$ from the bandwidth of conduction electrons and suggest that much longer time scales can be reached in the regime $W\sim \Gamma$.

Figure 6

(Color online) Quench dynamics in a Resonant Level Model after a sudden change of the energy level from ${\epsilon }_{d-}=0$ to ${\epsilon }_{d+}\ne 0$. Dashed lines is the exact solution for $n\left(t\right)$ as obtained by a standard methods. Points are diagMC results obtained at $T=0.1\Gamma$. We also add the dynamics for the trivial case ${\epsilon }_{d+}={\epsilon }_{d-}=0$ to show that unitarity is actually preserved by diagMC.

Figure 7

(Color online) Nonequilibrium dynamics of double occupancy $D\left(t\right)$ in the Anderson impurity model after a local quantum quench of the interaction strength at $T=0.1\Gamma$ and particle-hole symmetry. In the upper panel we start from an initial state with $U_{-}=10\Gamma$ and a very low double occupation and quench to different values of $U_{+}/\Gamma =0,2.5,5$ from top to bottom. In the lower panel, the opposite protocol is considered, namely, we start from different values of $U_{-}/\Gamma =2,4,6$ from top to bottom and quench to the same final $U_{+}=0$. In both cases we see that after a rather short transient the system relaxes to a new equilibrium state.

Figure 8

(Color online) Nonequilibrium dynamics of the impurity density after a quench of the energy level from ${\epsilon }_{d-}=0$ to ${\epsilon }_{d+}=2.5\Gamma$. We compare the dynamics for different values of the Coulomb repulsion $U=0,2.5,5,7.5,10$ from bottom to top, revealing a much faster thermalization in the correlated case due to Coulomb blockade. Inset: thermal value of the impurity density as a function of the level position for $U=0,2.5,5,7.5,10$ from top to bottom. For $U⪢\Gamma$, the curve is almost flat around ${\epsilon }_d=0$ hence departure from equilibrium is suppressed after a short transient.

Figure 9

(Color online) Spin dynamics in the Anderson model after a sudden quench of a local magnetic field. We start from a partially polarized impurity, with $h_{-}=5.0$, $T=0.1\Gamma$ and different values of local interaction $U$ (top panel). At time $t>0$, the magnetic field is switched off and the magnetization is allowed to relax toward an unpolarized steady state. We see that upon increasing the interaction $U$, the dynamics slows down. Due to fine-time resolution we cannot follow the decay of the spin toward zero magnetization. However, a large magnetic field in the final state gives rise to a rather fast relaxation.

Figure 10

(Color online) Nonequilibrium dynamics for an Anderson Impurity coupled to a gapped fermionic reservoir. We plot the real-time dynamics for double occupancy at the impurity site after a quench of the interaction from $U_{-}=10\Gamma$ to $U_{+}=0$, at particle-hole symmetry ${\epsilon }_{d+}={\epsilon }_{d-}=0$ and for $T=0.1\Gamma$. Different values of the gap in the bath $E_g=0,1,2.5,5.0$ are considered, resulting into very different dynamics. Contrarily to the gapless case ($E_g=0$, black points) which quickly approaches the thermal plateau fixed by PH symmetry and indicated by an arrow, $D_{\text{therm}}=1/4$, we see that due to the finite gap in the spectrum the real-time dynamics slows down thus preventing us to conclude on the long time behavior of $D\left(t\right)$. However, for very large values of $E_g$ (see $E_g=5\Gamma$) the dynamics seems actually to reach a steady state where the double occupation is different from $D_{\text{therm}}$.

Figure 11

(Color online) Nonequilibrium dynamics for an Anderson Impurity coupled to a gapped fermionic reservoir. We set the gap in the spectrum equal to $E_g=1.0$ and plot the real-time dynamics for the double occupancy after a quench of the interaction, at $T=0.1\Gamma$ and PH symmetry. Two kind of processes are considered, namely, a quench from $U_{-}=0$ to $U_{+}\ne 0$ and the reverse process from $U_{-}\ne 0$ to $U_{+}=0$, which differ among each other for the average work done during the quench, see Eq. (60) in the main text. Top panel shows data for $U_{-}=10\Gamma$, $U_{+}=0$ and vice versa, while bottom panel data for $U_{-}=2.5\Gamma$, $U_{+}=0$ and vice versa. We see that, provided the average work done is sufficiently larger than the gap $2E_g$, a fast thermalization can occur also in a gapped model (see top panel, black points). As opposite, when the amount of energy pumped into the system is too small the dynamics slows down and we cannot conclude with present data whether thermalization takes place or not.

Figure 12

(Color online) Nonequilibrium dynamics for an Anderson Impurity coupled to a pseudogapped fermionic reservoir. We consider the model at PH symmetry and $T=0.1\Gamma$. We fix the quench parameters, namely, the initial and final value of the interaction, equal, respectively, to $U_{-}=10$ and $U_{+}=0$, and tune the pseudogap exponent from $r=0$ (gapless metallic state) to $r=4$. The depletion of low energy states in the DoS reflects into a much slower dynamics, which eventually, for large enough $r$, seems to prevent the system from reaching the value of $D_{\text{therm}}=1/4$, which is fixed by PH symmetry and by the choice of $U_{+}=0$. However, due to finite time resolution we cannot conclude with present data whether thermalization occurs or not on a longer time scale.

Figure 13

(Color online) Nonequilibrium dynamics for an Anderson Impurity coupled to a pseudogapped fermionic reservoir. We consider the PH symmetric point, at $T=0.1\Gamma$ and for $r=1$. We study how the dynamics changes while tuning the average work $\overline{W}$ done during the quench (see text). We fix the initial interaction $U_{-}=0$ and perform a quench to $U_{+}=10\Gamma$ (top panel). We compare the results with the reverse process, from $U_{-}\ne 0$ to $U_{+}=0$. Despite the power-law DoS, the dynamics in the case of a large quantum quench (large average work) turns to be quite fast. On the contrary, in the low work regime (bottom panel) the dynamics is much slower and we cannot see whether thermalization take place or not at longer time scales.