The auxiliary Hamiltonian approach and its generalization to non-local self-energies

0.8 MB·

by  Karsten Balzer

The auxiliary Hamiltonian approach and its
generalization to non-local self-energies
6
1
0 Karsten Balzer
2
Rechenzentrum, Christian-Albrechts-Universita¨t zu Kiel, Ludewig-Meyn-Strasse 4, 24118 Kiel,
n
Germany
a
J E-mail: [email protected]

4
2 Abstract. The recently introduced auxiliary Hamiltonian approach [Balzer K and Eckstein

M 2014 Phys. Rev. B 89 035148] maps the problem of solving the two-time Kadanoff-Baym
]
l equationsontoanoninteractingauxiliarysystemwithadditionalbathdegreesoffreedom. While
e
the original paper restricts the discussion to spatially local self-energies, we show that there
-
r existsaratherstraightforwardgeneralizationtotreatalsonon-localcorrelationeffects. Theonly
t
s drawbackisthelossoftimecausalityduetoacombinedsingularvalueandeigendecomposition
. ofthetwo-timeself-energy,theapplicationofwhichinhibitsonetoestablishtheself-consistency
t
a directlyonthetimestep. Forderivationandillustrationofthemethod,weconsidertheHubbard
m
modelinonedimensionandstudythedecayoftheN´eelstateintheweak-couplingregime,using
- the local and non-local second-order Born approximation.
d
n
o
c
[
1. Introduction

1
ManynumericalapproacheswhicharebasedontheKeldyshformalismrequirethesolutionofthe
v
5 two-timeKadanoff-Baymequations(orthenonequilibriumDysonequation)whichareequations
1 ofmotionsfortheone-particlenonequilibriumGreen’sfunctiondefinedontheL-shapedKeldysh
4
contour[1,2,3]. Theparticulardouble-timeandnon-Markovianformoftheseintegro-differential
6
equations generally renders their integration a non-trivial and challenging task. Moreover, big
0
. efforts are needed to reach either long times, handling a large memory kernel, or to simulate
1
0 quantumsystemswithmanydegreesoffreedomandespeciallyspatiallyinhomogeneoussystems.
6 Despitetheseobstacles, muchcomputationalprogresshasbeenmadeovertherecentdecades,
1 driven by the widespread success of nonequilibrium Green’s functions in various fields, including
:
v nuclear [4, 5, 6], plasma [7, 8], semiconductor [9, 10, 11], condensed matter [12, 13] and atomic
i
X andmolecularphysics[14,15]. FollowingthepioneeringworkofRefs.[4,5],whichfocusesonthe
numerical simulation of heavy-ion collisions, time propagation schemes of the Kadanoff-Baym
r
a equations (KBE) have been developed in different contexts to study the correlation build-up
in homogeneous fermion systems [16, 17]. A general Fortran code is presented in Ref. [18]
which efficiently evaluates the collision integrals (convolutions) by using fast Fourier transforms
of the momentum resolved Green’s function. Concerning the treatment of initial correlations
different pathways have been pursued which range from the use of adiabatic switching [19, 20]
and the embedding of additional collision terms derived from the equilibrium pair correlation
function [21] to the extension of the real-time Green’s function onto the imaginary branch of the
Keldysh contour [22, 23]. With the advance of computer power, also finite and inhomogeneous
systemshavebecomeaccessible. Acorrespondingwellrecognizedpredictor-corrector-basedtime
propagation code is described in Ref. [24] and is used in many applications. Further progress
is due to adapted basis sets which have been deployed to simplify the interaction (Coulomb)
matrix elements, e.g., [25], and due to efficient parallelization strategies [26, 27, 28], where the
Green’s function is distributed over the memory of multiple compute nodes in such a way that
the collision integrals can be calculated locally and a minimum communication is required to
perform the time stepping.
A feature, that is common to virtually all time propagation codes, is the direct evaluation of
the memory kernel, i.e., the computation of the collision integrals which cover the whole time
history. Togetherwiththedouble-timestructureoftheGreen’sfunction,thisleadstothetypical
n3-scalingofthenumericalalgorithmsinthenumberoftimesteps. Onlyfurtherapproximations
t
canovercomethisunfavorablescaling, suchastheapplicationofthegeneralizedKadanoff-Baym
ansatz[20,29,30],whichreconstructsthetwo-timeGreen’sfunctionfromitstime-diagonalvalue
andisexpectedtobeareasonableoptiontoextendstate-of-theartsimulations, e.g.,[31,32,33].
The prevailing similarity of the time propagation schemes raises the question whether there
exist alternative methods to solve the two-time KBE. Aside from matrix inversion techniques,
e.g. [34], a promising idea has been triggered by nonequilibrium dynamical mean-field theory
(DMFT)[13],whereHamiltonian-basedimpuritysolvershaverecentlybecomeavailablethrough
a two-time decomposition of the hybridization function [35, 36, 37, 38]. In Ref. [39] it
has been shown that a similar decomposition of the self-energy can be used to derive an
auxiliary Hamiltonian representation of the Kadanoff-Baym equations, where the one-particle
nonequilibrium Green’s function of the interacting many-body problem under consideration
is obtained from a noninteracting auxiliary system which couples to additional bath orbitals
the number of which is finite. Although progress has been made recently by showing that a
nonequilibriumself-energyhasauniqueLehmannrepresentation[40], theauxiliaryHamiltonian
method has so far been proven to be effective only for a local self-energy, which represents an
approximation becoming exact in the limit of infinite spatial dimensions (and is just the central
approximation in DMFT).
In this contribution, we show that the method can be extended rather easily to treat also
non-local self-energies and that it is thus applicable more generally. Particularly, it can improve
the n3-scaling in the same manner as in Ref. [39] as long as the computation of the self-energy
t
itself represents a sufficiently simple task. Note that this is a rather general prerequisite of the
approach, independent of the self-energy’s spatial form. We start the discussion with a brief
review of the auxiliary Hamiltonian approach for systems with a local self-energy (Sec. 2.1).
Thereafter, we derive an extended auxiliary system for one of the most simple test cases, the
(two-site)Hubbarddimer, andshowhowtheadditionallyemergingbathparametersfollowfrom
a singular value decomposition of the offdiagonal components of the self-energy (Sec. 2.2.1).
After numerical validation of the scheme in Sec. 2.2.2, we then generalize the approach to an
arbitrary number of sites and hence to a general quantum system (Sec. 2.3). Based on further
numerical results, which are presented in Sec. 3, we finally discuss additional simplifications of
the method.
2. Theory
For the discussion of the auxiliary Hamiltonian approach and its extension to non-local self-
energies, we consider the single-band Fermi-Hubbard model,
(cid:18) (cid:19)(cid:18) (cid:19)
(cid:88) (cid:88) 1 1
H = Tijc†iσcjσ +U ni↑− 2 ni↓− 2 , (1)
ijσ i
whereT arehoppingmatrixelementsconnectingthelatticesitesiandj,U isthelocalCoulomb
ij
repulsion, c†iσ and ciσ are creation and annihilation operators for particles of spin σ ∈ {↑,↓},
(i−1)s is (i+1)s
U U U
... T T ... =⇒ J(i−1)sσ Jisσ J(i+1)sσ
i 1 i i+1
−
... ...
T T
i 1 i i+1
−
interacting noninteracting
Figure 1. Construction of the auxiliary Hamiltonian in the case of a local self-energy for
a Hubbard chain with nearest-neighbor hopping T and on-site interaction U. While the filled
circles denote the original lattice sites, the open circles mark the bath orbitals i which are
s
coupledindependentlyofoneanothertoalatticesiteibyatime-and(ingeneral)spin-dependent
hopping matrix element J (t).
isσ
and niσ = c†iσciσ is the density. Using the nonequilibrium Green’s function Gijσ(t,t(cid:48)) =
(cid:82)
c−oin(cid:104)TtoCucriσ-o(tr)dce†jrσi(ntg(cid:48))(cid:105)o,pweritahto(cid:104)rTTC.,..t(cid:105)h=e ttirm[TeCeexvpol(uSt)io.n..]o/ftrt[hTeCesxypst(eSm)],(a1c)tiisongiSve=n −byi tChdesKHa(dsa)naonffd-
C
Baym equation
(cid:90)
(cid:88) (cid:88)
[i∂ δ T ]G (t,t) = δ (t,t)δ + dsΣ (t,s)G (s,t), (2)
t ik ik kjσ (cid:48) (cid:48) ij ikσ kjσ (cid:48)
− C
k k C
which is accompanied by its adjoint equation with t t. For simplicity, we incorporate the
(cid:48)
↔
Hartree potential V (t) = U( n (t) 1) into the hopping matrix, thus we define as Σ (t,t)
i (cid:104) iσ¯ (cid:105)− 2 ijσ (cid:48)
only the correlation part of the self-energy. As the unit of energy (time) we choose the (inverse)
hopping T.
2.1. Auxiliary Hamiltonian for a local self-energy
Following Ref. [39], the Kadanoff-Baym equations (2) can be mapped onto a noninteracting
auxiliary system when the self-energy is local in space, i.e., if Σ (t,t) δ Σ (t,t). In this
ijσ (cid:48) ij iσ (cid:48)
∝
case the retardation effects described by Σ can be mimicked by sets of bath orbitals which
iσ i
S
are coupled to each individual lattice site. The corresponding auxiliary Hamiltonian is given by
(cid:88) (cid:88)(cid:88)
Haux(t) = Tijc†iσcjσ + (cid:15)isσ(t)a†isσaisσ (3)
ijσ iσ s
(cid:88)(cid:88)(cid:16) (cid:17)
+ Jisσ(t)a†isσciσ +Ji∗sσ(t)c†iσaisσ ,
iσ s∈Si
where a†isσ and aisσ describe the creation and annihilation of particles with spin σ in a bath
orbital i with energy (cid:15) (t) that is assigned to the lattice site i, see Fig. 1. By comparing the
s isσ
KBE of the auxiliary system (3) and the original lattice system (1) [Eq. (2)], one finds that the
auxiliarysystemhasthesameGreen’sfunctionG (t,t)(withi,j beinglatticeindices)provided
ijσ (cid:48)
that for all times t and t on the Keldysh contour the bath parameters can be determined from
(cid:48)
(cid:88)
Σ (t,t) = J (t)g((cid:15) ;t,t)J (t), (4)
iσ (cid:48) isσ isσ (cid:48) i∗sσ (cid:48)
s∈Si
where g((cid:15);t,t) is the Green’s function of an isolated bath orbital,
(cid:48)
g((cid:15);t,t(cid:48)) = i[fβ((cid:15)(0)) θ (t,t(cid:48))]e−i(cid:82)tt(cid:48)ds(cid:15)(s), (5)
− C
with f ((cid:15)) = 1/(eβ(cid:15) +1) being the Fermi-Dirac distribution for an inverse temperature β, and
β
θ the Heavyside step function on the contour.
C
Analyzing Eq. (4) regarding the location of the time arguments on the contour, one can
(1) (2)
distinguish two classes = of bath orbitals that model different dynamical effects:
Si ∪Si Si
(1)
(i) Forthemixedcomponents(Σ(cid:101) )oftheself-energy,thebathorbitalsin describethetime
iσ Si
evolutionofcorrelationsthatarepresentintheinitialstate. Asthesecontributionstypically
decay as function of time, the corresponding bath orbitals become decoupled (J (t) 0)
isσ →
in the limit t . In general, the bath parameters of this class depend on the generalized
→ ∞
spectral function
1 (cid:16) (cid:17)
C(cid:101) (t,(cid:15)) = Σ(cid:101) (t,(cid:15)+i0) Σ(cid:101) (t,(cid:15) i0) . (6)
iσ 2π iσ − iσ −
Details on the construction of the parameters are given in Ref. [35].
(ii) For the real-time greater (Σ>) and lesser (Σ<) components of the self-energy, the bath
iσ iσ
(2)
orbitals in describe the dynamical build-up of correlations. If we choose the initial
Si
bath energies (cid:15) (0) [at time t = 0] such that the Fermi-Dirac distribution f ((cid:15) (0)) is
isσ ≷ β isσ
either 0 or 1, the Green’s function g evaluates to g ((cid:15) ;t,t) = i, and one can decompose
isσ (cid:48) ∓
both real-time components separately,
(cid:88)
iΣ<(t,t) = J (t)J (t), (7)
− iσ (cid:48) isσ i∗sσ (cid:48)
s <
∈Si
(cid:88)
iΣ>(t,t) = J (t)J (t), (8)
iσ (cid:48) isσ i∗sσ (cid:48)
s >
∈Si
where < > = (2).
Si ∪Si Si
(2)
The form of Eqs. (7) and (8) reveals that the bath parameters J (t) in can be obtained
isσ Si
directlyfromamatrixdecompositionofthetwo-timeself-energy. Astheleft-handsidesrepresent
Hermitian and positive definite matrices (in t and t), one can either apply an eigenvalue1 or
(cid:48)
a Cholesky decomposition, see [35, 39] for details. We emphasize that an exact representation
requires the number of bath orbitals to be equal to the number of time steps. To drastically
reduce the size of the baths one however can use low-rank representations, e.g., by performing
the sums in Eqs. (7) and (8) only over the important eigenvalues or by using a finite number of
bathorbitalstogetanexactCholeskydecompositiononthefirsttimestepswhilefittingthetime
dependency at later times [35]. Below, we consider only bath orbitals which fall into class (ii),
i.e., we start either from a noninteracting initial state, a mean-field state or the atomic limit. In
all these cases, the Matsubara and mixed components of the self-energy vanish (ΣM = Σ(cid:100) = 0),
iσ iσ
(1)
and consequently the bath orbitals in are redundant.
Si
As a result of the mapping procedure, the solution of the interacting many-body problem (1)
has been replaced by the determination of the Green’s function Gaux(t,t) of the noninteracting
ijσ (cid:48)
auxiliary system. The only assumption, at this stage, is the DMFT approximation of a local
1 Incaseofaneigenvaluedecomposition,thesquarerootsoftherealeigenvaluesareintegratedintothehopping
parameters.
self-energy. A self-consistent solution is further prepared by iteration, starting, e.g., from the
noninteracting Green’s function. Moreover, we can use the causal property of the Cholesky
decomposition (where an increase of the matrix dimension does not alter the decomposition at
earlier times) to get the self-consistent bath parameters stepwise directly on the time step, see
in particular Ref. [35]. In practice, we also do not need to store the Green’s function. Instead
it is sufficient to save the self-energy which enters the construction of the bath parameters. For
the numerical computation of the auxiliary Green’s function, we have developed a Krylov-based
time-propagation scheme (see Appendix of Ref. [39]), where the greater (G> (t ,t)) and lesser
ijσ s
component ([G< (t,t )] ) of the Green’s function are calculated for all times t t on the
ijσ s † ≤ s
time step t , i.e., on the edge of the lower and upper propagation triangle, compare, e.g., with
s
Ref. [24]. In addition, OpenMP-based parallelization of the code is achieved by performing the
time stepping for all vectors G>(t ,t) with components G> (t ,t) [for all i] in parallel, thus for
jσ s ijσ s
an infinitesimal time step δt and times t t (and analogously for the lesser component), we
s
≤
evaluate
G>(t ,t)
G>(t +δt,t) = G>(t ,t) exp[ ihaux(t)δt] jσ s , (9)
jσ s | jσ s | − σ G>(t ,t)
| jσ s |
where haux(t) denotes the one-particle Hamiltonian of the auxiliary system H (t). We note
σ aux
that the normalization of the vector G> on the right-hand side of Eq. (9) is a precondition to
jσ
apply the unitary time-evolution operator within the Krylov method [41].
2.2. Generalization of the approach to non-local self-energies
From the first point of view, it seems to be not a straightforward endeavor to take into account
arbitrary non-local self-energies in the auxiliary Hamiltonian approach. This is due to the fact
that any non-local, i.e., offdiagonal component of the self-energy does not describe retardation
effects on a single site which is in clear correspondence to a noninteracting site coupled to a
bath as discussed in Sec. 2.1. Moreover, offdiagonal components of Σ are less symmetric, which
σ
raises the question whether it is at all possible to construct an auxiliary system such that the
conservation laws for particle number, momentum and energy (which come with a conserving
approximation)areobeyed. Nevertheless,itisintuitivetoimaginethatcorrelationeffects,which
originate from non-local parts of the self-energy, should be representable by additional particle
motionthatconnectsdifferentsitesofthelattice. Inthefollowing,weshowthatavalidmapping
exists for the Hubbard dimer. How the findings can be generalized to more extended systems is
discussed afterwards (see Sec. 2.3).
2.2.1. Hubbard dimer The simplest lattice model, in which a non-local self-energy appears, is
the Hubbard dimer [42, 43], consisting of two connected sites 1 and 2. In order to account for
diagonal and offdiagonal components of the self-energy we extend the auxiliary Hamiltonian of
Eq. (3) by a set of bath orbitals [with energies (cid:15) (t)] which exchange particles with both sites
3sσ
of the dimer. If we refer to the corresponding hopping matrix elements as C (t) and D (t) and
sσ sσ
define J (t) = A (t) and J (t) = B (t) (compare with Fig. 2), the ansatz for the mapping
1sσ sσ 2sσ sσ
is given by
(cid:88)(cid:16) (cid:17) (cid:88) (cid:88)
Haux(t) =−T c†1σc2σ +c†2σc1σ + (cid:15)isσ(t)a†isσaisσ (10)
σ i∈{1,2,3}σs∈Si
(cid:88)(cid:16) (cid:17)
+ Asσ(t)a†1sσc1σ +B2sσ(t)b†sσc2σ +Csσ(t)a†3sσc1σ +Dsσ(t)a†3sσc2σ +H.c. ,
s
where in the last term s runs only over existing bath sites. To determine the parameters of
the baths 1 , 2 and 3 in the presence of all components of the self-energy, we formulate the
s s s
1 2
s s
3
s
U U
= A Csσ Dsσ B
T ⇒ sσ sσ
1 2
T
1 2
interacting noninteracting
Figure 2. Construction of the auxiliary Hamiltonian representation of the KBE for the
Hubbard dimer, including non-local correlation effects.
equations of motions for the one-particle Green’s function of the auxiliary system (10) and
compare them to the KBE of the Hubbard dimer (i,j 1,2 ),
∈ { }
(cid:90)
(cid:88)
[i∂ +Tδ ]G (t,t) δ δ (t,t) = dsΣ (t,s)G (s,t). (11)
t i j1 ijσ (cid:48) ij (cid:48) ikσ kjσ (cid:48)
|− | − C
k 1,2 C
∈{ }
For the auxiliary system we have
(i) the lattice components of the auxiliary Green’s functions, which evolve in time according to
i∂ Gaux δ δ (t,t) = TGaux(t,t)+(cid:88)(cid:2)A (t)Gaux (t,t)+C (t)Gaux (t,t)(cid:3) , (12)
t 1jσ − 1j C (cid:48) − 2jσ (cid:48) sσ 1sjσ (cid:48) sσ 3sjσ (cid:48)
s
i∂ Gaux δ δ (t,t) = TGaux(t,t)+(cid:88)(cid:2)B (t)Gaux (t,t)+D (t)Gaux (t,t)(cid:3) ,
t 2jσ − 2j C (cid:48) − 1jσ (cid:48) sσ 2sjσ (cid:48) sσ 3sjσ (cid:48)
s
where j 1,2 , and
∈ { }
(ii) the mixed bath-lattice terms, which obey
[i∂ (cid:15) (t)]Gaux = A (t)Gaux(t,t), (13)
t− 1sσ 1sjσ ∗sσ 1jσ (cid:48)
[i∂ (cid:15) (t)]Gaux = B (t)Gaux(t,t),
t− 2sσ 2sjσ s∗σ 2jσ (cid:48)
[i∂ (cid:15) (t)]Gaux = (cid:2)C (t)Gaux(t,t)+D (t)Gaux(t,t)(cid:3) .
t− 3sσ 3sjσ s∗σ 1jσ (cid:48) s∗σ 2jσ (cid:48)
AnalogouslytoRef.[39],themixedbath-latticetermscanbedeterminedbyusingtheGreen’s
function g((cid:15) ;t,t) of an isolated bath orbital (cf. Eq. (5)):
isσ (cid:48)
(cid:90)
Gaux (t,t) = dsg((cid:15) ;t,s)A (s)Gaux(s,t), (14)
1sjσ (cid:48) 1sσ ∗sσ 1jσ (cid:48)
(cid:90)C
Gaux (t,t) = dsg((cid:15) ;t,s)B (s)Gaux(s,t),
2sjσ (cid:48) 2sσ s∗σ 2jσ (cid:48)
(cid:90)C
Gaux (t,t) = dsg((cid:15) ;t,s)(cid:2)C (s)Gaux(s,t)+D (s)Gaux(s,t)(cid:3) .
3sjσ (cid:48) 3sσ s∗σ 1jσ (cid:48) s∗σ 2jσ (cid:48)
C
Inserting the results of Eq. (14) into Eq. (12) and comparing term by term to the Kadanoff-
Baym equation (11), we find that the offdiagonal components of the self-energy can indeed be
represented by the bath 3 . However, the expressions for the diagonal components do not only
s
involve the hopping matrix elements A (t) and B (t). Instead, we obtain
sσ sσ
(cid:88) (cid:88)
Σ (t,t) = A (t)g((cid:15) ;t,t)A (t)+ C (t)g((cid:15) ;t,t)C (t), (15)
11σ (cid:48) sσ 1sσ (cid:48) ∗sσ (cid:48) sσ 3sσ (cid:48) s∗σ (cid:48)
s 1 s 3
∈S ∈S
(cid:88) (cid:88)
Σ (t,t) = B (t)g((cid:15) ;t,t)B (t)+ D (t)g((cid:15) ;t,t)D (t),
22σ (cid:48) sσ 2sσ (cid:48) s∗σ (cid:48) sσ 3sσ (cid:48) s∗σ (cid:48)
s 2 s 3
∈S ∈S
(cid:88)
Σ (t,t) = C (t)g((cid:15) ;t,t)D (t),
12σ (cid:48) sσ 3sσ (cid:48) s∗σ (cid:48)
s 3
∈S
(cid:88)
Σ (t,t) = D (t)g((cid:15) ;t,t)C (t).
21σ (cid:48) sσ 3sσ (cid:48) s∗σ (cid:48)
s 3
∈S
In passing, we note that the result (15) can also be obtained by integrating out the bath orbitals
in a path integral, which can be done exactly if the bath orbitals appear only up to quadratic
order [44]. Furthermore, the right-hand sides of the last two lines in Eq. (15) are in line with
the Hermitian symmetry of the real-time quantities, i.e.,
≷ ≷ def. (cid:88) ≷
Σ (t,t) = [Σ (t,t)] = D (t)[g ((cid:15) ;t,t)] C (t) (16)
12σ (cid:48) − 21 (cid:48) ∗ − s∗σ (cid:48) 3sσ (cid:48) ∗ sσ
s ≷
(cid:88)∈S ≷
= C (t)g ((cid:15) ;t,t)D (t),
sσ 3sσ (cid:48) s∗σ (cid:48)
s ≷
∈S
≷ ≷
where we used g ((cid:15);t,t) = [g ((cid:15);t,t)] in the last equality.
(cid:48) (cid:48) ∗
−
For the discussion of the result (15) we assume that the initial state of the Hubbard dimer
is not correlated, i.e., ΣM = Σ = 0. To fix the bath parameters, it is instructive to begin with
(cid:100)
the decomposition of the offdiagonal components of the self-energy. This will yield the time-
dependent hopping parameters C (t) and D (t) and, subsequently, allows us to determine the
sσ sσ
parameters A (t) and B (t) from the first two lines in Eq. (15). As in Sec. 2.1, we choose the
sσ sσ
energies of the bath orbitals such that the initial occupations of the bath orbitals are either 0
(for the greater component of Σ) or 1 (for the lesser component of Σ), which leads to
≷ (cid:88)
iΣ (t,t) = C (t)D (t). (17)
± 12σ (cid:48) sσ s∗σ (cid:48)
s ≷
∈S
Asfunctionoftandt, thematricesiΣ> and iΣ< areneitherHermitiannorpositivedefinite
(cid:48) 12σ − 12σ
which generally rules out an eigenvalue and a Cholesky decomposition. This is in contrast to
the diagonal components of the self-energy discussed in Sec. 2.1 (Eqs. (7) and (8)). However,
one can perform a singular value decomposition of the objects as
≷ (cid:88)
iΣ (t,t) = U (t)S V (t), (18)
± 12σ (cid:48) sσ sσ sσ (cid:48)
≷
s
∈S3
where U and V denote unitary matrices and S is a real vector containing the non-negative
σ σ σ
singular values. As the singular values S typically decay rapidly in magnitude, we can further
sσ
use a low-rank decomposition by choosing the number of bath orbitals smaller than the number
of time steps. Incorporating the singular values symmetrically into the matrices U and V , we
σ σ
obtain the first set of bath parameters as
(cid:112) (cid:112)
C (t) = S U (t), D (t) = S V (t). (19)
sσ sσ sσ sσ sσ s∗σ
As a next step we can use the result of Eq. (19) to determine the bath parameters A (t)
sσ
and B (t). Taking into account the initial occupations, the first two lines of Eq. (15) can be
sσ
rewritten as (i = 1,2)
≷ ≷ ≷
∆ (t,t) = iΣ (t,t) Γ (t,t), (20)
iσ (cid:48) ± iiσ (cid:48) − iσ (cid:48)
with
≷ (cid:88) ≷ (cid:88)
∆ (t,t) = A (t)A (t), Γ (t,t) = C (t)C (t),
1σ (cid:48) sσ ∗sσ (cid:48) 1σ (cid:48) sσ s∗σ (cid:48)
≷ ≷
s s
∈S1 ∈S3
≷ (cid:88) ≷ (cid:88)
∆ (t,t) = B (t)B (t), Γ (t,t) = D (t)D (t). (21)
2σ (cid:48) sσ s∗σ (cid:48) 2σ (cid:48) sσ s∗σ (cid:48)
≷ ≷
s s
∈S2 ∈S3
≷ ≷
While the matrices iΣ and Γ are both Hermitian and positive definite by construction,
± iiσ iσ ≷
the positive definiteness does not hold in general for their difference ∆ , Eq. (20). For this
iσ
≷
reason, the partitioning of ∆ into the bath parameters A (t) and B (t), respectively, cannot
iσ sσ sσ
be performed through a Cholesky decomposition. Instead, we use an eigenvalue decomposition,
≷ (cid:88)
∆ (t,t) = P (t)E P (t), (22)
1σ (cid:48) sσ sσ s∗σ (cid:48)
≷
s
∈S1
≷ (cid:88)
∆ (t,t) = Q (t)F Q (t),
2σ (cid:48) sσ sσ ∗sσ (cid:48)
≷
s
∈S2
where the eigenvectors form the rows of the matrices P and Q and the corresponding
σ σ
eigenvalues are contained in the vectors E and F . Comparing Eqs. (21) and (22), we thus
σ σ
obtain the remaining bath parameters as
(cid:112) (cid:112)
A (t) = E P (t), B (t) = F Q (t). (23)
sσ sσ sσ sσ sσ sσ
2.2.2. Numerical validation The derivation presented in the previous section proofs that there
exists a procedure which maps the full Kadanoff-Baym equations of the Hubbard dimer onto
a noninteracting auxiliary model. In this section, we demonstrate the validity of the mapping
also numerically. To this end, we prepare the Hubbard dimer (at time t = 0) in the perfect N´eel
state |Ψ0(cid:105) = c†2 c†1 |0(cid:105) and use the second-order Born approximation of the self-energy,
↓ ↑
Σ2B,≷(t,t) = U2G≷ (t,t)G≷ (t,t)G≶ (t,t), (24)
ijσ (cid:48) jiσ (cid:48) jiσ¯ (cid:48) ijσ¯ (cid:48)
tomonitorthedecayofthelocalstaggeredmagnetizationm (t) = n (t) n (t) intheweak-
i i i
(cid:104) ↑ (cid:105)−(cid:104) ↓ (cid:105)
coupling regime. For convenience, we compute the Green’s function only for one spin direction.
The opposite component follows from the symmetry G (t,t) = G (t,t), where
ijσ¯ (cid:48) (L+1 i)(L+1 j)σ (cid:48)
L = 2 is the spatial dimension of the dimer. − −
Fig. 3a shows the self-consistent second Born result for an on-site interaction U = 0.5, where
the KBE has been solved with a time step of size δt = 0.01 (on a time window [0,6]) including
n = 20 orbitals for each bath i . The result is further compared to the Hartree dynamics
aux s
(given by the black dashed line) and to the exact dynamics (see the black solid line), which has
been obtained by exact diagonalization. On the short time scale, unlike the mean-field result,
correlations induce a damping of the oscillatory magnetization dynamics in the dimer. If we
U =0.5
1 0.2
Re
exact U =0.5 Im
H
0.5 2Bii 0.1
2Bij0
2Bij
) )
(t 0 (t 0
m0 X0↑
(a)
1
0.5 0.1 10 2
− − −
S
10 4 n↑ (b)
−
0 n 10
1 0.2
− 0 1 2 3 4 5 6 − 0 1 2 3 4 5 6
t [T 1] t [T 1]
− −
Figure 3. Hubbard dimer prepared in the N´eel state at t = 0. (a) Time evolution of the
staggered magnetization m (t) = n (t) n (t) on the left site of the dimer for U = 0.5 in
0 0 0
(cid:104) ↑ (cid:105)−(cid:104) ↓ (cid:105)
Hartree (H) and second-order Born approximation (2B), as computed from the self-consistent
auxiliary Green’s function Gaux(t,t). The line labeled 2Bij (2Bii) denotes the result for a non-
ijσ (cid:48)
local (local) self-energy. The curve denoted as 2Bij0 shows the result where Γ is additionally
iσ
set to zero in Eq. (20). (b) Time-dependent hopping parameters A (t) (red lines), B (t)
0σ 0σ
(orange lines), C (t) (blue lines) and D (t) (green lines). The inset shows the singular values
0σ 0σ
of the decomposition (18) of the non-local self-energy Σ (t,t).
12σ (cid:48)
consider nothing but the local parts of the self-energy, i.e., if we set to zero the bath parameters
C (t) and D (t) in the auxiliary system (Eq. (10)), the time evolution of the magnetization
sσ sσ
followsthebluesolidline(labeled2Bii). Theinclusionofthenon-localself-energyimprovesthis
result as expected, see the red dash-dotted line referred to as 2Bij. Particularly, we find good
agreement with the exact dynamics for times t (cid:46) 2.5. As additional result (denoted 2Bij0),
≷ ≷
we present the green dotted curve, where in Eq. (21) the quantities Γ and Γ have been
1σ 2σ
neglected in the determination of the bath parameters A (t) and B (t). We note that such an
sσ sσ
incomplete solution also leads to an acceptable dynamics, although there are discrepancies at
short times. Moreover, we emphasize that all schemes conserve the particle number as function
of time.
To further illustrate the method, Fig. 3b shows the time dependence of the complex hopping
parameters A (t) to D (t), i.e., for the first bath orbital s = 0 in each set <. Starting
0σ 0σ ≷ Si
from zero (at time t = 0 where the self-energy Σ vanishes), the time evolution of the bath
ijσ
parameters is rather complicated despite the simplicity of the Hubbard dimer. Moreover, the
inset in Fig. 3b displays the singular values S of Eq. (18) which decay exponentially in
n
↑
magnitude and thus allow one to use a low-rank representation of the bath 3 . Concerning
s
the baths 1 and 2 , we have included the 20 most important eigenvalues of each self-
s s
energy decomposition, i.e., those which have the largest absolute values. Note that an exact
representation of the self-energy using a time step of δ = 0.01 would require the dimension of
each individual bath to be n = 600 on the considered time interval.
aux
2.3. Generalization to an arbitrary lattice system
For a lattice system with more than two sites, the self-energy should have a similar auxiliary
representation as discussed above in Sec. 2.2.1. Without giving the proof, we state that the
i j
s s
(i j)
↔ s
noninteracting
Jsiσ Ksi→σj Ksjσ→i Jsjσ
... ... T ... ...
i j
Figure 4. Generalization of the auxiliary Hamiltonian approach to an extended lattice system
≷
with a hopping matrix T and an arbitrary non-local self-energy Σ (t,t). Together with the
ijσ (cid:48)
dashed lines, the open circles show the layout of the bath and how it is connected to the lattice
sites (filled circles) by complex hopping matrix elements Ksi→σ j(t), Ksjσ→i(t), Jsiσ(t) and Jsjσ(t).
offdiagonal components of the self-energy can always be expressed in the form (i = j)
(cid:54)
iΣ≷ = (cid:88) Ki j(t)[Kj i(t)] , (25)
± ijσ s→σ sσ→ (cid:48) ∗
≷
s
∈Si↔j
i j
withtime-dependenthoppingparametersKs→σ (t)connectingalatticesiteiandabathorbitals.
i j
In the course of this, it is important to remark that the parameters Ks→σ (t) are different from
j i
Ksσ→ (t). However, both sets follow from a single singular value decomposition of the left-hand
side of Eq. (25). On the other hand, the auxiliary representation for the diagonal components
of the self-energy is given by
Σ≷ = (cid:88) Ji (t)[Ji (t)] +(cid:88) (cid:88) Ki j(t)[Ki j(t)] , (26)
± iiσ sσ sσ (cid:48) ∗ s→σ s→σ (cid:48) ∗
≷ j=i ≷
s∈Si (cid:54) s∈Si↔j
where the bath parameters Ji (t) J (t) are defined as in Sec. 2.1 and can be obtained from
sσ ≡ isσ
an eigenvalue decomposition applied to Eq. (26). For an illustration of the resulting auxiliary
system see Fig. 4.
In summary, the one-particle auxiliary Hamiltonian which replaces the full Kadanoff-Baym
equation of Eq. (2) for an arbitrary non-local self-energy is of the form
 T J K 
σ σ
haσux(t) =  J†σ 0 0  , (27)
K 0 0
†σ
where the matrix elements of T are the hopping parameters T of the original Hubbard
ij
model (1), the matrix J is diagonal involving the bath parameters Ji (t), and the matrix
sσ
K includes the bath parameters Ksi→σ j(t) and Ksjσ→i(t). If we use naux bath orbitals for the
representation of each component Σ (with i j) of the self-energy and consider a Hubbard
ijσ
≤
model with L sites, the dimension of the auxiliary Hamiltonian haux is given by
σ
= L+2n L+n (L2 L) = (n +1)L+L2, (28)
aux aux aux aux
D −
i.e., the size of the generalized auxiliary system scales linearly with the size of the bath and
quadratically with the number of lattice sites. Moreover, we emphasize that the low-rank
representationofEqs.(25)and(26)[withfixedn n ]genericallyleadstoatimepropagation
aux t
(cid:28)
scheme which scales quadratically with the number of time steps n . This is one of the main
t
advantages of the method and allows the scheme to outperform standard KBE solvers.