Exact real-time dynamics of single-impurity Anderson model from a   single-spin hybridization-expansion

Patryk Kubiczek; Alexey N. Rubtsov; Alexander I. Lichtenstein

arXiv:1904.12582·cond-mat.str-el·August 7, 2019

Exact real-time dynamics of single-impurity Anderson model from a single-spin hybridization-expansion

Patryk Kubiczek, Alexey N. Rubtsov, Alexander I. Lichtenstein

PDF

1 Repo

TL;DR

This paper introduces a modified quantum Monte Carlo solver for the time-dependent single-impurity Anderson model, reducing computational complexity and enabling longer simulation times.

Contribution

The authors develop a semi-analytical hybridization-expansion quantum Monte Carlo method that improves efficiency by focusing on one spin channel, extending simulation capabilities.

Findings

01

Enables simulation of longer time scales in impurity models.

02

Reduces the dynamical sign problem in quantum Monte Carlo simulations.

03

Successfully applied to calculate electric current in impurity systems.

Abstract

In this work we introduce a modified real-time continuous-time hybridization-expansion quantum Monte Carlo solver for a time-dependent single-orbital Anderson impurity model: CT-1/2-HYB-QMC. In the proposed method the diagrammatic expansion is performed only for one out of the two spin channels, while the resulting effective single-particle problem for the other spin is solved semi-analytically for each expansion diagram. CT-1/2-HYB-QMC alleviates the dynamical sign problem by reducing the order of sampled diagrams and makes it possible to reach twice as long time scales in comparison to the standard CT-HYB method. We illustrate the new solver by calculating an electric current through impurity in paramagnetic and spin-polarized cases.

Equations84

H (t)

H (t)

H_{loc} (t)

H_{bath, σ} (t)

H_{hyb, σ} (t)

Z

Z

= k = 0 \sum \infty (- 1)^{k} \int_{0_{+}}^{- i β} d t_{1}^{c} \dots \int_{t_{k - 1}^{c}}^{- i β} d t_{k}^{c} \int_{0_{+}}^{- i β} d t_{1}^{a} \dots \int_{t_{k - 1}^{a}}^{- i β} d t_{k}^{a}

w_{loc + bath ↑} ({t_{m}^{c}}, {t_{n}^{a}}) w_{bath ↓} ({t_{m}^{c}}, {t_{n}^{a}}) .

w_{bath ↓} ({t_{m}^{c}}, {t_{n}^{a}}) = T r_{c_{↓}} [T_{C} e^{- i \int_{C} d t H_{bath ↓} (t)} C_{↓} (t_{k}^{c}) \dots C_{↓} (t_{1}^{c}) C_{↓}^{†} (t_{k}^{a}) \dots C_{↓}^{†} (t_{1}^{a})],

w_{bath ↓} ({t_{m}^{c}}, {t_{n}^{a}}) = T r_{c_{↓}} [T_{C} e^{- i \int_{C} d t H_{bath ↓} (t)} C_{↓} (t_{k}^{c}) \dots C_{↓} (t_{1}^{c}) C_{↓}^{†} (t_{k}^{a}) \dots C_{↓}^{†} (t_{1}^{a})],

w_{loc + bath ↑} ({t_{m}^{c}}, {t_{n}^{a}})

w_{loc + bath ↑} ({t_{m}^{c}}, {t_{n}^{a}})

= T r_{d_{↑}, d_{↓}, c_{↑}} [T_{C} e^{- i \int_{C} d t [H_{loc} (t) + H_{bath ↑} (t) + H_{hyb ↑} (t)]} d_{↓}^{†} (t_{k}^{c}) \dots d_{↓}^{†} (t_{1}^{c}) d_{↓} (t_{k}^{a}) \dots d_{↓} (t_{1}^{a})] .

