{'fraction_non_alphanumeric': 0.04197825919358569, 'fraction_numerical': 0.019273764551692236, 'mean_word_length': 3.4809569047413422, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 438, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Transfer-matrix summation of path integrals for transport through nanostructures

On the basis of the method of iterative summation of path integrals (ISPI), we develop a numerically exact transfer-matrix method to describe the nonequilibrium properties of interacting quantum-dot systems. For this, we map the ISPI scheme to a transfer-matrix approach, which is more accessible to physical interpretation, allows for a more transparent formulation of the theory, and substantially improves the efficiency. In particular, the stationary limit is directly implemented, without the need of extrapolation. The resulting new method, referred to as"transfer-matrix summation of path integrals"(TraSPI), is then applied to resonant electronic transport through a single-level quantum dot.


I. INTRODUCTION
Quantum-dot systems have been well studied both experimentally and theoretically for over thirty years. Their optical properties, namely the quantum size effect, make them useful for commercial applications within liquid crystal displays [1]. The tunability of their electrical properties allows to control single electrons [2] and gives rise to a number of effects such as Coulomb blockade, the Kondo effect, tunnel magnetoresistance or Andreev bound states [3], on which versatile electronic and spintronic quantum-dot devices such as a single-electron transistor [4,5], a quantum-dot spin valve [6,7] or a Cooper-pair splitter [8][9][10] are based.
In order to study resonant transport in interacting quantum-dot setups, a numerically exact method called "Iterative summation of path integrals" (ISPI), was developed [11,12,14]. ISPI is based on the truncation of correlations that decay exponentially in time. The method is well suited to study quantum-dot systems at finite temperature, including both equilibrium and nonequilibrium, in the regime in which various energy scales, e.g., associated with Coulomb interaction, temperature, or transport voltage, are of the same order of magnitude and, therefore, lack a clear separation. The ISPI scheme was first introduced to discuss nonequilibrium transport through the Anderson model [11,12]. Applying the method to the Anderson-Holstein model, where the quantum dot is coupled to a phonon mode, demonstrated the impact on the Franck-Condon blockade, when entering the quantum-coherent regime [13]. Recently, the ISPI method was applied to quantum-dot spin valves, demonstrating the importance of resonant effects in the tunnel-magnetoresistance as well as unveiling interaction-induced current asymmetries caused by an interaction-induced exchange field [14,15].
The purpose of this work is to develop the ISPI scheme further. We show that the necessary truncation of correlations motivates a mapping of ISPI to a transfer-matrix approach, which by construction is formulated in the stationary limit, such that extrapolation of finite-time results is not needed anymore. We develop the theoretical corner stones of this new method, referred to as "Transfer-matrix summation of path integrals" (TraSPI). In order to keep the discussion transparent, we exemplify the method for a system with relatively few degrees of freedom, namely the Anderson model describing a singlelevel quantum dot coupled to two normal metal leads. We, however, emphasize that all new concepts discussed in this work are not limited to this simple model, but can easily be transferred to other, more intricate setups.
Single-electron transistors that utilize a quantum dot as an island have gathered a lot of attention throughout the years and are still under heavy investigation, both experimentally and theoretically. To mention just a few recent examples, a single-electron transistor consisting of a quantum dot and normal metal leads was realized experimentally to demonstrate that shot noise in a single-electron transistor can be reduced significantly via feedback control, which should allow the construction of efficient, nanoscale thermoelectric devices [16]. Furthermore, if the quantum dot is periodically driven via a gate voltage, it is possible to accurately control the dot's emission time statistics [17]. For a system, in which the quantum dot is coupled to a single lead only, the quantum dot can be driven out of equilibrium via a plunger gate voltage, which allows the measurement of the free energy of a confined electron in order to study thermodynamics on the microscopic level [18]. For a superconducting single-electron transistor, an attractive interaction was found that survives even far beyond the superconducting regime [19,20]. Employing the method of full-counting statistics for a negative-U Anderson model, it was shown that this phenomenon is robust, even for fast spin relaxation [21].
On a theoretical basis, different approaches are used and actively developed to study different parameter regimes of quantum-dot systems. The method of Dynamical Mean Field Theory is advanced and combined with other methods, such as functional renormalization group theory, to increase the predictive power of the method, even for strong and nonlocal electron correlations [22]. Perturbation theory in the tunnel coupling strength Γ within a master-equation approach, have proven highly useful in the description of quantum-dot systems. This method is developed further in different directions, e.g. by introducing SU(N )-invariant kinetic equations to ef-fectively study multilevel quantum dots [23] or by improving on the commonly used rotating-wave approximation, leading to a so-called coherent approximation [24]. While perturbative methods often allow for at least a qualitative description, nonperturbative effects are, by construction, beyond their scope. To cover them, numerically exact methods are in high demand. Several approaches are known to tackle this problem. Quantum Monte Carlo simulations were advanced to reach the stationary regime for systems in nonequilibrium [25,26]. Different flavors of renormalization group theory (RG) have been applied to quantum-dot systems almost since the inception of their theoretical discussion [27]. Since then significant advances have been made to the formalism: A combination of numerical RG and time-dependent density-matrix RG allows to discuss the nonequilibrium steady state transport properties of quantum-dot systems [28], while within functional RG it was possible to approximate the flow of the Luttinger-Ward functional while maintaining conservation laws [29]. Finally, it was shown that density functional theory is able to study outof-equilibrium transport theories, even in strongly correlated systems, such as the Anderson model [30], while a quasiparticle Fermi-liquid theory can be used to work within the low-energy limit of such systems [31].
The article is structured as follows. In Sec. II we introduce the Anderson model's Hamiltonian and derive the path-integral formulation of the generating function. We demonstrate how interactions are decoupled via a discrete Hubbard-Stratonovich (HS) transformation and then solve the path-integral. In Sec. III we discuss the main ideas behind the ISPI and TraSPI schemes, namely the truncation of exponentially decaying correlations after a memory time t K , as well as the truncation of the correlations induced by the HS transformation. We give a short overview of the ISPI formulation that was introduced in earlier works in Sec. III B. We then demonstrate how it can be mapped to a transfer-matrix approach, and discuss the main benefits of this new formulation in Sec. III C. We finish this section with a description of the convergence procedure that ensures that the results obtained via the TraSPI scheme are numerically exact. In Sec. IV we discuss the results for the Anderson model, obtained via TraSPI. First, we discuss current based observables, like the current and the conductance, and after that the dot's occupation number. Finally, we conclude in Sec. V.

II. MODEL
We write the well-known Hamiltonian for an interacting, single-level quantum dot that is tunnel-coupled to two metallic leads in the form [11,32] (we set = 1 throughout this work) The on-site occupation-number operator is given byn σ = d † σdσ , whered † σ andd σ create or annihilate an electron on the quantum dot with spin σ = ↑, ↓, respectively. The Coulomb interaction strength is given by U . In Eq. (1) we made use of the operator identitŷ n ↑n↓ = 1 2 (n ↑ +n ↓ ) − 1 2 (n ↑ −n ↓ ) 2 and incorporated the terms linear inn σ by shifting the energy of the quantum dot's level, such that 0,σ = 0 +σB/2 with 0 = E 0 +U/2, where the bare energy level E 0 in the absence of magnetic field and Coulomb interaction can be tuned via a gate voltage. To write the interaction in terms ofn ↑ −n ↓ turns later out to be advantageous for the discrete Hubbard-Stratonovich transformation. An electron with energy αk = k − µ α in lead α = (L, R) and with momentum k is created or annihilated by the operatorsĉ † αkσ andĉ αkσ , respectively. Finally, t α denotes the tunnel coupling between lead α and the quantum dot. The tunnel-coupling strength between quantum dot and lead α is given by denotes the density of states of lead α at the Fermi level. We work in the wide-band limit, which usually is a good approximation in the stationary regime [33]. We also assume a symmetric setup, where Γ = Γ L = Γ R and where the chemical potential of the leads are given by the bias voltage µ α = ±eV /2, for the left and right lead, respectively.

A. Path-integral formulation
We are interested in discussing time-local observables, namely the current, the quantum dot's occupation number and its spin-projector expectation value in the stationary limit. The latter is achieved for t b → ∞, < l a t e x i t s h a 1 _ b a s e 6 4 = " d g 5 Y u I t J j M w R A g 2 f s D o b a r Q b U y w = " > A A A C A 3 i c b V D L S g N B E O z 1 G e M r 6 t H L Y h A i a N g V U Y 9 B L 5 4 k g n l A s o T Z y S Q Z M j u 7 z P Q K y 5 K j H + B V P 8 G b e P V D / A J / w 0 m y B 5 N Y 0 F B U d d P d 5 U e C a 3 S c b 2 t p e W V 1 b T 2 3 k d / c 2 t 7 Z L e z t 1 3 U Y K 8 p q N B S h a v p E M 8 E l q y F H w Z q R Y i T w B W v 4 w 9 u < l a t e x i t s h a 1 _ b a s e 6 4 = " Q q 9 E 7 z 7 6 r 1 0 e u 7 Y P j F K z + 6 H y p R E e 6 r + n U h J o H U S + K Y z I D j U i 9 5 E / M 9 r x 9 i / 8 V I u o x i Z p L N F / V j Y G N q T 3 + 0 e V 4 y i S A w h V H F z q 0 2 H R B G K J q G 5 L d E w 0 Z z q c c E k 4 y 7 m s E w a F x X 3 q u I + X J a q t 1 l G e T i C Y y i D C 9 d Q h X u o Q R 0 o j O A F X u H N e r b e r Q / r c 9 a a s 7 K Z Q 5 i D 9 f U L m q C X a g = = < / l a t e x i t >

· · ·
< l a t e x i t s h a 1 _ b a s e 6 4 = " f 9 N 6 P C p G s 3 m P J V N Y V a q O f S k K D a E = " > A A A B / X i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K q M e A l x w T M I m Q L G F 2 0 k m G z M 4 u M 7 P C s g Q / w K t + g j f x 6 r f 4 B f 6 G k 2 Q P J r G g o a j q p r s r i A X X x n W / n c L G 5 t b 2 T n G 3 t L d / c H h U P j 5 p 6 y h R D F s s E p F 6 D K h G w S W 2 D D c C H 2 O F N A w E d o L J / c z v P K H S P J I P J o 3 R D + l I 8 i F n 1 F i p a f r l i l t 1 5 y D r x M t J B X I 0 + u W f 3 i B i S Y j S M E G 1 7 n p u b P y M K s O Z w G m p l 2 i M K Z v Q E X Y t l T R E 7 W f z Q 6 f k w i o D M o y U L W n I X P 0 7 k d F Q 6 z Q M b G d I z V i v e j P x P 6 + b m O G d n 3 E Z J w Y l W y w a J o K Y i M y + J g O u k B m R W k K Z 4 v Z W w s Z U U W Z s N k t b 4 n G q O d P T k k 3 G W 8 1 h n b S v q t 5 N 1 W t e V 2 r 1 P K M i n M E 5 X I I H t 1 C D O j S g B Q w Q X u A V 3 p x n 5 9 3 5 c D 4 X r Q U n n z m F J T h f v 8 Y P l f w = < / l a t e x i t > t < l a t e x i t s h a 1 _ b a s e 6 4 = " U B I r 6 t v a r f L i P H d w 1 y N X S J 5 G X a W g 0 l 8 U P B 4 r 4 q q e l 4 k u A b b / r Q y K 6 t r 6 9 n c R n 5 z a 3 t n t 7 C 3 3 9 B h r C i r 0 1 C E q u U R z Q S X r A 4 c B G t F i p H A E 6 z p D W / H f v O B K c 1 D e Q 9 J x N y A 9 C X 3 O S V g p G Y J u s H Z 6 U m 3 U L T L 9 g R 4 m T g z U q z g 2 s 9 X L n t Z 7 R a + O 7 2 Q x g G T Q A X R u u 3 Y E b g p U c C p Y K N 8 J 9 Y s I n R I + q x t q C Q B 0 2 4 6 O X e E j 4 3 S w 3  / r Q y K 6 t r 6 9 n c R n 5 z a 3 t n t 7 C 3 3 9 B h r C i r 0 1 C E q u U R z Q S X r A 4 c B G t F i p H A E 6 z p D W / H f v O B K c 1 D e Q 9 J x N y A 9 C X 3 O S V g p G Y J u s H p 2 U m 3 U L T L 9 g R 4 m T g z U q z g 2 s 9 X L n t Z 7 R a + O 7 2 Q x g G T Q A X R u u 3 Y E b g p U c C p Y K N 8 J 9 Y s I n R I + q < l a t e x i t s h a 1 _ b a s e 6 4 = " f H H u X W e m n 4 8 9 W 9 d e G S 0 C g Y f F a 0 w = " > A A A C A n i c b V D L S s N A F L 3 x W e u r 6 t L N Y B F c l U R E X R Z E c F n B P q A N Z T K Z t G M n M 2 F m I o T Q n R / g V j / B n b j 1 R / w C f 8 N p m 4 V t P X D h c M 6 9 3 H t P k H C m j e t + O y u r a + s b m 6 W t 8 v b O 7 t 5 + 5 e C w p W W q C G 0 S y a X q B F h T z g R t G m Y 4 7 S S K 4 j j g t B 2 M b i Z + + 4 k q z a R 4 M F l C / R g P B I s Y w c Z K r R 4 J p d  where t b = t m − t 1 (before) is the time interval between the initialization and the measurement time t m , while t a = t N −t m (after ) denotes the time interval between t m and the Keldysh return time t N . Nonequilibrium properties are taken into account within a functional integral formulation on the Keldysh contour C (see Fig. 1) [14,34,35]. To take the forward and backward branch of the Keldysh contour into account, it is useful to define τ = (t, ν), with physical time t and Keldysh branch index ν = ±, with + and − representing the upper and lower Keldysh contour, respectively.
For any time-local observableÔ at measurement time t m , we introduce a source term ηO(t m ), with η ∈ R, in order to break time-translation symmetry on the level of the system's action S in a specific way. This allows us to calculate the expectation value of said observable via a derivative of the system's Keldysh generating functional, The Keldysh generating functional for the Hamiltonian given in Eq. (1) takes the form and fulfills Z[0] = 1 by construction. The functional integral is built in the basis of fermionic coherent states |Ψ(τ ) , being defined via the eigenvalue equations of the fermionic annihilation operatorŝ As a result, the action S and the source term ηO(t m ) in Eq. (3) are functions of the Grassmann fields d σ (τ ), d σ (τ ),c αkσ (τ ), c αkσ (τ ), and the functional integral runs over all of these degrees of freedom [12,14,34]. Dropping the explicit τ -dependency, the action takes the form Performing the integral in Eq. (3) does not pose any challenge on the terms contributed by the leads' and the tunneling Hamiltonian, since they are quadratic in the Grassmann fields. However, the quartic on-site interaction term has to be tackled via a discrete Hubbard-Stratonovich transformation [36][37][38]. For this, we first discretize the Keldysh contour into 2N time slices of length δ t (see Fig. 1), and then perform a HS transformation, on each of these 2N slices, keeping in mind that n σ can only assume the values 0 and 1. This decouples the interaction term, at the cost of introducing one Ising-like degree of freedom, s = ±1, per time slice. The HS parameter ζ ν is determined uniquely for 0 ≤ δ t U < π via [12,14] cosh ζ ν = e − 1 2 iνδtU .
As a result, we have introduced 2N new HS spins but are able to also integrate over the dot degrees of freedom, solving the functional integral in Eq. (3). We finď For this we used the HS spin vector such that the sum includes all 2 2N possible spin configurations along the discretized Keldysh contour. In addition, we introduced several new matrices: We identify the inverse time-discrete Green's function ∆ −1 = ∆ −1 0 − Σ T of the non-interacting setup, where ∆ 0 is the free dot's Green function, and Σ T = α Σ T,α is the tunneling self energy. In addition, we introduced the charging self energy SΣ C with the diagonal spin matrix S = diag( s )⊗σ z , which is the only part depending on the HS spins s, as well as the source self energy Σ O , which is included to account for the source term. We also made use of the short-handsΣ C = Σ C ∆ andΣ O = Σ O ∆ for later convenience. All of these matrices have dimensions 4N × 4N due to the Trotter slicing (see Fig. 1) and the spin degree of freedom.
Note that we modified the generating functional by absorbing the factor 1 2 per spin from Eq. (6). Additionally, we have multiplied S from the left and ∆ from the right in Eq. (9). These changes do not affect the expectation value of the observable, since det S = 1 and det ∆ = const, which cancels due to the logarithmic derivative in Eq. (2). However, multiplying by ∆ ensures that D[η, s ] decays exponentially, while through the multiplication with S the HS spins are located only on the diagonal of D[η, s ] as well as on the parts affected by the source self energy. The first property will be crucial when implementing the ISPI scheme, while the second is useful for an efficient implementation of the transfer-matrix formulation.

B. Form of the matrices
To specify the elements of the matrices in Eq. (9a) we employ the basis (n, ν, σ), where the Trotter index n = 1, . . . , N labels the time slice, the Keldysh index ν = ± distinguishes the upper from the lower Keldysh contour, and σ = ↑, ↓ denotes the spin. In this basis, the Green's function of the non-interacting quantum dot in the presence of leads is given as [54] and σ z and σ 0 are the Pauli matrices, acting on either Keldysh or spin space if they appear on the left or on the right of the tensor product, respectively. The Fermi function f α (ω) = [exp(β(ω − µ α )) + 1] −1 describes the equilibrium occupation distribution of lead α. Note that due to the symmetric discretization of the derivative ∂ t → ω −1 in frequency space, the discretized advanced and retarded Green's functions have the diagonal |∆ a,r | nn = 1 2 instead of 1, such that det(2∆) = 1. The charging self energy SΣ C is time-local and therefore a diagonal matrix, with The form of the source self energy Σ O is based on the observable of interestÔ, since it is derived from its source term ηO(t m ). The source self energies for the occupation number, the spin projection in z direction, and for the current were derived in earlier works [11,12,14,15]. As a result, we only present the results here, and refer to the aforementioned references for a more detailed derivation. Assuming we include the measurement on the Trotter slice m = t m /δ t , we find for the occupation numberN = σn σ and for the spin projectionŜ z = 1 2 (n ↑ −n ↓ ) [15] These are sparse matrices, containing only two nonzero elements. Finally, the current operator is given byÎ = −ie/2 αkσ αt α (ĉ † αkσd σ −d † σĉαkσ ), with its source term given by [11,12,14] The dotted line is an exponential fit for the envelope function according to |∆| ∝ exp(−|t − t |/ξ0). Inset: Correlation time ξ0 of the envelope function for the non-interacting system as a function of temperature kBT . Other parameters are the same as in the main panel.
Note that the source self energy for the current operator is a sparse matrix, too, with only one row m having nonvanishing elements.

III. METHOD
A. Truncation of matrix D The sum in Eq. (8) runs over all configurations of the 2N HS spins, with N usually being of the order of several hundreds. Summing over these 2 2N configurations is an insurmountable task, and approximations are in order. The approximation for the ISPI scheme is based on the fact that lead-induced correlations decay exponentially with time at finite temperatures [11,12]. As a consequence the system's non-interacting Green's function ∆and with it D[η, s ] -also decays exponentially with time |t − t |. This is shown in the main panel of Fig. 2, where the absolute value of ∆ +− is plotted against |t−t | for the parameter set k B T = 0.2 Γ, 0 = 0, eV = Γ, B = 0 and δ t → 0. While for low temperatures one finds large oscillations as seen in the figure, the enveloping function still decays exponentially, as long as k B T > 0. We demonstrate this temperature dependence of the correlations' decay in the inset of Fig. 2, where the correlation time ξ 0 of the enveloping function is plotted as a function of temperature.
This motivates us to truncate the interacting Green's function Eq. (9), such that where = 1, . . . , L and L = N/K. The number of Trotter slices within one block, K = t K /δ t , depends on a chosen memory time t K and the Trotter step size δ t . Since only row m in Eq. (9) is affected by the source self energy, see Eqs. (14), this also holds true for the block matrix, where only the block rowm containing row m is affected by the source self energy. We also introduced the block spins s , each consisting of a set of 2K HS spins, such that s = ( s 1 , . . . , s L ), see Eq. (10).
Using this truncation we are now able to formulate the ISPI scheme, which was first applied to the Anderson model in Ref. 11 and later employed to study the Anderson-Holstein model [13], and the quantum-dot spin valve [14,15]. We introduce the main ideas in the next section, and refer to any of the aforementioned references for a more complete discussion. After that we introduce the mapping of the ISPI method to a transfer-matrix method, which by construction directly addresses the stationary limit and drastically increases the computational performance.

B. ISPI formulation
As can be seen in Eq. (8) one has to calculate the determinant of the block-tridiagonal matrix from Eq. (15). According to Ref. 39 such a determinant can be evaluated iteratively viaŽ where s k: = { s k , . . . , s }, anď denotes the Schur complement in the -th step, witȟ C 1 ( s 1:1 ) = D 1 ( s 1 ) and = 2, . . . , L. Therefore, eacȟ C depends on all previousČ k , with k < , and since the spins are distributed diagonally, Eq. (17) connects s with all previous s k . To remain consistent with the idea of truncating correlations after the memory time t K , we approximateČ ( s 1: ) by [11,14] C that is, in Eq. (17) we replaceČ −1 −1 ( s 1: −1 ) with D −1 −1 ( s −1 ), which only depends on s −1 . Therefore, C only connects s −1 and s , effectively truncating interaction-induced correlations after the memory time t K . As a result, we are able to rewrite Eq. (16), findinǧ Note that we used the fact that each C depends only on the two block spins s −1 and s , allowing us to evaluate L sums over M = 2 2K spin configurations, instead of the much larger sum over 2 2N configurations. Since we can choose K N , this is a huge reduction of complexity, allowing us to evaluate the generating functionalŽ[η] for much larger values of N .
ISPI can be characterized as a finite-time implementation due to the finite length of the considered Keldysh contour as depicted in Fig. 1. In order to reach the stationary limit, the time t b before the measurement has to be chosen large enough such that the system has relaxed from any arbitrary initial state to its stationary state. This can be ensured by choosing Γt b 10, depending on the system under consideration. For the Anderson model, we find that a time interval of Γt b = 15 is sufficient.
The finite-time implementation has the disadvantage that, first, one needs to carefully choose for each calculation the proper time interval to make sure that the stationary limit has been achieved. Second, increasing the time interval (due to larger relaxation times) increases the computational cost. The TraSPI formulation, put forward in this paper, is motivated by the desire to directly perform the limit t b → ∞ analytically by using transfer matrices. This strategy does not only strongly decrease the computational cost but has also other benefits, as described in the next section.

C. TraSPI formulation
The term "transfer matrix" is reminiscent of the transfer-matrix method known from statistical physics, which was first introduced to solve the 1-dimensional Ising model [40,41] and later used by Onsager as the basis for his well-known exact solution of the two-dimensional Ising model [42,43], for a recent application see, e.g., [44,45]. In fact, Eq. (19) is the basis for a mapping of the ISPI scheme to a transfer-matrix formulation, which we derive here.
First, we reiterate that each C depends on two block spins s −1 and s , while D 1 depends only on s 1 . Therefore, the shape of Eq. (19) suggests rewriting it into a matrix product in the space of HS-spin configurations. For this, we enumerate the M different configurations of the block spins s by µ = 0, . . . , M −1, such that, e.g., µ = 0 corresponds to 2K HS spins pointing up, s = (1, . . . , 1).

d K d A y 8 X J S h h z 1 X u m n 2 5 c k i a g w h G O t O 5 4 b G z / D y j D C 6 b j Y T T S N M R n h A e 1 Y K n B E t Z 9 N T x 2 j U 6 v 0 U S i V L W H Q V P 0 7 k e F I 6 z Q K b G e E z V A v e h P x P 6 + T m P D a z 5 i I E 0 M F m S 0 K E 4 6 M R J O / U Z 8 p S g x P L c F E M X s r I k O s M D E 2 n b k t 8 T D V j O h x 0 S b j L e a w T J r n V e + y 6 t 1 f l G s 3 e U Y F O I Y T q I A H V 1 C D O 6 h D A w g M 4 A V e 4 c 1 5 d t 6 d D + d z 1 r r i 5 D N H M A f n 6 x d z n Z Z J < / l a t e x i t > (b)
< l a t e x i t s h a 1 _ b a s e 6 4 = " J l 9 g Z 8 7 p u U Q t b X z a N r / r t x f d Y 6 (27), respectively. The TM T, shown as orange arrows, first propagates one step of size tK backwards in time, followed by two steps in forward direction. The TM U on the other hand is represented by blue arrows, and propagates two steps of size tK forward in time, while V or V −1 are shown by green lines, propagating one step forward or backward, respectively.

T < l a t e x i t s h a 1 _ b a s e 6 4 = " V P 7 i p f S y J G K o r 0 / / A n M i 2 3 + 6 O E o = " > A A A C / X i c h V L L T t t A F D 2 Y P g K l J Y U l G 4 u o E i v L o Q 6 Q H e I l N k g g N Y B 4 q B q b S b D i l + w J U o g i f o I t 7 N i h b v s t f A U / w I I z U w e p i 4 i x x v f O u f e e O X d m / C w K C + W 6 T x P W 5 I e P n z 5 X p q a / z H z 9 N l v 9 P n d Y p L 0 8 k K 0 g j d L 8 2 B e F j M J E t l S o I n m c 5 V L E f i S P / O 6 m j h 9 d y b w I 0 + S X 6 m f y P B a d J G y H g V C E T s 5 i o S 7 9 9 q A 1 / F 2 t u Y 7 b a H q r y 7 b r N N x 6 0 / t J p 9 l c 8 1 Y a d t 1 x z a i h H P t p 9 R l n u E C K A D 3 E k E i g 6 E c Q K P i d o g 4 X G b F z D I j l 9 E I T l x h i m r U 9 Z k l m C K J d / j t c n Z Z o w r X m L E x 1 w F 0 i z p y V N n 5 w 7 h h G n 9 l 6 V 0 m / o H 3 h v D Z Y Z + w O A 8 O s F f Z p f T J O G c Y 9 4 g q X z H i v M i 4 z R 1 r e r 9 R d K b S x Z r o J q S 8 z i O 4 z e O P Z Y i Q n 1 j U R G 9 s m s 0 M O 3 6 y v e A I J b Y s K 9 C m P G G z T 8 Q W t M F Y a l q R k F O T L a f X p a z 3 j u 8 v Y Q 5 8 R f U s F 9 I M Y 3 b o 9 3 j l c d u o r j n f g 1 d Y 3 y q d R w Q I W s c T 7 X 8 U 6 d r F P x Q H 1 3 O I O 9 9 a N 9 W A 9 W n / + p V o T Z c 0 8 / h v W 3 1 e P t p z K < / l a t e x i t >
In addition, we write f (µ) instead of f ( s ) for any function that depends on the HS spins, in order to keep the notation compact. With this, we define the M ×M transfer matrix (TM) Thus, each row corresponds to one of the M configurations µ of the block spin s −1 and each column to one of the M configurations µ of s . If we additionally define the two vectors then the Keldysh generating functional (19) takes the simple matrix-product form and each multiplication with T −1, corresponds to one term s det C ( s −1 , s ) in Eq. (19).
In the next step, we symmetrize Eq. (22) by introducing two new transfer matrices U −1, and V . For this, we reinspect the definition (18) of C , whose determinant provides the elements of T −1, , noting that it has the well-known form of the Schur complement of a part of the tridiagonal block matrix D[η, s ], given by This is a 2 × 2 block matrix, affected by the block spins s −1 and s . Thus, we are able to rewrite det C as the quotient of two determinants, Returning to the transfer-matrix formalism, this corresponds to the matrix product where U −1, is a dense M × M transfer matrix and V is a M × M diagonal matrix, Building the generating functional with these, we get a symmetrized version of Eq. (22), where we used the fact that v| = 1| V 1 . We emphasize that to this point the Eqs. (19), (22), and (27) are synonymous, and no additional assumptions have been made. In Fig. 3 we show a sketch of the elements of D[η, s ] that are taken into account for the case L = K = 5. Each of the small 1 × 1 boxes represents a single time slice of D[η, s ] (thus a 4 × 4 block, to take spin and Keldysh into account.). The number of matrices entering Eq. (27) increases linearly with L. To achieve the stationary limit, we have to choose L large. We now improve upon this: (i) we perform the stationary limit analytically, (ii) we optimize the position of the measurement, (iii) we calculate the derivative of ∂ ∂ηŽ L [η] at η = 0 analytically, and (iv) we implement a differential measurement method.

Stationary limit
From this point on we are interested in the case that the system reached the stationary limit, i.e., we perform the limit t b → ∞. In order to perform this limit, we start with an index shift bym in the block indices ( ,m) old → ( ,m) new = ( −m, 0), such that the measurement block becomesm = 0 and = −L b − 1, . . . , L a + 1, with L b + L a = L − 3. Note that L b and L a denote the number of TMs T before and after the measurement, respectively.
In the next step, we make use of the fact that the system is symmetric under time-translation as long as the source term is not included, i.e., as long as η = 0. As a consequence, the non-interacting Green's function ∆ obeys ∆ n,n = ∆ n+1,n +1 . Thus, for the tridiagonal matrix D[0, s ] we equally find that D (µ) and D ± are independent of . With this and Eqs. (18), (23), we find that due to time-translation symmetry the transfer matrices fulfill T −1, = T, U −1, = U, and V = V for all = 0, 1. However, the source term was implemented to specifically break this time-translation symmetry, and therefore whenever an element of D[η, s ] is affected by the source self energy Σ O , above relations do not hold anymore. We mentioned before that for time-local observables, the source self energy affects only the matrices D ± 0 and D 0 ( s 0 ) at most. Consequently, only the transfer matrices acting at positionm = 0, i. e., T −1,0 , T 0,1 , U −1,0 , U 0,1 and V 0 , are affected by the source term and are therefore different from the others, see Fig. 3. We finď To reach the stationary limit t b → ∞, we now can take the limit L → ∞, meaning we let both L b → ∞ and L a → ∞. For the high powers of T we use the identity where λ k are the eigenvalues of T = V −1 U, with |λ 0 | > |λ 1 | ≥ . . ., while λ k | and |λ k are the respective left and right eigenvectors, fulfilling λ k |λ k = δ kk . Likewise, we define λ k | = λ k | V −1 and |λ k = V |λ k as corresponding left and right eigenvectors of UV −1 (with the same eigenvalues,λ k = λ k ). Plugging this back into Eq. (28) and performing the stationary limit, we find In Fig. 4(a) we plot the absolute value of the largest five eigenvalues of the TM T as a function of the Coulomb interaction strength U for the parameter set k B T = 0.5 Γ, 0 = 0, eV = 0.1 Γ, and B = 0 memory time Γt K = 1.5, and K = 5. We find that the largest eigenvalue is λ 0 = M at U = 0, as would be expected, since in the noninteracting limit the transfer matrix T is just a M × M matrix, where each element is equal to 1. For finite Coulomb interaction it scales with δ u = δ t U as for small δ u . However, once U becomes large, we find occasionally peaks where one of the lower eigenvalues diverges. These divergencies are not physical but are a consequence of the discretization and truncation of the Green's function. We can track them down to elements of the diagonal matrix TM V that vanish, resulting in divergencies in V −1 , and hence in T. Nonetheless, the physical correct eigenvalue λ 0 is still present, and we choose that one (instead of the peaks) to calculate the correct generating functional. To choose the correct eigenvalue λ 0 we make use of the fact that the corresponding right eigenvector is |λ 0 = |1 + O(δ 3 u ) (see below). Thus, we identify the correct eigenvalue by maximizing the overlap of the corresponding right eigenvector with |1 . In Fig. 4(b) we plot this physically correct eigenvalue λ 0 as a function of U for different K. As can be seen, only for U 3 Γ the data gets noisy, but increasing K allows for the calculation of larger U .
The TraSPI formulation is a significant reduction in complexity compared to the traditional, i.e., finite-time ISPI implementation: Instead of L−1 dense transfer matrices, we only have to evaluate the three dense TMs U, U −1,0 and U 0,1 , as well as the two diagonal TMs V and V 0 . An additional benefit of the UV-decomposition (26) is that it allows for an analytic evaluation of the derivative with respect to η. Before, however, we demonstrate that we are able to shift the measurement time to the end of the Keldysh contour in order to reduce the number of necessary transfer matrices even further.

Position of the measurement
In Eqs. (30), we assumed that the measurement is placed somewhere in the middle of the Keldysh contour (see also Fig. 1). However, based on causality, whatever happens at physical times after the measurement must not have an impact on the outcome of the measurement itself. In other words, the time propagation t a along the Keldysh contour from t m to t N and back to t m is unitary and should therefore cancel out, see Fig. 1. As a consequence, it should not be detrimental to shift the measurement time forward in time on the Keldysh contour, until it is located at the rightmost point. In the words of the TM formulation from Eq. (28), it should be sufficient to let L b → ∞ and set L a = 0, or when taking a look at the results from Eqs. (30), the vector |1 should be a right eigenvector of the TM T.
However, this exact unitary present in the continuum limit (3) is violated by the Trotter discretization. Nevertheless, the error in the right eigenvector is quite small, |λ 0 = |1 + O(δ 3 u ), and can be safely neglected, while the other eigenvectors show a stronger dependency on δ u . This means that we are able to rewrite Eqs. (30) as Note that with this equation, the measurement takes place in the second to last 4K × 4K block. If we actually measure on the last possible Trotter slice, t m = t N , only a single TM remains that is affected by the source term, resulting in the even simpler expression This means, we are able to further reduce the number of necessary TMs to the two fully occupied TMs U and U −1,0 , and one diagonal TM V. We now turn to deriving the analytic derivative with respect to η, necessary to calculate observables from Eqs. (33).

Analytic derivative
To calculate expectation values of observables via Eq. (33b), we make again use of Eq. (2). We explicitly calculate the derivative with respect to η here. If for some reason we do not wish to position the measurement at the end of the Keldysh contour, meaning we use Eq. (30b) instead of Eq. (33b), calculations for the derivatives of the two remaining TMs work analogous as presented here. Performing the derivative yields where we wrote A for ∂ ∂η A| η=0 . The derivative of a matrix acts on its elements, which in this case are themselves determinants of D −1, . Therefore, we make use of the identity (ln det A) = (det A) / det A = tr(A −1 A ), and extend the notation introduced for D −1, , Eq. (23), to other matrices, meaning A −1, denotes a 2 × 2 part of a L × L block matrix A. We further simplify the notation by writing A [2] = A −1, if = 0, 1, which are 2 × 2 block matrices not affected by the source term. With this we find for the derivative [55] U −1,0 = det D [2] (µ, µ ) tr where we plugged in D [2] (µ, µ ) and D −1,0 (µ, µ ) from (9), and pulled the determinant inside the trace. Inserting Eq. (35b) back into Eq. (34) allows for an analytic expression of the derivative of the generating functional, and with it for the expectation value. When placing the measurement somewhere in the center of the Keldysh contour, c.f. Eq. (30b), one would find that where Ô −1,0 is given by Eq. (34), while Ô 0,1 = λ −1 0 λ 0 |U 0,1 |λ 0 and Ô 0 = λ 0 |V 0 |λ 0 . The respective derivatives are then calculated in analogy to Eq. (35b).
Using the analytic derivative instead of a numeric derivative reduces numerical errors and of course allows for a more straightforward implementation. We now continue to introduce a differential measurement, to decrease the impact of numerical errors.

Differential measurement
In a final step, we minimize discretization errors further by only calculating a differential form of the observable of interest. This means, instead of the U -dependent observable we only calculate the difference between the interacting case and the noninteracting limit numerically, and get an improved estimate for the considered observable. The analytic expectation values O (0) for U = 0 are calculated by employing the Meir-Wingreen formula [49,50]   number and z projection of the spin, leading to where we defined p(ω) = [Γ 2 + (ω − 0 ) 2 ]/(ω − 0 ). Since it is expected that the numerical calculation of the observable at U = 0 has similar errors as in the interacting case, using the differential measurement cancels these errors. This leads to a significant reduction of numerical errors. The remaining errors are then effectively eliminated during a convergence procedure, which is discussed in the next section.

D. Convergence
Throughout the derivation of Eq. (34), we introduced two systematic errors; one being the Trotter error, caused by the finite discretization length δ t [46], the other is the error introduced by the truncation at memory time t K . However, both can be eliminated using a convergence procedure which allows us to provide numerically exact data. In earlier works [11,15], different approaches were used for this convergence procedure, the most common being a two-step regression, first for δ t → 0, and then with the resulting values for 1/t K → 0. We refer to the aforementioned sources for a detailed discussion of these procedures. The main problem of these procedures is that the second regression 1/t K → 0 using a power series ansatz is difficult to motivate. Therefore, we employ, in this work, a more sophisticated convergence procedure, that also starts with a power series regression of δ t → 0 but then uses Aitken extrapolations [47].
The complete process is shown in Fig. 5 for the example of the current as an observable and for the parameter set k B T = 0.5 Γ, 0 = Γ, eV = 0.5 Γ, and U = 1.5 Γ.
First, we calculate the expectation value O of the ob-servableÔ for a fixed memory time t K and varying step size δ t . The result is a set of different realizations of the same observable for this specific parameter set and memory time, O(t K , δ t = t K /K) with K = 3, . . . , 8. It is well known that Trotter errors are of the order δ 2 t [46]. Consequently, we fit O(t K , δ t ) against a polynomial expression such that the observable with eliminated Trotter error O(t K ) is given by the constant c 0 (see Fig. 5(a)). Note that in the equation above it is sufficient to stop at n = 2, thus only taking up to second order in δ 2 t into account [11,15].
Having eliminated the Trotter error, we now turn to eliminate the truncation error. For this, we repeat the first step for different values of the memory time t K , with 1 ≤ Γt K ≤ 2. This leads to a set of realizations of the desired observable O(t K ). As can be seen in Fig. 5(b), the observable as a function of Γt K converges exponentially against a limiting value for t K → ∞, which is the numerically exact result of the observable for one specific parameter set. This behavior can be understood from the exponential decay of the Green's function ∆, cf. Fig. 2. For such exponentially converging sequences the Aitken extrapolation works exceptionally well, accelerating the convergence, which eventually leads to an approximately constant sequence at the limiting value. The Aitken extrapolation for a sequence f n is given by [47] with forward differences [48] ∆ n f n = f n+1 − f n . For f n we use the set of realizations of the desired observable O(t K ) and find, that after two Aitken extrapolations the data is approximately constant around the numerically exact result for sufficiently large values of t K . We use the mean value of this sequence as the final result and its standard deviation as an error estimate.

IV. RESULTS
Having introduced and extensively discussed the TraSPI method, we now use it to calculate current-based and occupation number-based observables of a singlelevel, interacting quantum dot coupled to two normal leads. To make interaction-induced effects more visible, we show derivatives of the current and the occupation number with respect to the bias voltage eV . These derivatives are calculated numerically via central differences, where we choose eδ V ≤ 0.01 Γ.
For all data sets we took memory sizes 1.0 ≤ Γt K ≤ 2.5 and K = 3, . . . , 7 into account. Whenever this was not sufficient to reach convergence, we included K = 8.

A. Differential conductance
We start the discussion of the conductance in the linear-response regime with Fig. 6. There, the dI/dV is shown as a function of the gate voltage 0 for a small bias voltage eV = 0.1 Γ. The other parameters are given by k B T = 0.2 Γ and B = 0. Each curve represents a different value of the Coulomb interaction strength U , starting at U = 0, which is calculated analytically, see Eqs. (38). For such small bias voltages, we are able to reach Coulomb interaction strengths of U ≤ 2.5 Γ, with the data sets still converging. For larger U it would be necessary to take memory times Γt K > 2.5 and thus larger K into account. For the noninteracting case, we find a peak at where the single level is within the transport window. The peak height is, for the chosen parameters, at dI/dV ≈ 1.81 e 2 /h. The deviation from the maximally possible value of 2 e 2 /h is due to finite temperature. As 0 moves away from 0, the dot's energy level is pushed out of the transport window and the conductance drops significantly. When increasing the interaction strength U , one would expect a level splitting for single and double occupation of the quantum dot, and as a result two peaks at 0 = ±U/2 in the dI/dV curves. We do not reach high enough values of U to clearly resolve this peak splitting, but we see that the central peak becomes significantly broader. In addition, at 0 = 0 the conductance is suppressed with increasing U , with the maximal value for U = 2.5 Γ being dI/dV = 1.623 ± 0.002 e 2 /h. Next, we turn to the nonlinear-response regime. In Fig. 7, we show the differential conductance for a large bias voltage of eV = 3 Γ. The other parameters are as in Fig. 6. Due to the large bias voltage, the 0 dependence of the conductance is more complex. For U = 0, there are two peaks at 0 ≈ ±eV /2 = ±1.5 Γ, reflecting the two resonance conditions of the dot level matching the Fermi level of the left and right electrode, respectively. The peaks are not perfectly centered around ±eV /2 due to the finite width of the peaks. Since the two peaks overlap, their maxima get pulled towards each other.
At finite Coulomb interaction, there are, in principle, four resonance conditions, given by 0 ≈ ±eV /2 ± U/2. This explains that with increasing the Coulomb interaction U , the two peaks are pushed away from each other. Simultaneously, the peak heights decrease significantly. This is a consequence of the reduced overlap of the peaks as they away from each other.
The most remarkable feature, however, is the appearance of a third peak around 0 = 0 for U > 2 Γ. This is due to the fact that here U ≈ eV , and therefore at 0 = 0 the singly occupied state is in resonance with the right lead and the doubly occupied state is in resonance with the left lead. Since a noninteracting-electron picture only predicts two peaks, the appearance of additional peaks is clear indication of Coulomb interaction. However, a third (and ultimately a fourth) peak can only be resolved for sufficiently large values of U . This, we could not achieve by using the finite-time implementation of ISPI, but the TraSPI formulation now allows us to enter this regime.
To address the crossover from the linear-to the nonlinear-response regime, we discuss in Fig. 8 the differential conductance as a function of the bias voltage at 0 = Γ. The temperature is again k B T = 0.2 Γ and B = 0. For vanishing Coulomb interaction, we find, as expected, two peaks located at eV ≈ ± 0 = ±Γ. With increasing Coulomb interaction U , the differential conductance increases in the linear-response regime but decreases in the nonlinear-response regime, in accordance to our findings in Figs. 6 and 7.
We emphasize that the system under consideration is particle-hole symmetric, such that the differential conductance dI/dV is an even function of both 0 and eV . Both symmetries are fulfilled perfectly by the TraSPI formalism up to numerical accuracy.

B. Occupation number
We now address the average number of electrons on the quantum dot. For low-lying 0 , the quantum dot is occupied with two electrons. With increasing 0 , the occupation number is reduced and ultimately becomes empty.
In Fig. 9, we show the derivative of the occupation number as a function of gate voltage in the nonlinearresponse regime. We choose the same parameters as in   7. Similarly as the differential conductance is better suited to resolve detailed structures, the derivative of the occupation number shows more structure than the occupation number itself. In that respect, Fig. 7 for the current and Fig. 9 for the occupation number are analogous to each other. The peaks of dN/dV are at the same positions as the peaks of dI/dV , indicating that a large change of the current is accompanied with a large change of the occupation number. With increasing U , the peaks move away from each other, and the TraSPI formalism is able to accurately describe the full region in between. The numerical calculations of the conductance and the occupation number are independent of each other. The information contained in these two quantities is partially the same but in most transport regimes differes from each other in detail. There is, however, one limit in which they carry identical information. This occurs at zero temperature, T = 0, and vanishing bias voltage V = 0. In this case, the electrons transversing the quantum dot scatter only elastically, such that Friedel's sum rule can be applied, which leads to the Langreth formula [52,53] This remarkable result is valid for any value of the interaction strength U , covering all regimes from noninteracting electrons to strong correlations, e.g., in the Kondo regime. The Langreth formula is, in general, violated for any approximation scheme that only includes a certain class of transport contributions. However, for a numerically exact treatment, which includes TraSPI, the Langreth formula serves as a consistency check to assess the quality of the method. In Fig. 10, we compare the linear conductance calcu- lated in two different ways: once directly and once via the calculation of the occupation number and Langreth's formula Eq. (42). As stated above, Langreth's formula holds at zero temperature. This regime can be accessed by TraSPI only away from resonance, whereas close to the resonance condition 0 = 0, numerical convergence is too slow. Therefore, perform our calculations at finite temperature, choosing k B T = 0.2, in accordance with the previous figures. We find very good agreement away from resonance, where the influence of finite temperature on the conductance is small. The deviation of the Langreth result Eq. (42) from the direct calculation at resonance 0 is fully understood as a finite-teperature effect. In conclusion, Fig. 10 gives us strong confidence in the quality of TraSPI as a numerically-exact method.

V. CONCLUSIONS
In the literature, there is a phlethora of methods to describe quantum transport through nanostructures, which all have their advantages and disadvantages in different regimes. Some are restricted to the linear-response regime while others can cover strong nonequilibrium situations. In scenarios with a clear hierarchy of the involved parameters, a perturbation expansion in one of them can be used. However, in real experiments very often many of these parameters characterizing, e.g., temperature, tunnel-coupling strength, Coulomb interaction as well as gate and bias voltage, are of the same order of magnitude. Then, numerically exact methods are desirable.
In this paper, we presented TraSPI as such a numerically exact method. It is based on an iterative sum-mation of path integrals, referred to as ISPI. The virtue of ISPI (and thus also TraSPI) is that naturally takes into account all orders in tunneling of electrons between quantum dot and leads, allows for arbitrary bias voltages that drive the system out of equilibrium, is not restricted to either low or high temperature, and is able to include finite Coulomb interaction. While the method is numerically exact, stronger correlations increase the convergence time, such that IPSI (and TraSPI) are best suited for small to intermediate strengths of the Coulomb interaction. In previous applications of the ISPI method [14,15], we concentrated on spin-dependent phenomena, which show an interaction dependence already for moderate Coulomb interaction strengths. To increase the range of applicability towards stronger Coulomb interaction strengths, the efficiency of the method needs to be increased.
This we do in the present paper by mapping the ISPI scheme to a transfer-matrix approach, which results in TraSPI. The major virtue of involving transfer matrices is that the stationary limit is implemented by construction. This avoids the numerically-costly extrapolation of the results of the finite-time formulation as done in ISPI. In addition, the use of transfer matrices allow for further improvements that enhance the efficiency of the method. Numerical effort is reduced by optimally choosing the position of the measurement in time. Furthermore, to minimize numerical errors, it is advantageous to analytically implement derivatives, e.g., of the current with respect to bias voltage to get the differential conductance, instead of performing a numerical derivative of the current. And finally, we reduce numerical errors by making use of the possibility to calculate the noninteracting limit analytically and to numerically calculate only the difference between the interacting and the noninteracting case.
To illustrate the performance of the TraSPI method, we analyzed the differential conductance through a single-level quantum dot in both the linear and nonlinear regime. We were able to reach values of the Coulombinteraction strength that are sufficient to resolve in the non-linear response regime a third conductance peak instead of only two peaks that are expected for noninteracting electrons. Finally, we were able to perform a quality check of TraSPI by demonstrating that Langreth's formula, which connects the zero-temperature, linear conductance with the dot's occupation, is fulfilled in the regime of its applicability range. Therefore, we are confident that the TraSPI formulation enables us to address systems, transport regimes and effects that could not be covered by previous methods.

VI. ACKNOWLEDGEMENTS
We thank S. Weiss for fruitful discussions on the ISPI scheme. Financial funding of the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under project 278162697 -SFB 1242 is acknowledged.