Quantum transport simulation scheme including strong correlations and its application to organic radicals adsorbed on gold

2.1 MB·

by  Andrea Droghetti

Quantum transport simulation scheme including strong correlations and its
application to organic radicals adsorbed on gold
Andrea Droghetti1,∗ and Ivan Rungger2,†
1Nano-Bio Spectroscopy Group and European Theoretical Spectroscopy Facility (ETSF),
Universidad del Pais Vasco CFM CSIC-UPV/EHU-MPC and DIPC,Av.Tolosa 72 ,20018 San Sebastian, Spain
2National Physical Laboratory, Hampton Road, TW11 0LW, United Kingdom
We present a computational method to quantitatively describe the linear-response conductance
ofnanoscaledevicesintheKondoregime. Thismethodreliesonaprojectionschemetoextractan
Andersonimpuritymodelfromtheresultsofdensityfunctionaltheoryandnon-equilibriumGreen’s
7 functions calculations. The Anderson impurity model is then solved by continuous time quantum
1 Monte Carlo. The developed formalism allows us to separate the different contributions to the
0 transport, including coherent or non-coherent transport channels, and also the quantum interfer-
2 encebetweenimpurityandbackgroundtransmission. Weapplythemethodtoascanningtunneling

n microscope setup for the 1,3,5-triphenyl-6-oxoverdazyl (TOV) stable radical molecule adsorbed on
a gold. The TOV molecule has one unpaired electron, which when brought in contact with metal
J electrodes behaves like a prototypical single Anderson impurity. We evaluate the Kondo tempera-
9 ture,thefinitetemperaturespectralfunctionandtransportproperties,findinggoodagreementwith
2 published experimental results.

]
ll I. INTRODUCTION bridization and of its dynamic character55. Finally, the
a device transmission function and conductance are calcu-
h
lated within the Landauer framework46, by considering
- In recent years, much research effort has been ded-
s the Kohn-Sham (KS) eigenvalues as the single-particle
e icated to study the electronic transport through mag-
excitations. Although this assumption is formally not
m netic molecules and single atoms in order to combine
correct, since even in exact DFT only the highest occu-
. molecular electronics with spintronics1–3. Experiments pied molecular orbital (HOMO) can be rigorously asso-
t
a and theoretical works have demonstrated that molecular ciated to the negative of the ionization potential56–60,
m and atomic spin states can be inferred and sometimes it practically works well for metallic point-contacts61,
- switched through an electrical current4–18. Among the nanowires and nanotubes62,63, quasi two-dimensional
d many interesting phenomena arising in devices compris- systems64, and tunnel junctions65,66. In contrast, the
n ing magnetic molecules, there is the Kondo effect19–22.
o calculation of transport properties via the KS eigen-
Below the so called Kondo temperature, θ , character-
c K values encounters some drastic limitations in case of
[ istic of each system, the coupling between the electrons molecules. The fundamental gap between the HOMO
fromtheelectrodesandthespinofthemoleculepromotes and the lowest unoccupied molecular orbital (LUMO)
1
the formation of a many-body state with a fully or par-
v is often severely underestimated by the KS gap com-
tially quenched magnetic moment. This results in a new
5 puted with standard (semi-)local exchange-correlation
resonant transport channel at the electrodes Fermi level.
0 density functionals. Furthermore, non-local correlation
4 To date the Kondo effect has been studied in a number effects, such as the dynamical response of the electronic
8 of molecular devices8,9,14,17,23–44, and exotic manifesta- system to the addition of an electron or hole67–71, are
0 tions, such as the orbital Kondo effect, have also been
not captured. These shortcomings hinder the ability of
1. reported45. DFT+NEGF to predict the correct energy level align-
0
Ab-initio computational studies play a prominent role mentbetweenamoleculeandtheelectrodes,whichoften
7
in molecular electronics and spintronics. In particu- resultsinoverestimatedvaluesfortheconductance. Con-
1
: lar density functional theory (DFT) combined with the sequently, different improvements have been proposed,
v non-equilibrium Green’s functions (NEGF) formalism46, such as corrections for self-interaction error72,73, scissor
Xi known as DFT+NEGF for short, has become the dom- operator schemes74–80 and constrained-DFT80–84. Fur-
inant method to address electronic transport47–54. As thermore, in the recent years, there have been several
r
a initial step in a typical DFT+NEGF study of a molecu- attempts to move beyond the DFT+NEGF method by
lar device, one optimizes the atomic configuration of the using the GW approximation of the many-body pertur-
deviceactiveregion,usuallycalled“scatteringregion”or bation theory85–88.
“extended molecule”. Then the scattering region is joint The description of the Kondo effect, even at the quali-
to two semi-infinite left- and right-hand side leads (elec- tativelevel,stillrepresentsachallenge. Ontheonehand,
trodes),whoseeffectonthestatesinthescatteringregion electron correlations leading to the Kondo effect are be-
istakenintoaccountbytheso-calledleads’self-energies. yond the GW perturbative scheme. On the other hand,
This approach allows to treat the extended molecule as theDFTKSspectrumfailstodisplayanyKondo-related
an open system, with the leads’ self-energies that pro- feature, and so does the conductance computed via the
vide quantitative estimates of the molecule-electrode hy- Landauer approach. We note that in principle DFT is
2
able to capture the Kondo effect in one-orbital lattice
models89–91 if the exchange-correlation potential has the
correctderivativediscontinuityatintegernumberofelec-
trons, and if the conductance is computed from the den-
sity through the Friedel sum rule92,93 and not from the
KS states.
In order to describe the Kondo effect in real molecular
systemsandtoovercomethelimitationsofDFT+NEGF
and GW, recent studies have proposed to combine DFT
with model calculations, thus extending to molecular FIG. 1. Schematic of the 1,3,5-triphenyl-6-oxoverdazyl
electronics theoretical schemes originally proposed for (TOV) molecule.
the study of strongly correlated solid state materials.
These schemes include the dynamical mean-field the-
ory (DMFT)94. The combination of DFT and models calculationoftheelectron-electroninteractionenergyfor
is typically achieved by partitioning the system of in- the Anderson impurity. Here we suggest that a partially
terest in two coupled subsystems, a weakly correlated screened value should be used in order to reproduce the
one, whose electronic structure is well accounted for by experimental Kondo temperature.
DFT, and a strongly correlated one. Mathematically, The article is separated in two main parts. In the
thismeansthatpartofthesystemofinterestisprojected first part (Sec. II) we describe the theoretical methods,
onto the correlated sub-system, and that the rest is in- while in the second (Secs. III and IV) we present the
tegrated out as an effective bath. In case of molecular results for the TOV molecule on Au. In the first part
devices, this approach ultimately leads to the reformula- weinitiallysummarizetheprojectionschemeandoutline
tion of the electronic structure and transport problem in how the different contributions to the transmission are
terms of an effective Anderson impurity model (AIM), computed (Secs. IIA to IIE), and then we describe the
which then has to be solved either exactly or within employedCTQMCalgorithm(Sec. IIF)andtheanalytic
someapproximations. Thepotentialofthisapproachhas continuation of the CTQMC data (Sec. IIG). In the
beenfirstlydemonstratedbycomparingtophotoemission secondpart,westartbylistingthecomputationaldetails
experiments95 the computed spectral properties of sin- oftheDFTandNEGFcalculations(Sec. III).Wepresent
gle magnetic atoms96–99 and molecules100,101 on metal- the DFT results for the gas phase TOV, for the TOV
lic surfaces. Then, the linear-response (i.e. zero-bias) on gold (Sec. IVA) and the DFT+NEGF results for
transportpropertieshavebeenaddressedforexampleby the transport properties (Sec. IIE). We then introduce
Smogunov,Tosattiandco-workers102–105,andJacoband the AIM and discuss its solution within the mean-field
co-workers106–111. We have recently contributed to fur- approximation (Sec. IVC) and by CTQMC (Sec. IVD).
ther develop the DFT+NEGF scheme including DMFT Finally, we present the transport properties computed
tostudyspintransportinsolidstatedevicessuchasmul- viaDFT+NEGF+CTQMC(Sec. IVE)andweconclude
tilayered heterostructures112. (Sec. V).
In this article, we present our scheme to project out
from DFT+NEGF an AIM, which is then solved nu-
merically using continuous time quantum Monte Carlo II. METHOD AND IMPLEMENTATION
(CTQMC)113. Thedevelopedmethodallowsustoevalu-
atethetemperature-dependentlinear-responsetransport The first step of the method is the separation of the
properties of magnetic molecules on metal surfaces. An Anderson impurity sub-system from the full system cal-
important class of such molecules is formed by the sta- culated in the DFT+NEGF setup. This is done by an
ble organic radicals, which are are paramagnetic com- appropriate projection scheme, outlined in subsections
pounds presenting an unpaired electron in a singly oc- IIA, IIB, and IIC, and schematically summarized in
cupied molecular orbital (SOMO)114–122. In particular, Appendix A6. Then, the AIM with an effective inter-
hereweconsiderthe1,3,5-triphenyl-6-oxoverdazyl(TOV, action term is introduced and solved. This means that
Fig. 1). This is the first organic radical for which the themany-bodyGreen’sfunctionandself-energyarecom-
Kondo effect was experimentally observed in a scanning puted (subsection IID) so that the transport properties
tunnelingmicroscope(STM)setupwithagoldsubstrate, can be obtained (subsection IIE). The method can, in
andtheKondotemperaturewasreportedtobeabout37 principle, address both zero- and finite-bias transport
K118. SincethesystemiswellcharacterizedbySTM,the problems provided the availability of a computationally
comparisonofthecalculateddensityofstatesandKondo efficientout-of-equilibriumsolverfortheAIM.Neverthe-
temperature with the experiment serves as a stringent less, in this work we only consider the zero-bias case,
test for the theory. Additionally, we demonstrate how and therefore the solution of the AIM is achieved by us-
thedifferentcontributionstotheconductancethatorigi- ingCTQMCforquantumsystemsinequilibriumatfinite
natefromelastic,non-coherentandquantuminterference temperatures, as outlined in subsection IIF. The use of
effects can be disentangled. An open issue concerns the CTQMCrequiresaschemetocarryouttheanalyticcon-
3
The setup used in the calculations is schematically il-
DFT: Relaxed atomic structure
lustratedinFig. 3,wherethebasissetdescribedaboveis
usedtoexpandthechargedensity52,124. Likeinstandard
DFT+NEGF: Extended molecule Green’s function electron transport simulations, the system is first split
𝐺(𝐸)= 𝐸𝑆−𝐻−Σ 𝐸 −Σ 𝐸 −1 intoasemi-infiniteleftlead, aso-calledscatteringregion
L R
(or extended molecule (EM)), and the semi-infinite right
lead52. We denote the number of basis orbitals within
• Define the Anderson impurity (AI) within the extended molecule by a
the EM by N. We can then introduce the Hamiltonian
set of interacting orbitals: ψ𝑖 matrixoftheEMas(H) =(cid:104)φ |Hˆ|φ (cid:105),wherethebasis
• Construct the projection matrix to decouple AI from bath (B): 𝑊(ψ,𝑆,𝐻) µν µ ν
𝑖 orbitals span all orbitals within the EM, and where Hˆ is
𝐻 =𝑊†𝐻𝑊= 𝐻 AI 𝐻 AI,B 𝐺 =𝑊−1𝐺𝑊−1†= 𝐺 AI 𝐺 AI,B the Hamiltonian operator. Since in general the basis set
𝐻 𝐻 𝐺 𝐺
B,AI B B,AI B is non-orthogonal, we also need to introduce the overlap
• Evaluate AI hybridization function, Δ AI, so that matrix of the EM, S, given by
𝐺 (𝑧)= 𝑧−𝐻 −Δ (𝑧)−1
AI AI AI
(S) =(cid:104)φ |φ (cid:105). (2)
µν µ ν
Δ AI(𝑧), 𝐻 AI We note here that we shift all energies of leads and EM
in such a way to set the Fermi energy, E , equal to 0.
Solve AI model F
Such a global shift of the spectrum does not affect the
CTQMC + analytic continuation
properties of the system.
Σ MB(𝐸)
AI The EM is then further subdivided in a total of 3
subsystems. A set of pre-determined wave-functions ψ
i
𝐺 MB(𝐸)= 𝐸𝑆 −𝐻 −Σ L 𝐸 −Σ R 𝐸 −Σ MB(𝐸)−1 (i ∈ [1,N ]), defines the Anderson impurity (AI), with
AI
𝑇=Tr Γ 𝐺 MB†Γ 𝐺 MB =𝑇 +𝑇 +𝑇 NAI being to the number of interacting states. The in-
L R B AI I teracting region (IR) includes all basis orbitals inside the
EM that contribute to the AI, which corresponds to the
FIG. 2. Schematic representation of the method for the cal- collection of those basis orbitals of the EM where any
culationofthetransportpropertiesinpresenceofinteracting of the {ψi} is non-zero. The number of basis orbitals in
Anderson impurities. the IR, N , is usually much larger than N . We re-
IR AI
quire that the overlap and Hamiltonian matrix elements
between orbitals within the IR and the orbitals of the
tinuation of the many-body Green’s function from the leads outside the EM is zero, so that the EM has to be
discrete imaginary Matsubara frequencies onto real en- chosen large enough to ensure this. As last subspace we
ergies This is presented in subsection IIG. A schematic introducetheextendedinteractingregion(ER),whichin-
summary of the whole method is shown in Fig. 2. cludes all basis orbitals inside the EM that have a finite
overlaporHamiltonianmatrixelementwiththebasisor-
bitals within the IR. We denote the number of orbitals
A. NEGF transport setup
within the ER as N , with N ≥ N . Ordered by
ER ER IR
decreasing size, the 4 subsystems are then: EM, ER, IR,
The method has been implemented in the NEGF code and AI. The AI is a subsystem of the IR, the IR is a
Smeagol52,whichusesalinearcombinationofatomicor- subsystem of the ER, and the ER is a subsystem of the
bitals(LCAO)basisset{φµ},andobtainstheKSHamil- EM (see Fig. 3).
tonian from the DFT package Siesta124. Note however The first step is therefore to define the set of molecu-
that the method is general and can be readily used for lar orbitals, {ψ }, inside the EM that constitute the AI.
i
any code based on the LCAO approach. Each basis or- The method outlined here is applicable for any arbitrary
bital|φµ(cid:105)inSmeagol ischaracterizedbyitsintegerindex set of wave-functions, although in general the choice of
µ, which is a collective index that includes the atom, I, the {ψ } is based on physical intuition. In the simplest
i
the orbital n, and the angular momentum (l,m) indices. caseonecansimplytakeacombinationofdorf orbitals
Theorbitalindexncanrunoverdifferentradialfunctions for correlated magnetic atoms in the system, or else the
corresponding to the same angular momentum, accord- wave-functionsoftheHOMOorLUMOofamoleculeat-
ingtothemultiple-zetasscheme123. AnyoperatorOˆ can tached to metal electrodes. A practical way to construct
be expressed in this basis by using its matrix form O, such a set of interacting molecular orbitals is to extract
with the matrix elements given by fromthefullEMHamiltonian,H,andoverlap,S,matri-
(O) =(cid:104)φ |Oˆ|φ (cid:105). (1) cessub-blockswithbasisorbitalsonthemolecule, which
µν µ ν we denote as H and S , respectively. We can then
M M
Asamatterofnotationweremarkthatingeneralupper calculate the eigenvalues (cid:15) and eigenvectors ψ of these
i i
case symbols represent matrices, and to distinguish op- sub-systems by solving H ψ =(cid:15) S ψ . Note that each
M i i M i
erators from matrices we explicitly add a hat on top of of the ψ is a vector of dimension N .
i IR
the symbol for operators. Once the set of wave-functions that defines the AI as
4
Left lead Extended molecule (EM) Right lead within the EM on the left (right) of the ER, so that
Extended interacting region (ER) N =N +N +N . Thedimensionsoftheoff-diagonal
α β ER
Interacting region (IR) Anderson Impurity (AI) matrixblocksaredeterminedfromtheonesofthediago-
nalblocks,andarethereforenotexplicitlygiven. Wesub-
dividetheHamiltonianmatrixinananalogousway. Note
thatbyconstructionthefollowingimportantrelationsare
satisfied: S U =0, S U =0, H U =0,
α,ER ER β,ER ER α,ER ER
and H U =0.
H H H β,ER ER
L H H R The Green’s function (GF) matrix of the EM is then
LM RM
given by the standard form52
FIG. 3. (Color online) Schematic representation of the two-
terminaldevicecomprisingtheleftlead,theextendedmolecule G(z)=[zS−H −Σ (z)−Σ (z)]−1, (6)
L R
(EM), and the right lead. We denote the Hamiltonian ma-
trix of the semi-infinite left (right) lead by HL (HR), the which has the block-matrix structure analogous to the
one of the EM as H, and the coupling Hamiltonian matrix
one of the H and S matrices in Eq. (5); z is an arbitrary
as H (H ). The EM is further subdivided in the inter-
LM RM complex number, and Σ (z) and Σ (z) are the left and
acting region (IR) and the extended interacting region (ER). L R
rightleads’self-energiesthatdescribethecouplingofthe
The IR includes a set of interacting molecular orbitals {ψ },
i
EM to the leads. These are computed according the al-
with i∈[1,N ] (N is the number of interacting molecular
AI AI
orbitals;intheschematicfigureN =4),formingtheAnder- gorithm in Ref. 53. Importantly, G(z) can be either the
AI
son impurity (AI). The ER consists of all orbitals in the EM retardedGFifz =E+iδ,withE therealenergyandδ a
that have finite overlap with the IR orbitals, so that the IR vanishinglysmallpositiverealnumber,ortheMatsubara
orbitalsarealwaysalsopartoftheER.Intotalthesystemis GF if z = iω , with ω = (2n+1)π/β, for n ∈ Z and
n n
thereforesubdividedinto4subspacesofdecreasingsize: EM, β = 1/kθ the inverse temperature (k is the Boltzmann
ER, IR, and AI. constant, θ is the temperature). Spin indices are omit-
ted in above equations and in the following subsections
to simplify the notation, but they will be explicitly re-
{ψ ,ψ ,...,ψ }ischosen(herewehaveadded
1,IR 2,IR NAI,IR introduced when required to emphasize spin-dependent
the subscript ”IR” to indicate explicitly that the wave-
relations.
functionsextendovertheIR),theprojectionmatrixonto
the AI inside the IR, UIR, is then defined as
U =(cid:0)ψ ψ ψ ... ψ (cid:1), (3) B. Projection to the Anderson impurity subsystem
IR 1,IR 2,IR 3,IR NAI,IR
which is of dimension N × N . We then construct
IR AI Here we outline how to separate explicitly the AI sub-
the ER by adding to the IR all basis orbitals within the
system from the rest of the system, which we refer to as
EM that have finite overlap with the IR basis orbitals,
the “bath”, and that includes the orthogonal subspace
and construct the set of wave-functions {ψ }. Each
i,ER totheAIwithintheEMaswellasthesemi-infiniteelec-
ψ is a vector of length N , and is equal to ψIR for
i,ER ER i trodes. To this aim we introduce a basis transformation
its elements within the IR, and 0 for the others. These
matrix, W, which transforms the overlap, Hamiltonian,
vectorsthereforedefinetheAIintheER,andleadtothe
and self-energy matrices as
ER projection matrix onto the AI, U , given by
ER
(cid:0) (cid:1) S¯=W†SW, H¯ =W†HW, (7)
U = ψ ψ ψ ... ψ , (4)
ER 1,ER 2,ER 3,ER NAI,ER Σ¯ (z)=W†Σ (z)W, Σ¯ (z)=W†Σ (z)W, (8)
L L R R
which is of dimension N ×N . Note that the simula-
ER AI
tionsetupinSmeagol requiresthatoneleads’unitcellis where we denote the transformed matrices with a bar
included at both the left and right ends of the EM. For on top of the symbol. This transformation is required
the projection onto the AI we then further require that to bring any general matrix extending over the orbitals
the orbitals of those cells cannot be part of the IR, while of the EM, M, into a transformed form, M¯, in which
they are allowed to be part of the ER. the top left corner describes the AI, the bottom right
If the basis orbitals indices are approximately ordered corner describes the part of the bath included in the EM
from left to right in the EM, then the overlap matrix (B), and the off-diagonal blocks describe the connection
of the EM52,124 can be written in the following general terms. We note that the full bath includes orbitals from
block-matrix structure both the EM and the semi-infinite electrodes, while the
N ×N matrix M¯ contains only the N = N −N
 S S S  B B B B AI
αα α,ER αβ bath orbitals within the EM. In mathematical terms we
S =Sα†,ER SER Sβ†,ER , (5) therefore require M¯ to have the following block-matrix
S† S S structure:
αβ β,ER ββ
where S is a N × N matrix. S (S ) is a (cid:18) M¯ M¯ (cid:19)
ER ER ER αα ββ M¯ = AI AI,B . (9)
square matrix that includes all the N (N ) orbitals M¯ M¯
α β B,AI B
5
We also require the projection to lead to a zero overlap Thegeneralstructureoftheresultingmatriceshasthe
between the AI orbitals and the bath orbitals. required shape given in Eq. (9), with
The form of W to achieve this for the transport setup
is derived in Appendix A1, and is given by  Sαα S¯α,NI Sαβ 
 0 1 0 0  S¯B =S¯α†,NI S¯NI S¯β†,NI , (22)
Nα S† S¯ S
W =WAI 0 WNI 0 . (10) αβ β,NI ββ
0 0 0 1 S¯ =(cid:0)0 0 0(cid:1), (23)
Nβ AI,B
Here and in the following we denote an identity matrix
and
of dimensions m×m as 1 . We have introduced
m
 H H¯ H 
W =U W , (11) αα α,NI αβ
AI ER 2,AI H¯B =H¯α†,NI H¯NI H¯β†,NI , (24)
which is equal to the projection matrix UER [Eq. (4)], H† H¯ H
multipliedbytheN ×N matrixW ,whichingen- αβ β,NI ββ
eral is constructed iAnI suchAIa way to o2r,tAhIogonalize both H¯AI,B =(cid:0)0 H¯AI,NI 0(cid:1). (25)
the overlap and Hamiltonian matrices of the AI itself.
Note that the transformed self-energy matrices extend
In principle this second transformation is arbitrary and
only over the block of the bath inside the EM, and are
can also be omitted, depending on how the AI problem
zero for the other matrix blocks. This is ensured by the
is solved. The N ×N matrix W spans the orthog-
ER NI NI
onalspacetotheAIinsidetheER,sothatW† W =0, requirement that the orbitals of the left and right leads’
AI NI unit cells included at the boundaries of the EM are not
and therefore projects onto the non-interacting part of
part of the IR.
the ER. We denote the number of non-interacting states
The transformation for the GF is given by
within the ER as N , so that N =N −N . There
NI NI ER AI
issomefreedomintheconstructionofthematricesW
2,AI G¯(z)=W−1G(z)W−1†, (26)
and W , since they are not uniquely defined. We con-
NI
struct them in such a way to leave the non-interacting
and can also be evaluated directly in the transformed
part close to the original system, and the detailed rela-
system as
tions to construct W and W are given in Appendix
2,AI NI
A1, Eqs. (A14) and (A18). G¯(z)=(cid:2)zS¯−H¯ −Σ¯ (z)−Σ¯ (z)(cid:3)−1. (27)
The explicit form of the final transformed overlap ma- L R
trix is then evaluated to
It is possible to evaluate the required inverse of W by
1N0AI S0 S¯0 S0  blocks (Appendix A2), and the result is
S¯= 0 S¯α†α,NαI S¯αN,NII S¯β†α,NβI , (12)  0 WiAI 0 
0 Sα†β S¯β,NI Sββ W−1 =1Nα 0 0 , (28)
 0 W 0 
iNI
which as required has zero overlap between AI and bath 0 0 1
orbitals, and where Nβ
S¯ =W† S W , (13) where
NI NI ER NI
S¯α,NI =Sα,ERWNI, (14) WiAI =W2−,A1IWi,ER,ψSER, (29)
S¯ =S W . (15)
β,NI β,ER NI
with
The final general form of the projected Hamiltonian ma-
trix is (cid:16) (cid:17)−1
W = U †S U U †. (30)
i,ER,ψ ER ER ER ER
 (cid:15) 0 H¯ 0 
AI,D AI,NI
0 H H¯ H
H¯ =H¯A†I,NI H¯α†α,NαI H¯αN,NII H¯β†α,NβI , (16) TinhAepfopremndoixf tAhe2,NENqI.×(AN2E6R). block-matrix WiNI is given
0 H† H¯ H
αβ β,NI ββ
where
C. Hybridization function
(cid:15) =W† H W , (17)
AI,D AI ER AI
H¯ =W† H W , (18) By removing the AI orbitals from the system we can
AI,NI AI ER NI
introduce the GF of only the bath orbitals within the
H¯ =W† H W , (19)
NI NI ER NI EM, g¯(z), as
H¯ =H W , (20)
α,NI α,ER NI
H¯ =H W . (21) g¯(z)=(cid:2)zS¯ −H¯ −(cid:0)Σ¯ (cid:1) (z)−(cid:0)Σ¯ (cid:1) (z)(cid:3)−1, (31)
β,NI β,ER NI B B L B R B
6
which can be written in a block-matrix structure analo- oftheleadsisaccountedforinthebathGreen’sfunction
gous to the one of S¯ [Eq. (22)] as [Eq. (31)], and therefore in the hybridization function
B
[Eq.34]throughtheprojectedself-energies(cid:0)Σ¯ (cid:1) (z)and
  L B
g¯αα(z) g¯α,NI(z) g¯αβ(z) (cid:0)Σ¯ (cid:1) (z),sothatthefullH¯ neverneedstobeexplic-
g¯(z)=g¯NI,α(z) g¯NI(z) g¯NI,β(z). (32) itlyRcoBmputed. TB
g¯ (z) g¯ (z) g¯ (z)
βα β,NI ββ In order to account for many-body correlation effects,
we supplement the AI with an effective Coulomb inter-
By using also Eqs. (27), (12) and (16) the GF on the
action, expressed by the operator
AI is then obtained as
Hˆ =Hˆ −Hˆ . (41)
G¯ (z)=(cid:2)z−(cid:15) −∆¯ (z)(cid:3)−1. (33) I C dc
AI AI,D AI
Then the interacting AIM Hamiltonian operator is de-
Herewehaveintroducedtheso-calledhybridizationfunc-
fined as
tion, ∆¯ (z), which is given by
AI
Hˆ =Hˆ +Hˆ . (42)
IAIM AIM I
∆¯ (z)=H¯ g¯ (z)H¯† . (34)
AI AI,NI NI AI,NI Hˆ is typically chosen to have the form of a general-
C
We remark that in general ∆¯AI(z) is a dense NAI×NAI ized Hubbard-like interaction, while Hˆdc is the double-
matrix. Therefore, optionally, as alternative possibility counting correction94. This double-counting correction
one can modify the transformation matrix W to an is required in order to subtract the correlation effects in
2,AI
energy-dependent form that diagonalizes ∆¯ (z) rather theAIthatarealreadyincludedattheKSlevel. Theex-
AI
than H¯AI. In Appendix A3 we derive an equivalent, act form of Hˆdc is unfortunately not known, and several
commonlyusedexpressionforthehybridizationfunction. approximationshavebeenintroduced,themostcommon
Furthermore, in Appendix A4 we show that ∆¯ (z) ex- one being the so-called “fully localized limit”94. This is
AI
hibits the correct physical decay for large z. the one that we also employ.
Once G¯ (z) and g¯(z) are known, then the remaining For a single impurity one-orbital Anderson model
AI
block matrices of the GF within the EM [see Eq. (9) for (SIAM) one has
the general matrix structure] can be evaluated to
Hˆ =Unˆ nˆ , (43)
C ↑ ↓
G¯ (z)=g¯(z)+g¯(z)H¯† G¯ (z)H¯ g¯(z), (35)
B AI,B AI AI,B withtheHubbardU beingarealnumber,andnˆ (nˆ )the
↑ ↓
G¯ (z)=−G¯ (z)H¯ g¯(z), (36) occupationoperatorforup-spin(down-spin)electronson
AI,B AI AI,B
G¯ (z)=−g¯(z)H¯† G¯ (z). (37) theAI.Thedouble-countingcorrectioninthefullylocal-
B,AI AI,B AI ized limit has the simple expression
Hˆ =U(n− 1)(cid:88)nˆ , (44)
D. Effective Coulomb interaction and many-body dc 2 σ
self-energy σ
where n is the DFT total occupation of the impurity.
The AI is fully characterized at the KS-level through
the “on-site” upper-block (cid:15) of the Hamiltonian ma- ThesolutionoftheinteractingAIMleadstothemany-
AI,D
trix in Eq. (16), which is coupled to the bath via the body GF on the AI, G¯MB(z). One can then define the
AI
off-diagonal block H¯ . The AIM Hamiltonian opera- many-body self-energy on the AI, Σ¯MB(z), by means of
AI,NI AI
tor can be rewritten in its standard form as an operator the Dyson equation
in second quantization Hˆ = Hˆ¯ +Hˆ¯ +Hˆ¯ ,
AIM AI,D TB AI,NI Σ¯MB(z)=G¯−1(z)−(cid:0)G¯MB(cid:1)−1(z). (45)
with AI AI AI
Hˆ¯ =(cid:80)NAI ((cid:15) ) (cid:80) dˆ† dˆ , (38) Since we are considering the case of a single interacting
AI,D i,j=1 AI,D ij σ iσ jσ region, the full many-body self-energy matrix of the EM
Hˆ¯ =(cid:80)∞ (H¯ ) cˆ† cˆ , (39) is zero for all elements, except for those of the AI
TB p,q=1 TB pq pσ qσ
Hˆ¯AI,NI =(cid:80)Ni=A1I(cid:80)Np=N1I(cid:104)(H¯AI,NI)ipdˆ†iσcˆpσ+c.c(cid:105). (40) Σ¯MB(z)=(cid:18)Σ¯MAIB(z) 0 (cid:19). (46)
0 0
NB
Here dˆ(†) is the annihilation (creation) operator for an
iσ Here we use the notation 0 to denote a m×m matrix
m
electron of spin σ in the orbital i of the AI, while cˆ(p†σ) with all zeroes; the sizes of the off-diagonal blocks are
is the annihilation (creation) operator for an electron of determinedbytheonesofthediagonalblocks. IfΣ¯MB(z)
spin σ in an orbital p of the total bath, whose Hamilto- is known, it can be used to directly calculate the many-
nian matrix is formally written as H¯TB. We note that body GF as
H¯ includes theelementsof H¯ obtained fromthe pro-
TB B
jection of the EM Hamiltonian with the Eq. (7), as well G¯MB(z)=(cid:2)zS¯−H¯ −Σ¯ (z)−Σ¯ (z)−Σ¯MB(z)(cid:3)−1.
L R
astheelementsofthesemi-infiniteelectrodes. Theeffect (47)
7
Incaseswhenthemany-bodyself-energyisrequiredin and Γ¯ (E)=W†Γ (E)W. Note that if the Mat-
{L,R} {L,R}
the original basis, this can be obtained by applying the subara GF is computed by using Eqs. (6) or (47) with
inverse transformation as z = iω , one has to perform the analytic continuation
n
ΣMB(z)=W−1† Σ¯MB(z)W−1. (48) from the complex Matsubara energies to the real axis in
order to obtain the retarded GF prior to the calculation
By using Eqs. (28) and (46) the explicit structure of of the transmission. The procedure that we have chosen
ΣMB(z) becomes tocarryoutsuchoperationisdescribedinSec.IIG.Note
also that we can equally express the transmission and
 
0 0 0
Nα current in terms of the quantities in the original basis
ΣMB(z)= 0 ΣMERB(z) 0 . (49) and in the transformed one. In what follows we work in
0 0 0
Nβ the transformed basis.
The transport properties of a molecular device in
The non-zero block of ΣMB extends over the ER, and is
DFT+NEGFareusuallyexpressedintermsofthetrans-
given by
mission coefficient for the KS system52, denoted as
ΣMERB(z)=Wi†AIΣ¯MAIB(z)WiAI, (50) T0(E), which is given by52
with WiAI given in Eq. (29). After calculating ΣMERB(z) T0(E)=Tr(cid:2)Γ¯L(E)G¯†(E)Γ¯R(E)G¯(E)(cid:3). (57)
one can then obtain the full many-body GF in the origi-
nal basis as It has the analogous structure to T(E), but with the
many-body GF replaced by the KS one. Therefore T(E)
GMB(z)=(cid:2)zS−H −Σ (z)−Σ (z)−ΣMB(z)(cid:3)−1.
L R accounts for the renormalization of the coherent trans-
(51)
port properties via many-body effects not included at
the KS DFT level, and in the following, we refer to I
R,L
as the coherent, or elastic, component of the current. In
E. Current and transmission
contrast, I represents the incoherent component of
R,AI
the current.
The current flowing out of the EM into the right elec-
The incoherent component of the current flowing from
trode, I , can be written as sum of a component trans-
R the AI to the right electrode is given by (see Appendix
mitted into the right electrode from the left electrode,
A5)
I ,andacomponentflowingfromtheAIintotheright
eRle,cLtrode, IR,AI,125,126 I = e (cid:90) dE Tr(cid:2)(cid:0)F¯MB(E)−f (E)(cid:1)
R,AI h R
I =I +I . (52)
R R,L R,AI (cid:16)Γ¯MB(E)G¯MB†(E)Γ¯ (E)G¯MB(E)(cid:17)(cid:105), (58)
A similar expression holds for the current flowing from R
the left electrode into the EM, I , as outlined in Ap-
L where we have introduced the occupation matrix of the
pendix A5. At steady state the current conservation AI, F¯ (E),126 and
MB
condition implies that I =I .
L R
The value of IR,L can be evaluated using the left to Γ¯MB(E)=i(cid:104)Σ¯MB(E)−(cid:0)Σ¯MB(E)(cid:1)†(cid:105), (59)
right energy dependent many-body transmission coeffi-
cient, T (E), which we will simply denote as transmis-
R,L whichisdefinedinasimilarwaytothecouplingmatrices
sion,T(E),fromnowoninordertosimplifythenotation.
inEqs.(55)and(56),butwiththemany-bodyself-energy
I is then given by
R,L replacing the the leads’ self-energies. The general struc-
e (cid:90) ture of the occupation matrix is analogous to the one of
IR,L = h dE (fL(E)−fR(E))T(E), (53) Σ¯MB in Eq. (46), and can be written as
where e is the electron charge, h is the Planck constant, (cid:18)F¯MB(E) 0 (cid:19)
f (f ) is the Fermi Dirac distribution at the chemical F¯MB(E)= AI . (60)
L R 0 0
potential of the left (right) electrode, and NB
T(E)=Tr(cid:104)Γ (E)GMB†(E)Γ (E)GMB(E)(cid:105) The quantity inside the trace of Eq. (58) can
L R be interpreted as a transmission matrix, given
=Tr(cid:2)Γ¯ (E)G¯MB†(E)Γ¯ (E)G¯MB(E)(cid:3). (54) by (cid:16)Γ¯MB(E)G¯MB†(E)Γ¯ (E)G¯MB(E)(cid:17), times the
L R R
T includes all interference effects between the possible difference in the distribution matrices, given by
transport channels across the system, including the ef- (cid:0)F¯MB(E)−fR(E)(cid:1). It therefore has the analogous
fects of interactions on the AI. The so called coupling structure to IR,L in Eq. (53), with the only difference
matrices Γ (E) are defined as that since F¯MB is a matrix it cannot be moved outside
L,R
(cid:104) (cid:105) the trace. The analogous equations for the current from
ΓL(E)=i ΣL(E)−Σ†L(E) , (55) the left lead into the EM, IL, are given in Appendix A5.
ΓR(E)=i(cid:104)ΣR(E)−Σ†R(E)(cid:105), (56) orImf woreeagsesnuemraelltyhtahtatthtehmeiartvraicluesesHaαreβsamndallSeαnβouargehztehraot,
8
they can be neglected, and that Γ (E) (Γ (E)) extends which depends only on the GF of the AI, and on the
L R
onlyovertheregionα(β), thenwecanobtainthetrans- hybridization matrices of the AI, given by
mission function from the GF block matrix of the NI as
γ¯ (E)=H¯ g¯ (E)Γ¯ (E)g¯† (E)H¯† ,(71)
(cid:104) (cid:105) L,AI AI,NI NI L,NI NI AI,NI
T(E)=Tr Γ¯L,NI(E)G¯MNIB†(E)Γ¯R,NI(E)G¯MNIB(E) , (61) γ¯R,AI(E)=H¯AI,NI g¯N†I(E)Γ¯R,NI(E)g¯NI(E)H¯A†I,NI.(72)
(cid:104) (cid:105) We can perform the analogous transformations also for
where Γ¯L,NI(E) = i Σ¯L,NI(E)−Σ¯†L,NI(E) and IR,AI from Eq. (58), and obtain
(cid:104) (cid:105)
Γ¯R,NI(E) = i Σ¯R,NI(E)−Σ¯†R,NI(E) . The matri- I = e (cid:90) dE Tr(cid:2)(cid:0)F¯MB(E)−f (E)(cid:1)
ces Σ¯ (E) are given by R,AI h AI R
{L,R},NI
Σ¯ (E)=K¯ (cid:2)K¯ −(cid:0)Σ¯ (cid:1) (E)(cid:3)−1K¯ ,(62) Γ¯MAIB(E)G¯MAIB†(E)γ¯R,AI(E)G¯MAIB(E)(cid:105), (73)
L,NI NI,α αα L αα α,NI
Σ¯ (E)=K¯ (cid:104)K¯ −(cid:0)Σ¯ (cid:1) (E)(cid:105)−1K¯ ,(63) where we have defined
R,NI NI,β ββ R ββ β,NI Γ¯MB(E)=i(cid:104)Σ¯MB(E)−Σ¯MB†(E)(cid:105). (74)
AI AI AI
with K¯ = ES¯−H¯ having the analogous block-matrix
structuretoH¯ inEq. (16),andforrealenergiesK¯ = In general therefore the calculation of the current re-
NI,β
K¯† and K¯ = K¯† . Note that K¯ = K , quires the knowledge of the impurity occupation matrix
K¯β,N=I K ,aNnId,αwiththαe,aNsIsumptionsmadeαinαthispaαrαa- FMB(E), which in turn requires the solution of the non-
grβaβph weβhβave (cid:0)Σ¯ (cid:1) =(Σ ) and (cid:0)Σ¯ (cid:1) =(Σ ) . equilibrium problem. In Refs.127,128 a specific shape is
L αα L αα R ββ R ββ implicitly assumed for this matrix126, which allows sim-
With Eq. (35) the GF block matrix of the NI can be
plifiedestimatesofthecurrents. Theapproachproposed
obtained from the one of the AI as
in those references requires the inversion of the coupling
G¯MB(E)=g¯ (E)+∆G¯MB(E), (64) matrices ΓL and ΓR, which is however not defined in the
NI NI NI general case53.
We now consider the special case where γ¯ (E) =
with L,AI
λγ¯ (E), with λ a constant125. This rather strict con-
R,AI
∆G¯MB(E)=g¯ (E)H¯† G¯MB(E)H¯ g¯ (E). (65) dition can usually only be fulfilled for a SIAM, so that
NI NI AI,NI AI AI,NI NI
in the remainder of this section we consider only the
The transmission can then be decomposed in three com- SIAM. In this case it is possible to avoid the calculation
ponents of F¯MB(E) by using the current conservation condition
I =I 125, and one obtains
L R
T(E)=TB(E)+TAI(E)+TI(E), (66) e (cid:90)
I = dE (f (E)−f (E))T (E), (75)
wherewehaveintroducedthebackgroundorbathtrans- R,AI h L R R,AI
mission with
(cid:104) (cid:105)
TB(E)=Tr Γ¯L,NI(E)g¯N†I(E)Γ¯R,NI(E)g¯NI(E) , (67) TR,AI =
Γ¯MB(E)G¯MB†(E) γ¯L,AI(E)γ¯R,AI(E) G¯MB(E),(76)
the transmission through the AI AI AI γ¯ (E)+γ¯ (E) AI
L,AI R,AI
(cid:104) (cid:105) wherewehaveusedthefactthatfortheSIAMallquanti-
T (E)=Tr Γ¯ (E)∆G¯MB†(E)Γ¯ (E)∆G¯MB(E) ,
AI L,NI NI R,NI NI tiesintheequationarejustnumbersinsteadofmatrices.
(68) Using Eqs. (53) and (75) for the special case of a
and the interference term of the transmission SIAM,andγ¯ (E)=λγ¯ (E),wecanwritethetotal
L,AI R,AI
(cid:104) (cid:105) current as
T (E ) =Tr Γ¯ (E)∆G¯MB†(E)Γ¯ (E)g¯ (E)
I L,NI NI R,NI NI e (cid:90)
+Tr(cid:104)Γ¯ (E)g¯† (E)Γ¯ (E)∆G¯MB(E)(cid:105). (69) IR = h dE (fL(E)−fR(E))Tt(E), (77)
L,NI NI R,NI NI
withthetotaleffectivetransmission,whichincludeselas-
Often the background and AI transmission do not inter- tic and incoherent terms, given by
fere significantly, so that it is useful to give estimates for
theirseparatevalues. Inthatcasetypicallythetransmis- Tt(E)=TB(E)+TAI(E)+TI(E)+TR,AI(E). (78)
sion is composed of a background value, onto which the
We can collect the terms that describe the effective total
peaks due to transport through the AI are added. The
transmission across the AI, T (E), as125
transmission through the AI can be rewritten as t,AI
T (E)=T (E)+T (E)=
(cid:104) (cid:105) t,AI AI R,AI
TAI(E)=Tr γ¯L,AI(E)G¯MAIB†(E)γ¯R,AI(E)G¯MAIB(E) , = γ¯L,AI(E)γ¯R,AI(E) Im(cid:2)−G¯MB(E)(cid:3).(79)
(70) γ¯ (E)+γ¯ (E) AI
L,AI R,AI
9
Eqs. (77-79) extend the results of Ref.125 to the more teractiontermHˆ [Eq.(43)],inCT-HYBoneimposes136
C
general case including background transmission and in-
terferenceterms,astypicallyfoundforSTMexperiments Hˆ1 =Hˆ¯AI,D+Hˆ¯TB+HˆI (81)
of molecules on surfaces. Hˆ =Hˆ¯ (82)
2 AI,NI
In this article we only consider applications to linear-
response transport, so that the conductance G = dI/dV with Hˆ¯ , Hˆ¯ , Hˆ and Hˆ¯ defined in subsection
AI,D TB I AI,NI
is given by
IID. Therefore the perturbation Hamiltonian in CT-
HYB is represented by the bath-AI coupling, while the
e2 (cid:90) (cid:18) df (cid:19) reference Hamiltonian is that of the decoupled atomic-
G = dE − T (E). (80)
h dE t,AI like correlated AI and of the isolated bath.
Having divided the Hamiltonian in two parts, one is
able to introduce the interaction picture, where an op-
This equation generalizes the Landauer formula used in
erator Oˆ depends on the imaginary time τ as Oˆ(τ) =
DFT+NEGF to the case of an interacting EM, with the
only difference that T (E) (used in the Landauer for- eτHˆ1Oˆe−τHˆ1, with 0<τ <β, and the partition function
0
mula) is replaced by T (E). Therefore, following the is written as the standard time-ordered exponential
t,AI
satnaanlydzairndgptrhaectizceerou-sbeidasintrDanFsTpo+rNt EpGroFpetrhtiaetscboynsmistesanins Z =Tr(cid:2)e−βHˆ1Tτe−(cid:82)0βdτHˆ2(τ)(cid:3)
opfretsheenttrthanesrmesisuslitosnfoartzeeqrou-ibliibarsiutmranTs0p(oErt),inheprreeswenecewiolfl =(cid:80)∞n=0(cid:82)0βdτ1...(cid:82)τβn−1dτnwn, (83)
many-body effects by plotting Tt,AI(E). with Tτ the time-ordering operator, and
w =Tr(cid:2)e−(β−τn)Hˆ1(−Hˆ )...e−(τ2−τ1)Hˆ1(−Hˆ )e−τ1Hˆ1(cid:3).
n 2 2
(84)
F. CTQMC impurity solver The partition function has the form of an integral over a
configurationspace. Insuchspace,anyparticularconfig-
uration is specified by the expansion order n, the times
In the present work the AIM is solved by us-
{τ ,...,τ } and a set of discrete variables, for instance
1 n
ing CTQMC for quantum systems in thermodynamic
the spin, and it is characterized by the probability dis-
equilibrium113, since we only address linear-response tribution p = w (cid:81)n dτ . It is this integral, which is
transport. In this case the method is well-established. n n k=1 k
ultimately evaluated by Monte Carlo techniques.
We note however that recently there have been a num- By using the definitions of Hˆ and Hˆ in Eqs (81) and
1 2
ber of developments towards the extension of CTQMC
(82), w in Eq. (84) becomes113,136
n
to out-of-equilibrium problems, both in steady-state and
time-dependent frameworks129–134.
w =
In this work, two different algorithms have been con- n
(cid:20) (cid:21)
stiimdeereadu:xtilhiaerywefiaekld-co(CupTl-inAgUaXp)p1r3o5,acahn,dctahlleedhycbornitdiinzuaotiuosn- ZBTr e−βHˆlocTτ(cid:89)dˆσ(τnσσ)dˆ†σ(τ(cid:48)σnσ)...dˆσ(τ1σ)dˆ†σ(τ(cid:48)σ1)
σ
expansion,strongcouplingapproach(CT-HYB)136. CT- ×(cid:89)detD−1(τσ,...,τσ ;τ(cid:48)σ,...,τ(cid:48)σ ). (85)
AUX scales as the product of the interaction U and of σ 1 nσ 1 nσ
theinversetemperature,sothatitsapplicationturnsout σ
too computationally demanding for Kondo systems at Here Z is the bath partition function, Hˆ = Hˆ¯ +
temperatures of the order of only a few Kelvin. In con- B loc AI,D
trast, we found CT-HYB able to provide quite accurate HˆI, dˆ(σ†)(τ) = eτHˆlocdˆ(σ†)e−τHˆloc and the matrix Dσ−1 has
results at a reasonable computational cost for the spe- elements
cific Au/TOV system considered in the following. Since (D−1) =∆˜f (τσ−τ(cid:48)σ), (86)
this system is described as a SIAM, here we present the σ ij AI,σ i j
method only for this case. However, we remark that our
which are the Fourier transforms of the hybridization
implementation can treat multi-orbital systems as well,
function
although only for density-density interaction terms.
∞
CTQMCisrestrictedtofinite-temperatures,wherethe ∆˜f (τ)= 1 (cid:88) e−iωnτ∆¯ (iω ), (87)
partition function is given by Z = Tr(cid:2)e−βHˆIAIM(cid:3). The AI,σ β AI,σ n
n=−∞
starting step, which is common to all algorithms, is to
separate the interacting AIM Hamiltonian [Eq. (42)] in evaluatedfortheimaginary-timeintervalτσ−τ(cid:48)σ,which
i j
two parts, a reference Hamiltonian Hˆ and a perturba- separate pairs of operators dˆ (τσ) and dˆ†(τσ) in Eq.
1 σ i σ j
tion Hamiltonian Hˆ , so that Hˆ = Hˆ +Hˆ . Each (85). In the function w the trace accounts for the im-
2 IAIM 1 2 n
CTQMC algorithm differs in the exact definition of Hˆ purity that fluctuates between different quantum states
1
and Hˆ 113. While in CT-AUX and other weak-coupling as electrons jump in and out, while the determinants re-
2
approaches Hˆ is set equal to the effective Coulomb in- sum the bath evolutions which are compatible with the
2
10
sequence of quantum fluctuations in the impurity. Note mation (89) can be written as a unitary transformation.
that here we have explicitly re-introduced the spin-index Second, only few Legendre coefficients are needed to
thatwasdroppedfromtheequationsintheprevioussub- expresstheMatsubaraGF.ThisisbecausetheLegendre
sections. However, if the bath is non-magnetic like in coefficients for a Matsubara GF decay faster than any
most cases, the hybridization function is the same for power of the Legendre expansion index l137. A valuable
spin-up and down, and can therefore be obtained from effect of this is that the statistical noise is filtered out
non spin-polarized calculations. due to the cutoff at a certain expansion order. After ex-
In the specific case of a one-orbital Hubbard interac- tensive tests for the specific system studied in this work,
tion, CT-HYB can be efficiently implemented by using we found that an appropriate choice for the cutoff of
the so-called “segment representation”136. This means the Legendre polynomials is around 100. However, such
that each configuration is depicted by segments, which cutoffmustgenerallybedeterminedforeachspecificcase.
represent time intervals τσ −τ(cid:48)σ during which an elec-
tron of a given spin resides on the impurity. Notably, Themany-bodyself-energyisobtainedeitherfromthe
with such representation the trace in w can be evalu- Dyson Eq. (45), or by using the expression138
n
ated in polynomial time. New configurations are then
obtained by either adding or removing segments. This Σ¯MB (iω )=UF¯AMIB,σ(iωn), (90)
isenoughtoensureergodicity,althoughotheroperations AI,σ n G¯MB (iω )
AI,σ n
(such as shifting the segments’ end-points) and global
updatesmustbeimplementedtoensureanefficientsam- where F¯MB(iω ) is the Matsubara representa-
AI,σ n
pling. Each update is accepted or rejected according to tion of the correlation function F˜MB(τ − τ(cid:48)) =
the Metropolis algorithm. The acceptance probability is AI,σ
−(cid:104)T dˆ (τ)dˆ† (τ(cid:48))nˆ (τ(cid:48))(cid:105), which can be easily com-
efficiently computed with standard fast matrix update τ −σ −σ σ
methods113. puted in the Legendre polynomial basis at no extra
cost139. This last approach usually provides much more
As seen in Eq. (86), CT-HYB requires the Fourier
accurate results than the calculation through the Dyson
transform of the hybridization function, which in prin-
equation. The inversion of the GF in Eq. (45) amplifies
ciple is calculated through a summation over an infi-
the statistical noise, in particular at high Matsubara
nite number of frequencies. However, in practice only
frequencies, when the difference between interacting
a finite number of frequencies smaller than a certain
and non-interacting GF is very small. In contrast, such
cutoff N can be inevitably summed up, although the
ω
high-frequencylimitof∆¯ (iω )determines∆˜f (τ)close problem is not present when using Eq. (90).
AI n AI Finally, we note that, although in this subsection we
to τ = 0. Therefore, in order to accurately calculate
have made explicit the dependence of the GF on the
∆˜f (τ) for any arbitrary τ, we use a standard approach
AI spin-index for notation completeness, the spin up and
that consists in adding and removing from Eq. (87) the
down GFs are equal as long as the substrate is non-
high frequency limit of ∆¯ (iω ). This limit is given by
AI n magnetic, and there is no Zeeman-like term in the AIM
−iM /ω , with M given in Eq. (A36) of the Appendix
1 n 1 Hamiltonian. Since this is the case for the application
A4. Hence we evaluate ∆˜f (τ) as
AI presented in this work, we will once again drop this
index in the following.
∆˜f (τ)= 2M1 (cid:88)Nω (cid:104)(cid:60)(cid:110)∆ (iω )(cid:111)cosω τ+
AI,σ β AI n n
n=0
G. Analytic continuation
(cid:16) (cid:110) (cid:111) M (cid:17) (cid:105) M
+ (cid:61) ∆ (iω ) − 1 sinω τ + 1, (88)
AI n ω n 2
n In order to compute the transport properties, as out-
where the summation is restricted to positive fre- linedinsubsectionIIE,weneedtheretardedmany-body
quencies, since (cid:60)(cid:110)∆ (iω )(cid:111) = (cid:60)(cid:110)∆ (−iω )(cid:111) and GFontheAI,G¯MB(E),whoseimaginarypartdefinesthe
AI n AI n AI
(cid:110) (cid:111) (cid:110) (cid:111) spectral function
(cid:61) ∆ (iω ) =−(cid:61) ∆ (−iω ) .
AI n AI n
1
During Monte Carlo sampling the many-body Matsub- AMB(E)=− (cid:61)G¯MB(E), (91)
AI π AI
ara GF can be directly estimated. However, following
Boehnkeet al.137,weuseanexpansionoftheG¯MAIB,σ(iωn) normalized as
in terms of Legendre polynomials, P [(2τ/β)−1],
l (cid:90) ∞
dE AMB(E)=1, (92)
G¯MB (iω )= AI
AI,σ n −∞
√
(cid:88)G¯MB (l) 2l+1(cid:90) βdτeiωnτP (2τ/β−1), (89) and which is related to the real part via the Kramers-
AI,σ β l Kronig relation
l≥0 0
andweestimatetheexpansioncoefficientsG¯MB (l). The (cid:90) ∞ AMB(E)
AI,σ (cid:60)G¯MB(E)=−P dE(cid:48) AI . (93)
advantage of this choice is twofold. First, the transfor- AI E(cid:48)−E
−∞