Δ_{σ} (t, t^{'}) = α p \sum V_{p σ}^{*} (t) g_{α p σ} (t, t^{'}) V_{p σ} (t^{'}),

Δ_{σ} (t, t^{'}) = α p \sum V_{p σ}^{*} (t) g_{α p σ} (t, t^{'}) V_{p σ} (t^{'}),

g_{α p σ} (t, t^{'}) = ⟨ T_{C} c_{α p σ} (t) c_{α p σ}^{†} (t^{'})⟩ = - i [θ_{C} (t, t^{'}) - \frac{1}{e ^{β ϵ_{p σ} (0)} + 1}] e^{- i \int_{t^{'}}^{t} d \overline{t} ϵ_{p σ} (\overline{t})},

g_{α p σ} (t, t^{'}) = ⟨ T_{C} c_{α p σ} (t) c_{α p σ}^{†} (t^{'})⟩ = - i [θ_{C} (t, t^{'}) - \frac{1}{e ^{β ϵ_{p σ} (0)} + 1}] e^{- i \int_{t^{'}}^{t} d \overline{t} ϵ_{p σ} (\overline{t})},

w_{bath ↓} ({t_{m}^{c}}, {t_{n}^{a}}) = i^{k} Z_{bath ↓} det [Δ_{↓} (t_{m}^{c}, t_{n}^{a})]_{m, n = 1, \dots, k},

w_{bath ↓} ({t_{m}^{c}}, {t_{n}^{a}}) = i^{k} Z_{bath ↓} det [Δ_{↓} (t_{m}^{c}, t_{n}^{a})]_{m, n = 1, \dots, k},

w_{loc + bath ↑} ({t_{m}^{c}}, {t_{n}^{a}}) = φ ({t_{m}^{c}}, {t_{n}^{a}}) \cdot det [\mathbbm 1 + U (- i β, 0_{+}; {t_{m}^{c}}, {t_{n}^{a}})] .

w_{loc + bath ↑} ({t_{m}^{c}}, {t_{n}^{a}}) = φ ({t_{m}^{c}}, {t_{n}^{a}}) \cdot det [\mathbbm 1 + U (- i β, 0_{+}; {t_{m}^{c}}, {t_{n}^{a}})] .

φ ({t_{m}^{c}}, {t_{n}^{a}}) = {1 (- 1)^{k} if t_{1}^{c} > t_{1}^{a} else} \cdot e^{- i \int_{C} d t E_{d ↓} (t) n_{↓} (t)} .

φ ({t_{m}^{c}}, {t_{n}^{a}}) = {1 (- 1)^{k} if t_{1}^{c} > t_{1}^{a} else} \cdot e^{- i \int_{C} d t E_{d ↓} (t) n_{↓} (t)} .

h_{↑} (t) = E_{d ↑} (t) + U (t) n_{↓} (t) V_{1 ↑} (t) ⋮ V_{N ↑} (t) V_{1 ↑}^{*} (t) ε_{1 ↑} (t) ⋮ 0 \dots \dots ⋱ 0 V_{N ↑}^{*} (t) 00 ε_{N ↑} (t) .

h_{↑} (t) = E_{d ↑} (t) + U (t) n_{↓} (t) V_{1 ↑} (t) ⋮ V_{N ↑} (t) V_{1 ↑}^{*} (t) ε_{1 ↑} (t) ⋮ 0 \dots \dots ⋱ 0 V_{N ↑}^{*} (t) 00 ε_{N ↑} (t) .

U (t, t^{'}; {t_{m}^{c}}, {t_{n}^{a}}) = T_{C} e^{- i \int_{t}^{t^{'}} d \overset{ˉ}{t} h_{↑} (\overset{ˉ}{t})}

U (t, t^{'}; {t_{m}^{c}}, {t_{n}^{a}}) = T_{C} e^{- i \int_{t}^{t^{'}} d \overset{ˉ}{t} h_{↑} (\overset{ˉ}{t})}

Tr_{f} [T_{C} e^{- i \int_{C} d t \sum_{ab} f_{a}^{†} h_{ab} (t) f_{b}}] = det [\mathbbm 1 + T_{C} e^{- i \int_{C} d t h (t)}] .

Tr_{f} [T_{C} e^{- i \int_{C} d t \sum_{ab} f_{a}^{†} h_{ab} (t) f_{b}}] = det [\mathbbm 1 + T_{C} e^{- i \int_{C} d t h (t)}] .

φ ({t_{m}^{c}}, {t_{n}^{a}}) = {1 (- 1)^{k} if t_{1}^{c} > t_{1}^{a} else} \cdot segments \prod e^{- i \int_{seg.} d t E_{d ↓} (t)}

φ ({t_{m}^{c}}, {t_{n}^{a}}) = {1 (- 1)^{k} if t_{1}^{c} > t_{1}^{a} else} \cdot segments \prod e^{- i \int_{seg.} d t E_{d ↓} (t)}

U (- i β, 0_{+}; {t_{m}^{c}}, {t_{n}^{a}}) = {u_{1} (- i β, t_{k}^{c}) u_{0} (t_{k}^{c}, t_{k}^{a}) u_{1} (t_{k}^{a}, t_{k - 1}^{c}) \dots u_{0} (t_{1}^{c}, t_{1}^{a}) u_{1} (t_{1}^{a}, 0_{+}) u_{0} (- i β, t_{k}^{a}) u_{1} (t_{k}^{c}, t_{k}^{a}) u_{0} (t_{k}^{c}, t_{k - 1}^{a}) \dots u_{1} (t_{1}^{a}, t_{1}^{c}) u_{0} (t_{1}^{c}, 0_{+}) if t_{1}^{c} > t_{1}^{a} else,

U (- i β, 0_{+}; {t_{m}^{c}}, {t_{n}^{a}}) = {u_{1} (- i β, t_{k}^{c}) u_{0} (t_{k}^{c}, t_{k}^{a}) u_{1} (t_{k}^{a}, t_{k - 1}^{c}) \dots u_{0} (t_{1}^{c}, t_{1}^{a}) u_{1} (t_{1}^{a}, 0_{+}) u_{0} (- i β, t_{k}^{a}) u_{1} (t_{k}^{c}, t_{k}^{a}) u_{0} (t_{k}^{c}, t_{k - 1}^{a}) \dots u_{1} (t_{1}^{a}, t_{1}^{c}) u_{0} (t_{1}^{c}, 0_{+}) if t_{1}^{c} > t_{1}^{a} else,

u_{n_{↓}} (t, t^{'}) = T_{C} exp - i \int_{t^{'}}^{t} d \overset{ˉ}{t} E_{d ↑} (\overset{ˉ}{t}) + U (\overset{ˉ}{t}) n_{↓} V_{1 ↑} (\overset{ˉ}{t}) ⋮ V_{N ↑} (\overset{ˉ}{t}) V_{1 ↑}^{*} (\overset{ˉ}{t}) ε_{1 ↑} (\overset{ˉ}{t}) ⋮ 0 \dots \dots ⋱ 0 V_{N ↑}^{*} (\overset{ˉ}{t}) 00 ε_{N ↑} (\overset{ˉ}{t}) . .

u_{n_{↓}} (t, t^{'}) = T_{C} exp - i \int_{t^{'}}^{t} d \overset{ˉ}{t} E_{d ↑} (\overset{ˉ}{t}) + U (\overset{ˉ}{t}) n_{↓} V_{1 ↑} (\overset{ˉ}{t}) ⋮ V_{N ↑} (\overset{ˉ}{t}) V_{1 ↑}^{*} (\overset{ˉ}{t}) ε_{1 ↑} (\overset{ˉ}{t}) ⋮ 0 \dots \dots ⋱ 0 V_{N ↑}^{*} (\overset{ˉ}{t}) 00 ε_{N ↑} (\overset{ˉ}{t}) . .

w_{loc + bath ↑} ({t_{m}^{c}}, {t_{n}^{a}}) = φ ({t_{m}^{c}}, {t_{n}^{a}}) \cdot Z_{bath ↑} \cdot Z_{d ↑} ({t_{m}^{c}}, {t_{n}^{a}}),

w_{loc + bath ↑} ({t_{m}^{c}}, {t_{n}^{a}}) = φ ({t_{m}^{c}}, {t_{n}^{a}}) \cdot Z_{bath ↑} \cdot Z_{d ↑} ({t_{m}^{c}}, {t_{n}^{a}}),

Z_{d ↑} ({t_{m}^{c}}, {t_{n}^{a}}) = \int D (d_{↑}, d_{↑}^{*}) e^{i \int_{C} \int_{C} d t d t^{'} d_{↑}^{*} (t) G_{↑}^{- 1} (t, t^{'}) d_{↑} (t^{'})}

Z_{d ↑} ({t_{m}^{c}}, {t_{n}^{a}}) = \int D (d_{↑}, d_{↑}^{*}) e^{i \int_{C} \int_{C} d t d t^{'} d_{↑}^{*} (t) G_{↑}^{- 1} (t, t^{'}) d_{↑} (t^{'})}

G_{↑}^{- 1} (t, t^{'})

G_{↑}^{- 1} (t, t^{'})

G_{0 ↑}^{- 1} (t, t^{'})

Z_{d ↑} ({t_{m}^{c}}, {t_{n}^{a}}) = det [i G_{↑}^{- 1}] .

Z_{d ↑} ({t_{m}^{c}}, {t_{n}^{a}}) = det [i G_{↑}^{- 1}] .

Z_{d ↑} ({t_{m}^{c}}, {t_{n}^{a}}) = z_{d ↑} ({t_{m}^{c}}, {t_{n}^{a}}) \cdot det [\mathbbm 1 - G_{0 ↑} \circ Δ_{↑}],

Z_{d ↑} ({t_{m}^{c}}, {t_{n}^{a}}) = z_{d ↑} ({t_{m}^{c}}, {t_{n}^{a}}) \cdot det [\mathbbm 1 - G_{0 ↑} \circ Δ_{↑}],

z_{d ↑} ({t_{m}^{c}}, {t_{n}^{a}}) = det [i G_{0 ↑}^{- 1}] = 1 + e^{- i \int_{C} d t [E_{d ↑} (t) + U (t) n_{↓} (t)]} .

z_{d ↑} ({t_{m}^{c}}, {t_{n}^{a}}) = det [i G_{0 ↑}^{- 1}] = 1 + e^{- i \int_{C} d t [E_{d ↑} (t) + U (t) n_{↓} (t)]} .

Z = C \sum w (C) .

Z = C \sum w (C) .

⟨ A ⟩ \equiv \frac{1}{Z} Tr [T_{C} e^{- i \int_{C} d t H (t)} A] = \frac{\sum _{C} A ( C ) w ( C )}{\sum _{C} w ( C )} = \frac{⟨ A sgn ⟩ _{MC}}{⟨ sgn ⟩ _{MC}},

⟨ A ⟩ \equiv \frac{1}{Z} Tr [T_{C} e^{- i \int_{C} d t H (t)} A] = \frac{\sum _{C} A ( C ) w ( C )}{\sum _{C} w ( C )} = \frac{⟨ A sgn ⟩ _{MC}}{⟨ sgn ⟩ _{MC}},

sgn (C) = \frac{w ( C )}{∣ w ( C ) ∣}

sgn (C) = \frac{w ( C )}{∣ w ( C ) ∣}

⟨ A ⟩_{MC} = \frac{\sum _{C} A ( C ) ∣ w ( C ) ∣}{\sum _{C} ∣ w ( C ) ∣} .

⟨ A ⟩_{MC} = \frac{\sum _{C} A ( C ) ∣ w ( C ) ∣}{\sum _{C} ∣ w ( C ) ∣} .

n_{↓} (t) = ⟨ n_{↓} (t; {t_{m}^{c}}, {t_{n}^{a}}) ⟩ .

n_{↓} (t) = ⟨ n_{↓} (t; {t_{m}^{c}}, {t_{n}^{a}}) ⟩ .

n_{↑} (t) = ⟨ 1 - [\mathbbm 1 + U (t, 0_{+}) U (- i β, t)]_{00}^{- 1} ⟩ .

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Code & Models

Repositories

patryk-kubiczek/noneq-cthalfhyb
noneOfficial

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Exact real-time dynamics of single-impurity Anderson model from a single-spin hybridization-expansion

P. Kubiczek1*, A. N. Rubtsov2,3, A. I. Lichtenstein1

1 I. Institute for Theoretical Physics, University of Hamburg, Hamburg, Germany

2 Russian Quantum Center, Skolkovo, Moscow Region, Russia

3 Department of Physics, Lomonosov Moscow State University, Moscow, Russia

[email protected]

Abstract

**In this work we introduce a modified real-time continuous-time hybridization-expansion quantum Monte Carlo solver for a time-dependent single-orbital Anderson impurity model: CT-1/2-HYB-QMC. In the proposed method the diagrammatic expansion is performed only for one out of the two spin channels, while the resulting effective single-particle problem for the other spin is solved semi-analytically for each expansion diagram. CT-1/2-HYB-QMC alleviates the dynamical sign problem by reducing the order of sampled diagrams and makes it possible to reach twice as long time scales in comparison to the standard CT-HYB method. We illustrate the new solver by calculating an electric current through impurity in paramagnetic and spin-polarized cases.

**

1 Introduction
2 Method
2.1 Single-spin hybridization-expansion
2.2 Dressed local weight from discretized bath
2.3 Discretizing contour-time instead of bath
2.4 Measurement of observables
3 Results
3.1 Discrete representation of the bath
3.2 Current through impurity
3.3 Effect of magnetic field
3.4 Computational performance
4 Conclusion
A Code

1 Introduction

Continuous-time quantum Monte Carlo methods (CT-QMC) are widely used to solve impurity models in equilibrium [1]. Despite their successful application for uncovering short-time nonequilibrium dynamics [2, 3, 4, 5, 6, 7, 8], they suffer from an inherent dynamical sign problem preventing an access to longer times. There have been numerous attempts to tame the sign problem in real-time QMC simulations, such as explicit summation over Keldysh indices [9, 10, 11, 12], bold-line methods [13, 14, 15, 16, 8], or inchworm algorithm [17, 18, 19, 20, 21]. Alternatively, a deterministic iterative summation of path integrals is possible [22, 23].

Here we propose a method that alleviates the dynamical sign problem in case of the single-orbital Anderson impurity model (AIM). It is based on the continuous-time hybridization-expansion QMC (CT-HYB-QMC [24]) and can be viewed as a kind of a bold-line QMC. In our method, called CT-1/2-HYB, only the diagrams for one out of two spins need to be sampled (hence the name). The contribution from the other spin can be summed up semi-analytically for each expansion term. The advantage of such an approach, which to our best knowledge has not yet been applied neither to equilibrium nor non-equilibrium calculations, is that the Monte Carlo average expansion order is decreased leading to the reduction of the sign problem. The problem is however that the semi-analytic evaluation of dressed diagrams for a continuous-bath is a non-trivial task. The solution we come up with in this work is to discretize the bath, which limits the spectral resolution and the timescale obtainable within our method.

The paper is organized as follows. In Section 2 we describe CT-1/2-HYB method in detail introducing two possible ways of semi-analytical evaluation of the bold diagrams. In Section 3 we present the strategy for discretizing the bath and results for a current through impurity in a paramagnetic and a spin-polarized case, as well as discuss the computational performance of the algorithm. In section 4 we summarize our findings and draw conclusions.

2 Method

2.1 Single-spin hybridization-expansion

The system for which CT-1/2-HYB is designed is the single-orbital AIM with general time-dependent parameters. Since the test case we consider in this paper is electric current calculation, it is convenient to introduce a setup with two baths, where each bath (lead) is labelled by index $\alpha$ . The Hamiltonian is composed of the local, bath and hybridization parts and reads

[TABLE]

The method of solution is based on the perturbative expansion of the dynamical partition function $Z$ with respect to hybridization strength on the Keldysh contour $\mathcal{C}$ with an additional imaginary time branch, defined by inverse temperature $\beta$ and maximum evolution time $t_{\mathrm{max}}$ (Fig. 1). The imaginary time branch enables us to start the time-evolution from an equilibrium thermal state at $t=0$ . An accessible introduction to this formalism containing imaginary- and real-time evolution can be found e.g. in [25].

Contrary to CT-HYB, the expansion is performed only for one spin species. We will be referring to spin-down as the expansion channel, which is an arbitrary choice made for a convenience of the discussion. In the following we work in the Schrödinger picture: contour-time arguments of operators serve only as contour labels (they do not denote Heisenberg-like evolution). Thus, we expand $Z$ in $H_{\mathrm{hyb}\downarrow}$ obtaining

[TABLE]

The contribution of each diagram, specified by expansion order $k$ along with sets of impurity-electron creation $\{t^{\mathrm{c}}_{m}\}$ and annihilation times $\{t^{\mathrm{a}}_{m}\}$ , can be represented as a product of the spin-down bath weight

[TABLE]

where $C^{\dagger}_{\downarrow}(t)\equiv\sum_{p}V_{p\downarrow}(t)c^{\dagger}_{p\downarrow}$ , and the local weight dressed by spin-up hybridization

[TABLE]

The spin-down bath weight can be readily evaluated following the procedure known from CT-HYB method [1]. Defining the hybridization function

[TABLE]

where $g_{\alpha p\sigma}(t,t^{\prime})$ is a free-electron Green function [25]

[TABLE]

we obtain

[TABLE]

where $Z_{\mathrm{bath}\downarrow}=\mathrm{Tr_{c_{\downarrow}}}\left[\mathcal{T}_{\mathcal{C}}\,e^{-i\int_{\mathcal{C}}dtH_{\mathrm{bath}\downarrow}(t)}\right]$ is a diagram-independent constant. Namely, for each diagram one needs to evaluate a determinant of a $k\times k$ spin-down hybridization matrix.

2.2 Dressed local weight from discretized bath

The novel aspect of CT-1/2-HYB is the presence of the spin-up hybridization dressing in the local weight, Eq. (7), which incorporates both local and spin-up bath degrees of freedom. This auxiliary system is governed by the Hamiltonian $H_{0}=~{}H_{\mathrm{loc}}+~{}H_{\mathrm{bath}\uparrow}+~{}H_{\mathrm{hyb}\uparrow}$ and is effectively non-interacting since $H_{0}$ conserves $n_{\downarrow}=\langle d^{\dagger}_{\downarrow}d_{\downarrow}\rangle$ . The contour-time dependence of $n_{\downarrow}$ is specified by the the locations of operators $d^{\dagger}_{\downarrow}$ : $\{t^{\mathrm{c}}_{m}\}$ and $d_{\downarrow}$ : $\{t^{\mathrm{a}}_{n}\}$ . Because $n_{\downarrow}$ can equal only 0 or 1, the only non-vanishing contributions to the dressed local trace come from the following orderings of the spin-down operators:

•

$-i\beta>t^{\mathrm{c}}_{k}>t^{\mathrm{a}}_{k}>t^{\mathrm{c}}_{k-1}>\ldots>t^{\mathrm{a}}_{2}>t^{\mathrm{c}}_{1}>t^{\mathrm{a}}_{1}>0_{+}$ (permutation sign: $1$ ),

•

$-i\beta>t^{\mathrm{a}}_{k}>t^{\mathrm{c}}_{k}>t^{\mathrm{a}}_{k-1}>\ldots>t^{\mathrm{c}}_{2}>t^{\mathrm{a}}_{1}>t^{\mathrm{c}}_{1}>0_{+}$ (permutation sign: $(-1)^{k}$ ).

In the following we assume that only those two orderings of contour-times $\{t^{\mathrm{c}}_{m}\},\{t^{\mathrm{a}}_{n}\}$ are allowed. Note that for $k=0$ both contributions have to be summed, corresponding to constant $n_{\downarrow}=0$ or 1 on the whole contour.

The dressed local weight can be represented as

[TABLE]

The factor $\varphi$ takes into account the permutation sign and the phases resulting from tracing over impurity spin-down occupations. If we denote the effective contour-time-dependent spin-down impurity occupation by $n_{\downarrow}(t)\equiv n_{\downarrow}(t;\{t^{\mathrm{c}}_{m}\},\{t^{\mathrm{a}}_{n}\})$ then

[TABLE]

The determinant part represents the trace of the evolution operator of the spin-up impurity and the spin-up bath subject to an effective contour-time-dependent single-particle Hamiltonian $h_{\uparrow}(t)$

[TABLE]

The single particle evolution operator $\mathcal{U}$ is given by

[TABLE]

and the relation between the trace of the many-body evolution operator and the single-particle evolution operator is provided by the identity valid for arbitrary fermions $f$ and a non-interacting Hamiltonian $h(t)$

[TABLE]

Note that this identity is also used in the determinant QMC [26], which is similarly based on the summation of determinants of matrices in the single particle basis. It is known that care must be taken while evaluating determinants of the form of Eq. (15) for large $\beta$ due to exponentially large and small matrix elements. We followed the advice from [27] to stabilize the necessary matrix multiplications and evaluation of determinants.

Because the contour-time dependence of $n_{\uparrow}(t)$ is step-wise constant we can take advantage of a convenient segment picture [1] to obtain a more practical representation of $\phi$ and $\mathcal{U}$ . We call a segment a time interval during which $n_{\downarrow}=1$ , i.e. after a creation time and before an annihilation time. Then we have

[TABLE]

and

[TABLE]

where $u_{n_{\downarrow}}$ is a single-particle time-evolution matrix

[TABLE]

2.3 Discretizing contour-time instead of bath

Another approach to evaluate the dressed local weight is to discretize the contour-time instead of the bath. Within the path integral-formalism one can integrate out the spin-up bath from Eq. (7) and work directly with the hybridization function $\Delta_{\uparrow}$ defined by Eq. (8). After calculating traces over the bath and spin-down impurity degrees of freedom we obtain

[TABLE]

where $\varphi$ is the same as in Eq. (12) and (16), $Z_{\mathrm{bath}\uparrow}$ is a constant and $Z_{d\uparrow}$ is given by the path integral

[TABLE]

with

[TABLE]

and where $n_{\downarrow}(t)$ is determined by the sets of contour-times $\{t^{\mathrm{c}}_{m}\},\{t^{\mathrm{a}}_{n}\}$ . When $\mathcal{G}_{\uparrow}$ is understood as a matrix in time domain, the path integral $Z_{d\uparrow}$ is given by the formula [28]

[TABLE]

Because it is undesirable to discretize the time derivative operator we can recast this formula using Eq. (20)

[TABLE]

where $\circ$ denotes contour-time convolution and $z_{d\uparrow}$ is a generalized partition function for a system with time-local action which can be explicitly evaluated

[TABLE]

For each diagram a generalized Green function $\mathcal{G}_{0\uparrow}$ should be calculated as a solution to Eq. (21). However, this has to be done for several contour-time spacings $\Delta t$ as an extrapolation of the result down to $\Delta t=0$ is required. Provided the convolution in Eq. (23) is done carefully as described in [29] the error of the discrete approximation should scale as $\Delta t^{2}$ .

We have decided not to follow this method due to its high computational cost. Nevertheless, this approach might be of a future interest since it allows one for an introduction of a retarded interaction $U_{\uparrow\downarrow}(t,t^{\prime})$ by adding an additional term $-U_{\uparrow\downarrow}(t,t^{\prime})n_{\downarrow}(t)$ to the left hand side of Eq. (20), which is only a small numerical effort.

2.4 Measurement of observables

The partition function $Z$ from Eq. (5) can be viewed as a sum of the weights $w(C)$ of diagrams $C=\{k;\{t^{\mathrm{c}}_{m}\},\{t^{\mathrm{a}}_{n}\}\}$ for all expansion orders $k$ and contour-times $\{t^{\mathrm{c}}_{m}\},\{t^{\mathrm{a}}_{n}\}$

[TABLE]

The space of diagrams $C$ is sampled within a Markov-chain Monte Carlo simulation. We have implemented the segment removal and addition moves, which guarantee the ergodicity of the sampling algorithm, and the shift moves (shifting the beginning or the end of a segment), which decrease the autocorrelation time. A detailed description of those moves along with Metropolis-Hastings acceptance rates can be found in [1]. We use fast updates for determinants of the hybridization matrices (bath weight) but calculate the dressed local weight explicitly for each Monte Carlo diagram.

The diagram weights discussed above can be in general any complex numbers. To work around this issue, we sample the diagrams with probabilities given by the absolute value of weights while simultaneously measuring their average (complex) sign. The expectation value of an observable $A$ is given by

[TABLE]

where $A(C)$ is the expectation value of $A$ in diagram $C$ , the complex sign is defined as

[TABLE]

and the Monte Carlo average follows the prescription

[TABLE]

It is known that CT-HYB-QMC for a single-impurity Anderson model does not suffer from the fermionic sign problem, i.e. the average sign on the imaginary axis is always one. However, the presence of real-time branches causes the dynamical sign problem. For each diagram having operators on the real branch there exists another one which has exactly the same weight but the opposite sign. This results from the fact that operator closest to $t_{\mathrm{max}}$ on the real axis can be moved freely from the upper to the lower branch (or conversely) and this operation changes only the sign of the diagram. This means that the average complex sign $\left<\mathrm{sgn}\right>$ is a real quantity which equals the probability of having no operators on real-time branches. Practically, the average sign decreases exponentially with $t_{\mathrm{max}}$ which effectively prohibits an accurate evaluation of observables using Eq. (26) for large $t_{\mathrm{max}}$ . To overcome the inaccuracy of the results we have to scale the total number of Monte Carlo samples exponentially with $t_{\mathrm{max}}$ until it is feasible.

In the following we present explicit measurement formulas for certain observables of interest.

The time-dependent impurity spin-down occupation is the average of the diagram-dependent function $n_{\downarrow}(t;\{t^{\mathrm{c}}_{m}\},\{t^{\mathrm{a}}_{n}\})$ where $t$ can lie either on the upper or the lower real-time branch

[TABLE]

The measurement of spin-up observables requires the evaluation of the elements of occupation matrix for the system of coupled spin-up impurity and bath. Namely, it can be shown

[TABLE]

This result follows from the general formula analogous to the one found in [26]

[TABLE]

The spin $\sigma$ electron current from lead $\alpha$ to the impurity is given by

[TABLE]

which can be obtained by solving the Heisenberg equation of motion for the impurity particle number operator $d^{\dagger}_{\sigma}d_{\sigma}$ and employing the continuity equation

[TABLE]

While evaluation of the expectation value of the current operator requires a reference solution within the standard CT-HYB algorithm (see e.g. Eq. (60) in [5]), it is possible to extract the full spin-up single-particle occupation matrix directly in CT-1/2-HYB and by that all the relevant matrix elements for the spin-up current. Similarly to (30)

[TABLE]

where $\alpha p$ is the index of the column corresponding to $p$ th electronic level from bath $\alpha$ . In practice, one measures the whole spin-up single-particle occupation matrix as defined in (30) and (35), from which the observables of interest can be constructed.

3 Results

3.1 Discrete representation of the bath

Before presenting the results, we show in detail the parameters of the bath we have chosen for our simulation. We want to model a bath defined by the following time- and spin-independent coupling function

[TABLE]

which represents a flat band of half-width $D=4$ with a soft cutoff controlled by parameter $\nu=3$ , coupled to the impurity with the coupling strength $\Gamma=1$ . However, within the discrete bath approach, only a peak-structured coupling function can be constructed

[TABLE]

In order to approximate the continuous bath-impurity coupling function, Eq. (36), by a discrete one, Eq. (37), we apply the equal weight method from Ref. [30], obtaining equal hybridization parameters $V_{p}$

[TABLE]

The resulting discrete representations of the bath for several values of bath size $N$ are presented in Fig. 2.

3.2 Current through impurity

We measure the electric current in a simulation of a voltage quench from $\phi(t=0)=0$ to $\phi(t>0)=2$ . In order to analyze the convergence of result with increasing bath size, we perform QMC calculations for several numbers $N$ of bath sites. We set inverse temperature $\beta=5$ and Coulomb interaction on the impurity $U=4$ . We consider the half-filled case, such that the (time-independent) local energy level

[TABLE]

In Fig. 3 we plot the spin-up electric current through impurity as a function of time. Due to the particle-hole symmetry the left hand side of Eq. (34) vanishes, thus currents from the left and right are equal up to the sign. Also, in the paramagnetic case the total current can be recovered by multiplication by factor of 2. For reference we also provide a result of an exact-diagonalization calculation for $N=6$ which confirms the validity of our method.

As also observed in a similar study [18], the current relatively quickly (around 1 unit of $\Gamma^{-1}$ ) reaches its steady state value. We note a quick convergence of the result with $N$ . The difference between results for $N=18$ and $N=24$ is only minimal.

It is instructive to look at the comparison of the average expansion order and the average Monte Carlo sign between CT-1/2-HYB and CT-HYB, which we present in Fig. 4 for $N=12$ (the dependence of those quantities on $N$ is negligible). In the paramagnetic case CT-1/2-HYB leads to the reduction of the average expansion order by a factor of two since only spin-down diagrams have to be sampled while spin-down and spin-up diagrams are contributing to $Z$ with equal weights. Noteworthy, the expansion order grows linearly with time and the rate of growth is smaller by a factor of 2. According to the established exponential decrease of the average sign with the growing expansion order [5, 4] twice as long time scales can be reached within CT-1/2-HYB method as compared to CT-HYB, when the smallest average sign feasible for obtaining the result is fixed.

3.3 Effect of magnetic field

A valid research question is whether the CT-1/2-HYB method can lead to an even better performance improvement over CT-HYB in presence of a magnetic field breaking the SU(2) symmetry of the model. In such a case, the average expansion orders in the two spin channels can differ, and hence one can choose to expand $Z$ in the channel with the lower expansion order. In the following, we modify the definition of the local energies $\mathcal{E}_{d\sigma}$ such that it takes into account non-zero chemical potential $\mu$ and magnetic field $B$

[TABLE]

Fig. 5(a) shows the average expansion orders for the imaginary-time (equilibrium) CT-1/2-HYB as a function of $B$ for several values of $\mu$ . The maximum spin-down average order occurs when the spin-down impurity occupation is close to 0.5, which is evident from Fig. 5(b). For further analysis, it helps to realize that the following relation holds

[TABLE]

where $\left<k_{\uparrow}\right>$ is a hypothetical average expansion order in case the hybridization expansion was performed in the spin-up channel. Then, we define the improvement factor $\eta$ as

[TABLE]

where the sum $\left<k_{\downarrow}\right>+\left<k_{\uparrow}\right>$ is a hypothetical average expansion order in the standard CT-HYB algorithm. We plot the equilibrium improvement factor as a function of $|B|$ in Fig. 5(c). In the paramagnetic case the improvement factor $\eta$ is always 2. For $\mu=0$ the improvement factor also equals 2 and does not depend on $B$ , which results from the particle-hole symmetry. Only by simultaneous introduction of non-zero $\mu$ and $B$ the average expansion orders can be made different in the two spin channels, leading to $\eta>2$ .

While Fig. 5(c) proves that CT-1/2-HYB can more than halve the average expansion order on the imaginary axis, it is not yet clear if the same applies to the real-time evolution. Ultimately, it is the growth of the expansion order with the evolution time that is causing the dynamical sign problem. That is why we decided to estimate this quantity for the calculation of time-dependent current in the presence of chemical potential and magnetic field. We choose $\mu=2$ and $B=2,-2$ , which is a case where an improvement factor in equilibrium of around 2.5 can be achieved.

In the polarized case currents from the right and from the left bath are different. Fig. 6 presents the left and the right spin-up currents in the presence of magnetic fields $B=-2,2$ . The spin-up current is lower for $B=2$ since then the spin-up impurity level lies almost at the edge of the bath which hinders transport dynamics.

Fig. 7 shows the growth of the expansion order with time for the paramagnetic case compared with the two polarized cases. Contrary to the imaginary-time calculation, the expansion order grows at almost the same rate irrespective of chemical potential and magnetic field. Noteworthy, for $B=-2$ the growth rate of the average expansion order is slightly higher than for other cases, even though the initial average order (for $t=0$ ) is the lowest in this case.

The different behaviors of the average expansion orders in imaginary and real-time expansions were noted previously in [5]. The average expansion order on the Keldysh contour has no interpretation as a quantity proportional to the kinetic energy and is thus only weakly influenced by the position of the impurity level relative to the conduction bath.

3.4 Computational performance

The biggest computational burden of the method comes from the need to overcome the dynamical sign problem. For that reason we have scaled the number of Monte Carlo steps exponentially with $t_{\mathrm{max}}$ in order to keep the error approximately constant for all times. We have been able to generate enough samples such that the statistical errors are smaller than size of markers on all the plots presented in this work.

The dependence of the average sign on the physical parameters is the same as in the standard CT-HYB method on a Keldysh contour. Namely, keeping the hybridization strength as energy unit, the sign is largely unaffected by interaction strength, chemical potential or magnetic field but decreases with bath bandwidth [4].

Now we discuss the computational complexity of diagram evaluations. For the discretized bath approach the number of bath sites $N$ plays the most important role. Generally, the expected scaling of this approach is $\mathcal{O}(\left<k_{\uparrow}\right>N^{3})$ and results from multiplications single-particle evolution matrices of size $(N+1)\times(N+1)$ , whose average number scales linearly with $\left<k_{\uparrow}\right>$ . The actual scaling can deviate from this prediction as certain implementation-dependent optimizations might be possible.

For the discretized time approach the number of contour time steps $N_{t}$ determines the computational scaling. Both the solution of the Dyson equation [25] and the evaluation of the determinant in time-space will scale as $\mathcal{O}(N_{t}^{3})$ .

Moreover, it takes on average $\left<k_{\uparrow}\right>$ Monte Carlo moves to decorrelate subsequent measurements. This implies an additional scaling factor of $\left<k_{\uparrow}\right>$ [13]. The average order $\left<k_{\downarrow}\right>$ is linear in inverse temperature $\beta$ and $t_{\mathrm{max}}$ : $\left<k_{\downarrow}\right>=c_{\beta}\,\beta+c_{t}\,t_{\mathrm{max}}$ . Leaving the dynamical sign problem aside and neglecting the cost of bath weight evaluation, the overall theoretical complexity of the discretized bath approach is $\mathcal{O}(\left<k_{\downarrow}\right>^{2}N^{3})$ and that of the discretized contour-time approach is $\mathcal{O}(\left<k_{\downarrow}\right>N_{t}^{3})$ (where $N_{t}$ is the number of time-steps, including those on the imaginary branch).

The computational scaling of our implementation of diagram evaluations can be learned from Fig. 8. The benchmark problem is the evaluation of the paramagnetic current from section 3.2, where discretized bath approach has been utilized.

Fig. 8(a) compares CPU time per Monte Carlo (MC) step in the CT-HYB algorithm against the CT-1/2-HYB results for different numbers $N$ of bath sites. The cost of a single MC step in CT-HYB algorithm is lower for all the considered values of $N$ , pointing to a huge computational effort related to the evaluations of the dressed local weight. The approximately linear growth of the CPU time per step with $t_{\mathrm{max}}$ is due to an increasing average order $\left<k_{\downarrow}\right>$ .

It is also constructive to understand the scaling of the necessary computational resources with growing $N$ . In Fig. 8(b) we present the CPU time per MC step as a function of $N$ for a fixed $t=1$ . The fitted computational complexity is $\mathcal{O}(aN^{2}+bN^{3})$ with $a\gg b$ .

The most general version of our implementation of CT-1/2-HYB involves a precomputation of all $u_{0}(t,t^{\prime})$ and $u_{1}(t,t^{\prime})$ on a two-dimensional time-grid and subsequent interpolations during a Monte Carlo run. Unfortunately, this approach can quickly lead to very high memory demands, especially for large $N$ . However, when the time-dependence of Hamiltonian parameters is step-wise, which is the case for the voltage quench, the time-evolution generated by the operators $u_{n_{\downarrow}}$ can be performed in their eigenbases. For this purpose, one needs to store only a small set of eigenenergies and rotation matrices, which not only solves the problem of high memory requirement but also avoids the interpolation-related numerical inaccuracy and presumably improves performance. We have applied this strategy for the calculations presented here, which possibly explains the $\mathcal{O}(N^{2})$ contribution to the scaling due to presence of diagonal matrices.

4 Conclusion

CT-1/2-HYB-QMC allows one to solve the single-orbital Anderson impurity model exactly on twice as longer time scales than standard CT-HYB-QMC due to a reduction of the average expansion order. The growth rate of the average expansion order with $t_{\mathrm{max}}$ is independent of the location of impurity spin-up and spin-down levels in the considered time-evolution after the voltage quench, which rules out an even better performance gain in cases with non-zero chemical potential and magnetic field.

It is possible to extend this approach to multi-orbital impurity models with density-density interactions. If all the orbitals are interacting, the improvement factor will diminish with the number of fermionic flavors (combined spin and orbital numbers) $n$ as $\frac{n}{n-1}$ . This results from the fact that the semi-analytical evaluation of the diagrams will be possible for only one out of $n$ flavors. Moreover, in this case the method can only be applied to Hamiltonians conserving orbital and spin quantum numbers, which poses a serious limitation to the modelling of realistic materials. However, if the interaction matrix is sparse, also some single-particle mixing terms can be allowed between non-interacting flavors, while the improvement factor can be better than $\frac{n}{n-1}$ as more flavors can be treated semi-analytically.

We note that the cost of the CT-HYB inchworm algorithm [17] increases exponentially with number of orbitals due to an exponential increase in the size of many-body Hilbert space, which does not depend on the sparsity of interactions. On the other hand, in CT-1/2-HYB the increase in computational cost with the number of orbitals is in principle only polynomial. The main challenge of a real-time QMC is however the dynamical sign problem, which still leads to an exponential CPU-time scaling with $t_{\mathrm{max}}$ in CT-1/2-HYB and only to a polynomial scaling in case of inchworm algorithm. Even though CT-1/2-HYB-QMC may significantly extend accessible timescales for multi-orbital models with sparse interactions, it is still subject for future analysis to understand whether it will be enough to study any phenomena of interest.

It is worthwhile to consider whether the CT-1/2-HYB could be combined with bold or inchworm QMC. Their strength is that dressed many-body propagators can be used to evaluate Monte Carlo diagrams. Those propagators are actually of no use in CT-1/2-HYB as an attempt to include many-body corrections to CT-1/2-HYB would undermine the basic assumption about the effective single-particle dynamics. Nevertheless, the inchworm expansion could be done around propagators dressed by spin-up processes in the spirit of CT-1/2-HYB. In contrast to pure CT-1/2-HYB, some spin-up diagrams would be missed in such an inchworm procedure and would still have to be sampled separately. It is unclear at the moment how the efficiency of the inchworm algorithm would be affected by use of dressed propagators since no such attempts have been reported.

An interesting single-orbital application of CT-1/2-HYB method as an impurity solver for DMFT could be investigation of the cross-over behavior between Falicov-Kimball and Hubbard models, by tuning the spin-down hopping from zero to the value of the spin-up hopping. Our method is especially suitable for studying this problem since the zeroth order term of 1/2-HYB expansion corresponds to Falicov-Kimball-like impurity.

Acknowledgements

We acknowledge the use of Armadillo library [31] for linear-algebra operations in the core of CT-1/2-HYB code and of QuSpin package for ED benchmark [32]. The calculations were performed on the HLRN-4 supercomputer of the North-German Supercomputing Alliance.

Funding information

The work has been supported by SFB 925 “Light induced dynamics and control of correlated quantum systems” financed by DFG (P.K., A.I.L.) and by a joint DFG/RNF grant no. 16-42-01057/LI 1413/9-1 (A.N.R., A.I.L.).

Appendix A Code

Our open-source implementation of the CT-1/2-HYB-QMC impurity solver is available at https://github.com/patryk-kubiczek/noneq-cthalfhyb.

Bibliography32

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1[1] E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer and P. Werner, Continuous-time Monte Carlo methods for quantum impurity models , Rev. Mod. Phys. 83 , 349 (2011), 10.1103/revmodphys.83.349 . · doi ↗
2[2] L. Mühlbacher and E. Rabani, Real-Time Path Integral Approach to Nonequilibrium Many-Body Quantum Systems , Phys. Rev. Lett. 100 , 176403 (2008), 10.1103/Phys Rev Lett.100.176403 . · doi ↗
3[3] M. Schiró and M. Fabrizio, Real-time diagrammatic Monte Carlo for nonequilibrium quantum transport , Phys. Rev. B 79 , 153302 (2009), 10.1103/Phys Rev B.79.153302 . · doi ↗
4[4] M. Schiró, Real-time dynamics in quantum impurity models with diagrammatic Monte Carlo , Phys. Rev. B 81 , 085126 (2010), 10.1103/physrevb.81.085126 . · doi ↗
5[5] P. Werner, T. Oka and A. J. Millis, Diagrammatic Monte Carlo simulation of nonequilibrium systems , Phys. Rev. B 79 , 035320 (2009), 10.1103/Phys Rev B.79.035320 . · doi ↗
6[6] M. Eckstein, M. Kollar and P. Werner, Thermalization after an Interaction Quench in the Hubbard Model , Phys. Rev. Lett. 103 , 056403 (2009), 10.1103/physrevlett.103.056403 . · doi ↗
7[7] P. Werner, T. Oka, M. Eckstein and A. J. Millis, Weak-coupling quantum Monte Carlo calculations on the Keldysh contour: Theory and application to the current-voltage characteristics of the Anderson model , Phys. Rev. B 81 , 035108 (2010), 10.1103/Phys Rev B.81.035108 . · doi ↗
8[8] A. E. Antipov, Q. Dong and E. Gull, Voltage Quench Dynamics of a Kondo System , Phys. Rev. Lett. 116 , 036801 (2016), 10.1103/physrevlett.116.036801 . · doi ↗

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Code & Models

Videos

Abstract

Contents

1 Introduction

2 Method

2.1 Single-spin hybridization-expansion

2.2 Dressed local weight from discretized bath

2.3 Discretizing contour-time instead of bath

2.4 Measurement of observables

3 Results

3.1 Discrete representation of the bath

3.2 Current through impurity

3.3 Effect of magnetic field

3.4 Computational performance

4 Conclusion

Acknowledgements

Funding information

Appendix A Code