Auxiliary master equation approach within stochastic wave functions: Application to the Interacting Resonant Level Model

TL;DR

This paper advances the auxiliary master equation approach (AMEA) by integrating stochastic wave functions, enabling faster and more accurate simulations of the Interacting Resonant Level Model in and out of equilibrium.

Contribution

The paper introduces a stochastic wave functions implementation within AMEA, significantly reducing computation time and improving accuracy for many-body quantum system simulations.

Findings

Enhanced simulation speed with stochastic wave functions

Accurate current-voltage and correlation function results

Effective extrapolation to infinite auxiliary bath sites

Abstract

We present further developments of the auxiliary master equation approach (AMEA), a numerical method to simulate many-body quantum systems in as well as out of equilibrium, and apply it to the Interacting Resonant Level Model (IRLM) to benchmark the new developments. In particular, our results are obtained by employing the stochastic wave functions (SWF) method to solve the auxiliary open quantum system arising within AMEA. This development allows to reach extremely low wall-times for the calculation of correlation functions with respect to previous implementations of AMEA. An additional significant improvement is obtained by extrapolating a series of results obtained by increasing the number of auxiliary bath sites, $N_B$, used within the auxiliary open quantum system formally to the limit of $ N_B \rightarrow \infty$. Results for the current-voltage characteristics and for equilibrium correlation functions are compared with the one obtained by exact and matrix-product based approaches.

Figures

Equations

$$H=H_{\text{imp}}+H_{\text{Baths}}+H_{\text{Hyb}}.$$

$$\mathcal{L}\rho$$

$${\mathcal{L}}_D\rho$$

$${\mathcal{L}}_{\alpha }\rho$$

$$+$$

$$\frac{d}{dt}\rho (t)=\mathcal{L}\rho (t).$$

$${\mathbf{G}}_{\text{imp}}^{-1}(\omega )={\mathbf{g}}_{0,\text{imp}}^{-1}(\omega )-{\mathbf{\Delta }}_{\text{ph}}(\omega )-\mathbf{\Sigma }(\omega ),$$

$${\chi }^2({\mathbf{x}}_{\alpha })\equiv {\chi }_{\alpha }^2$$

$$\left|{\Delta }_{\mathrm{ph}}^{(\alpha )}-{\Delta }_{\mathrm{aux}}^{(\alpha )}\right|$$

$$H=({\vec{h}}_1,..,{\vec{h}}_{{\text{rank}}_H}),\Gamma =H{H}^{†},$$

$$|\psi ⟩=\sum _{\lambda }\psi (\lambda )|\lambda ⟩,$$

$$D\psi D{\psi }^{\ast }\equiv \prod _{\lambda }\frac{i}{2}d\psi (\lambda )d{\psi }^{\ast }(\lambda ),$$

$$⟨A(t)⟩=\int D\psi D{\psi }^{\ast }⟨\psi |A|\psi ⟩P[\psi (\lambda ),t].$$

$$⟨A(t)⟩=\frac{1}{n}\sum _i^n\frac{⟨{\psi }_i(t)|A|{\psi }_i(t)⟩}{||{\psi }_i(t)|{|}^2}$$

$$G_{BA}(t,t^{\prime })=⟨\psi (t_0)|B(t)A(t^{\prime })|\psi (t_0)⟩$$

$$\mathbf{\Theta }(t^{\prime })=\left(\begin{array}{c}\psi (t^{\prime })\\ \varphi (t^{\prime })\end{array}\right)$$

$$G_{BA}(t,t^{\prime })=\frac{1}{n}\sum _i^n\frac{||\mathbf{\Theta }(t^{\prime })|{|}^2}{||{\mathbf{\Theta }}_i(t)|{|}^2}⟨{\psi }_i(t)|B|{\varphi }_i(t)⟩$$

$$\Delta A(\chi )=A_{\text{ph}}-A_{\text{aux}}(\chi )=k_A\chi +O({\chi }^2)$$

$$H_{\text{IRLM}}=$$

$$H_{\text{L}}=$$

$$H_{\text{R}}=$$

$$H_{\text{dot}}=$$

$$+U\sum _{r=\pm 1}\left(c_r^{†}{c}_r-\frac{1}{2}\right)\left(c_0^{†}{c}_0-\frac{1}{2}\right),$$

$$I(V)={\frac{V}{2\pi }}_2{F}_3\left[\left\{\frac{1}{4},\frac{3}{4},1\right\},\left\{\frac{5}{6},\frac{7}{6}\right\};-{\left(\frac{V}{V_c}\right)}^6\right],$$

$$H_{\text{imp}}$$

$$H_{{\text{B}}_L}$$

$$H_{{\text{B}}_R}$$

$$H_{\text{Hyb}}$$

$$I_{i,i+1}$$

$$=2E_{i+1,i}\Im ⟨c_{i+1}^{†}{c}_i⟩,$$

$$\mathcal{L}={\mathcal{L}}_H+{\mathcal{L}}_D.$$

$${\mathcal{L}}_H\tilde{\varrho }=-i[H,\tilde{\varrho }],$$

$$\begin{array}{rl}{\mathcal{L}}_D\tilde{\varrho }& =2\sum _{ij}{\Gamma }_{ij}^{(1)}\left(\pm c_j^{}\tilde{\varrho }{c}_i^{†}-\frac{1}{2}\left\{\tilde{\varrho },c_i^{†}{c}_j^{}\right\}\right)\\ & +2\sum _{ij}{\Gamma }_{ij}^{(2)}\left(\pm c_i^{†}\tilde{\varrho }{c}_j^{}-\frac{1}{2}\left\{\tilde{\varrho },c_j^{}{c}_i^{†}\right\}\right).\end{array}$$

$$2{\Gamma }_{ij}^{(\beta )}=\sum _k{U}_{ik}^{(\beta )}{\gamma }_k^{(\beta )}{U}_{jk}^{(\beta {)}^{\ast }},$$

$${\mathcal{L}}_D\rho$$

Full text

I Introduction

Quantum impurity models have a long history in many-body quantum mechanics. Some prominent examples are, the Single Impurity Anderson Modelande.61 (SIAM), the (Anderson-) Holstein Modelho.59 ,the Kondo Modelkond.64 and the Interacting Resonant Level Modelwi.fi.78 (IRLM). They feature interesting, unconventional physics such as the Kondo effecthews or negative differential conductancebo.sa.08 and allow for experimental realizations in terms of quantum dotsgo.go.98 . Besides this, the solution of quantum impurity problems alone constitutes already a crucial task in Dynamical Mean Field Theoryme.vo.89 .

Over the last decade, there has been increasing interest in quantum impurities out of equilibrium and the development of numerical methods which are able to accurately simulate such systems poses a great challenge for contemporary condensed matter theory. Existing methods ec.ha.09 include, iterated perturbation theory sc.mo.02u , numerical renormalization group jo.fr.08 , real time quantum monte carlo (QMC) ec.ko.09 ; ec.ko.10 , noncrossing approximation and beyond okam.08 ; ar.ko.12 , or imaginary-time QMC supplemented by a double analytical continuation han.07 ; ha.he.07 ; di.we.10 ; ar.we.13 , scattering-states approaches me.an.06 ; ande.08 , perturbative and renormalization group(RG) methods me.wi.92 ; sc.sc.94 ; ro.pa.05 ; scho.09 , time-dependent density-matrix RG (DMRG) and related tensor-network approaches wh.fe.04 ; da.ko.04 ; pr.zn.09 , numerical RG an.sc.05 , flow equation kehr.05 , functional RG ge.pr.07 ; ja.me.07 , dual fermions ju.li.12 ; ch.co.18u . A method developed over the last years is the so-called auxiliary master equation approachar.kn.13 ; do.nu.14 ; do.so.17 (AMEA). The advantage of this approach is that, in contrast to approaches which simulate a closed Hamiltonian system, it allows to directly address the steady state. Also time-dependent correlation functions can be readily evaluated starting from the steady state or any arbitrary initial condition. AMEA was successfully used as impurity solver within steady state non-equilibrium DMFTti.do.15 ; do.ti.16 ; ti.do.17 ; ti.so.18 as well as to calculate highly accurate spectral functions of the SIAM under the influence of a bias voltagedo.ga.15 ; fu.do.18 .

AMEA is based upon mapping the physical system to an auxiliary open quantum system of Lindblad form. The dynamics of the resulting auxiliary system is described by the density matrix and is solved by numerical means. In previous works the Lindblad system was solved by using the so called super-fermion (SF) representationdz.ko.11 , which formulates the super operator problem in terms of a standard operator problem with twice as many sites. The operator problem was than solved by standard numerical many-body techniques such as Krylov-space methodskn.ar.11.ec ; do.nu.14 (ED) or time evolution with Matrix Product Statesdo.ga.15 (MPS).

In this work, we implement an alternative solution strategy which does not rely on the SF representation, namely Stochastic Wave functionsda.ca.92 ; br.ka.97 ; br.ka.98 (SWF). The new method is statistical in nature and most notably highly parallelizable. This makes it a very promising candidate to exploit the multi-core architecture of (future-) cluster facilities. In addition, we introduce the notion of finite-size scaling within AMEA and report on progress regarding the optimization problem arising when mapping to the auxiliary system.

To test and benchmark the new developments, we apply AMEA to the IRLM where we can compare to existing literature. This work is structured as follows.

We begin by describing the technical aspects in Sec. II-IV and present the results for the IRLM in Sec. V. In more detail, in Sec. II we outline AMEA for spinless one dimensional systems. Sec. III, and the corresponding appendix Sec. A , is devoted to the description of the SWF algorithm, the finite-size scaling is introduced in Sec. IV. In Sec. V we apply AMEA within SWF to the IRLM and test the finite-size scaling scheme as well as the capability to compute correlation functions against the literature. Finally, we present our conclusion together with a summary and outlook in Sec. VI

II Auxiliary Master Equation Approach

We briefly review AMEA to deal with fermionic impurity problems. We consider a generic interacting region - the impurity - of size $N_{\text{imp}}$ connected to a left and right bath of non-interacting fermions. Accordingly, we write the Hamiltonian as

[TABLE]

Here, $H_{\text{imp}}$ describes the interacting region, $H_{\text{Baths}}=\sum_{\alpha=L/R}H_{\text{B}_{\alpha}}$ corresponds to the leftover reservoirs and $H_{\text{Hyb}}$ contains the hopping terms connecting the baths to the impurity. In the following we will assume that an individual bath is connected only to a single site of the impurity.

The idea of AMEA is now to model the physical situation by an auxiliary open quantum system described by the Lindblad equation. It consists of the impurity and additional bath sites to approximate the action of the leftover Hamiltonian on the interacting region.

In more detail, the Lindblad super-operator (Liouvillian) defining the dynamics of the open quantum system of size $L=N_{\text{imp}}+2N_{B}$ reads111for simplicity we neglect spin

[TABLE]

where $\rho$ is the density matrix of the Lindblad system. The Liouvillian of the dissipative bath sites is given by

[TABLE]

where $\alpha$ denotes the left/right reservoir222Here, it is worth noting that the matrices $E^{(}\alpha),\Gamma^{(\alpha),1},\Gamma^{(\alpha),2}$ are only non-zero in the part of the system which describes the corresponding bath. and $c_{i}^{(\dagger)}$ are the creation (annihilation) operators of a fermion on site $i$ of the open quantum system. The time-evolution of the system is described by the Lindblad equation,

[TABLE]

For the steady state of the original system, Eq.15, the Dyson equation for the interacting region in the formulation of Keldysh Green’s functions reads,

[TABLE]

where all objects have a matrix structure in the physical sites and Keldysh-space. In Eq. 6 $g_{0,\text{imp}}$ is the Green’s Function (GF) of the interacting region when isolated from the baths and without interaction, $\Sigma$ is the unknown selfenergy, holding all information about the interaction and $\Delta$ is the so-called hybridization describing the effect of $H_{\text{Baths}}+H_{\text{Hyb}}$ on the impurity. For the case considered here, the hybridization has the spacial structure $\text{diag}(\Delta_{\text{ph}}^{(L)},0,...,0,\Delta_{\text{ph}}^{(R)})$.

The mapping from the physical to the auxiliary system is performed by fitting the parameters $E^{\alpha},\Gamma^{\alpha,(1,2)}$ in Eq. 4 such that the hybridization in the auxiliary system approximates the physical hybridization as close as possible, $\Delta_{\text{aux}}^{(\alpha)}\approx\Delta_{\text{ph}}^{(\alpha)}$, and this is the only approximation made within AMEA. The accuracy of the mapping can than be systematically improved by increasing the number of auxiliary bath sites $N_{B}$ and it becomes formally exact in the limit of $N_{B}\rightarrow\infty$. Once the mapping is performed, one can solve the auxiliary system by some appropriate numerical method and evaluate observables belonging to the impurity. Their accuracy in describing the corresponding exact quantities will be directly related to the difference between $\Delta_{\text{aux}}^{(\alpha)}$ and $\Delta_{\text{ph}}^{(\alpha)}$.

II.1 Mapping to the auxiliary System

Here, we briefly want to summarize the mapping procedure and mention key points that we need for the present work. For a thorough discussion of the mapping and technical details we refer to our previous workdo.so.17 . The mapping is performed for each individual bath $\alpha$ by minimizing a suitable cost function

[TABLE]

Here we have introduced a parameter vector $\boldsymbol{x}_{\alpha}$ that parametrizes the matrices $E^{\alpha},\Gamma^{\alpha,(1,2)}$ in (4), from which one evaluates the auxiliary hybridization $\Delta_{\mathrm{aux}}^{(\alpha)}$. It is important to note that the precise form of the cost function is very flexible and may be chosen differently for different physical situations. One important property of the mapping is that the cost function decreases exponentially with the number of fit parameters, $-\log\chi_{\alpha}\propto\text{dim}(\boldsymbol{x_{\alpha}})$, which typically leads to a rapid increase of accuracy when the number of bath sites $N_{B}$ is increased.

In previous works, Eq. 7 was minimized via a parallel tempering (PT) algorithm which is appropriate to find the global minimum. However, it should be noted that within AMEA it is not strictly necessary to find the global optimum333although desirable as it gives the best approximation for a given system size.. In general, the fit struggles to resolve sharp features such as band-edges in the retarded component or the fermi-jumps in the Keldysh component at zero temperature. Therefore, $T=0$ can not be reached exactly in practice and the auxiliary system always has some non-zero effective temperature.

II.1.1 Developments of the fit

With increasing dimensionality of the fitting problem, the PT algorithm gets computationally prohibitive and it is not able to find even good local minima anymore for444This corresponds to the case $N_{B}>6$ when allowing for the most general Lindblad couplings. $\text{dim}(\boldsymbol{x})=2N_{B}(N_{B}+1)\gtrsim 80$. Good minima should be such that they display an exponential decrease in the cost function when the number of bath sites is increased. To obtain good enough minima for $N_{B}=7,8$, we use the fact, which we observed empirically, that the $\Gamma$ matrices of obtained minima typically have very low rank. Utilizing a variable rank parametrization in terms of a corresponding matrix $H$

[TABLE]

where $\vec{h}_{i}$ denote column vectors of length $L$. Note that the maximal useful rank typically increases with the system size555For example, we observed that for a system with $N_{B}=9$ bath sites, increasing $\text{rank}_{H}>4$ was not fruitful in terms of the cost-function. With this procedure, we have reduced the dimensionality of the parameter vector to $\text{dim}(\boldsymbol{x})=2N_{B}(\text{rank}_{H}+1)$ extending the applicability of the PT algorithm to about $N_{B}=8$. To achieve an exponential decrease in the cost function for even more bath sites we have adopded an optimization algorithm which makes use of the gradient of the cost function, which can be evaluated directly. This information is not used in the PT algorithm. Suitable gradient-based approaches can be found in the area of machine learning, which provides algorithms tailored to find local minima in very high-dimensional problems utilizing variants of steepest descent. Here, we employ the ADAM ADAM optimizer as implemented in the python library tensorflowtensorflow .

Steepest descent approaches are obviously very sensitive to the starting point. In our case, it has proven to be very effective to first find the solution for a small auxiliary system (small $N_{B}$) and consequently add bath sites until the required $N_{B}$ is reached. For a fixed $N_{B}$ we start with the result of the previous system size and increase the rank stepwise until no significant decrease in the cost-function is observed. In addition to beeing applicable for larger $N_{B}$, the ADAM routine is faster than PT for a given $N_{B}$.666This is because the PT algorithm tries to explore the total phase space, wheres as ADAM only follows a certain path..

III Solution of the Linblad system with Stochastic Wave functions

The auxiliary open system is still correlated but due to its finite size can be addressed by numerical techniques. One route is to make use of the so-called Super-Fermion (SF) representation dz.ko.11 , which maps a super-operator problem to a standard, albeit non-hermitian, operator problem. The drawback of this approach is that the resulting SF problem is formulated on twice as many effective sites leading to a rapid increase in the numerical complexity. In previous works employing AMEA we have successfully used the SF representation together with established many-body techniques such as Krylov-space methodsar.kn.13 ; do.nu.14 or MPSdo.ga.15 to solve for steady state properties. A completely different route is to use Stochastic Wave functionsda.ca.92 ; br.ka.97 ; br.ka.98 (SWF), also referred to as “quantum jumps”, to solve the auxiliary many-body problem. The method is based on the stochastic nature of the Lindblad problem and is formulated in terms of wave functions instead of a density matrix and thus circumvents the need to square the Hilbert space. In the following, we will only give a brief introduction to the SWF method and focus more on a practical prescription to simulate the many-body Lindblad system arising within AMEA. For more details, mathematical definitions and background we refer to the literatureda.ca.92 ; br.ka.97 ; br.ka.98 .

The density operator $\rho(t)$ can be mapped onto a probability distribution $P[\psi(\lambda),t]$ for the quantum mechanical (many-body) wave function777To be consistent with quantum mechanics, $P[\psi,t]$ must not depend on the phase of the wavefunction and it is only non-vanishing for normalized states.

[TABLE]

where $\lambda$ indexes a complete set of (many-body) basis states888More generally, $\lambda$ indexes a complete set of quantum numbers. With the Hilbert space volume element,

[TABLE]

defining the needed probability measure999$P[\psi(\lambda),t]D\psi D\psi^{*}$ can then be interpreted as the probability to find the system within the volume element $D\psi D\psi^{*}$ around the state $\psi(\lambda)$ at time $t$., the expectation value of an observable can then be formally expressed as

[TABLE]

In short, instead of dealing with an evolution equation for the density matrix, one formulates a stochastic process on the Hilbert space.

For the specific case of a Lindblad system, the process is simulated according to a stochastic differential equation101010The mapping to a stochastic differential equation is possible also in a more general context, but the exact form of the latter is only known in special cases. leading to the algorithm presented in Fig. 1.

In this algorithm, a state vector $|\psi\rangle$ is evolved in time according to an effective, but non-hermitian Hamiltonian, $H_{\mathrm{eff}}$. $H_{\mathrm{eff}}$ comprises the Hamiltonian $H_{\text{imp}}$ as well as the particle-number conserving terms from the part describing the L/R baths, i.e. the terms proportional to $E^{(\alpha)}_{ij}$ as well as the terms containing the anticommutators in Eq. (4). This deterministic time evolution is interrupted by stochastic jump processes to different particle sectors, mediated by jump operators $L_{k}^{(\beta)}$, see appendix A for details. Observables are determined as the average over expectation values in independent realizations of $|\psi\rangle$. Such a stochastic unraveling of the Lindblad equation into a pure state description, as described above, only works for proper density operators $\rho$. When evaluating a Green’s function, one needs the stochastic time evolutions operators obtained by applying an operator $A$ to $\rho$. In order to compute two-time correlation functions,

[TABLE]

we follow the approach outlined in Ref. br.ka.97, and consider the stochastic time evolution of a doubled Hilbert space resulting in the algorithm in Fig. 2.

Here, a state vector $|\psi\rangle$ is evolved in time together with a corresponding vector $A|\psi\rangle$. A Green’s function is then proportional to the stochastic sample of off diagonal matrix elements of the second operator $B$, see Eq. 13. Notice that for single fermion Green’s functions, $A$ is a fermionic creation/annihilation operator. In that case one has to use the negative sign in front of the jump term for the lower part of the doubled Hilbert space, cf. Eq. 25, see appendix B in  sc.go.16 . Notice that generalizing the doubled Hilbert space to a multiple Hilbert space allows to sample different correlation functions at once, see Appendix Sec. A.2.2.

The SWF algorithm requires a routine which is able to time-evolve an initial vector with a non-hermitian generator for some (arbitrarily-) small time $dt$111111Although it is in favor of the algorithm if the routine is able to time evolve directly between consecutive jumps, that is for times $dt\approx\tau_{\text{jump}}<1/J$. In the present work we use the so-called Arnoldi algorithmkn.ar.11.ec for the time evolution which is the Lanczos method generalized to the non-hermitian case. For more details we refer to the Appendix Sec. A.

IV Extrapolation of Observables to the limit of vanishing cost function

As illustrated above, AMEA is a method which can be systematically improved by increasing the number of bath sites $N_{B}$ leading to an exponential decrease in the cost function, $\chi=\sum_{\alpha}\chi_{\alpha}$, which is a measure of the overall accuracy. Clearly, the best approximation for some quantity of interest for given $N_{B}$ is obtained within the auxiliary system with the smallest $\chi$. To improve on these results one can think of numerically extrapolating the results to the $\chi\to 0$ limit. This is equivalent to a scaling to the limit of an infinite number auxiliary bath sites $N_{B}\to\infty$. However, since the accuracy is directly related to $\chi$ rather than $N_{B}$, it is more convenient to use $\chi$ as an extrapolating parameter. For a given observable $A$ of interest we can assume for its deviation from its exact (physical) value

[TABLE]

with some constant of proportionality $k_{A}$. This suggests that given a series of value pairs $\{\chi_{i},A(\chi_{i})\}$ one can obtain an approximation to $A_{\text{aux}}(\chi=0)$ by performing a linear fit in the $(\chi,A)$ plane. Within AMEA a series of value pairs $\{\chi(N_{B}),A\}$ is naturally generated by the different possible auxiliary system sizes.

We want to emphasize that the extrapolation scheme presented here is not able to give a consistent error estimate of the extrapolated value as the uncertainty of the individual data points is unknown and not statistically distributes 121212Here, one has to distinguish between a purely statistical error stemming from the solution of the Lindblad system within stochastic wave functions, which is known and negligible, and the systematic error introduced by the mapping to the auxiliary system which is unknown. Further, the role of the higher order terms in Eq. 14 introduces another source of unknown error. To get a grip on the error due to the AMEA mapping, one could perform the extrapolation in some limit where the true value in the physical system is known, for example at zero interaction strength or for some other parameters where the value is known from the literature. One could then use the deviation from the extrapolation fit as approximation to the error of a data point. Since there is a lot of freedom in obtaining this error estimates - and it will thus be very situation dependent- we will not pursue this further in the current work where we are interested in an unbiased benchmark of the extrapolation scheme.. Nevertheless, this scheme provides a significant improvement, for example in the current, as can be seen in Fig. 5.

V Application to the Interacting Resonant Level Model

The IRLMwi.fi.78 is a commonly used non-equilibrium impurity model of spinless fermions. It features an impurity site connected to two semi-infinite tight-binding chains together with an interaction term coupling the particle densities of the impurity site to the neighboring chain sites, see Fig. 3. The Hamiltonian is defined as

[TABLE]

where $c_{r}^{\dagger}/c_{r}$ denote the fermionic creation/annihilation operators at site $r$. Here, $H_{\text{L/R}}$ describe the semi-infinite tight-binding chains of bandwidth $W=4J$ and $H_{\text{dot}}$ introduces the hopping to the impurity as well as the interaction term. A non-equilibrium steady state situation is induced in the system via an applied bias voltage $V$ simulated by shifting the chemical potentials of the leads symmetrically, that is $\mu_{l}=-\mu_{r}=\frac{V}{2}$. We use $J$ as unit of energy and work in units where $\hbar=e=k_{B}=1$.

The IRLM is known to be integrable wi.fi.78 and becomes equivalent to the continuum model in the so-called scaling regime where the bandwidth becomes the dominant energy scale in the system. Most notably, there is a closed form expression for the steady state current as a function of the bias voltage bo.sa.08 ; ca.ba.11 for the special value of the interaction $U=2$,

[TABLE]

with $V_{c}=r(J^{\prime})^{\frac{4}{3}}$ and $r\approx 3.2$131313In more detail, $V_{c}=\frac{\sqrt{3}}{4^{2/3}}\frac{4\sqrt{\pi}\Gamma(2/3)}{\Gamma(1/6)}T_{B}$ and $T_{B}=c(J^{\prime})^{\frac{4}{3}}$ with $c\approx 2.7$ from bo.sa.08 .. Here, ${}_{2}F_{3}(a,b;z)$ is the generalized hypergeometric function. The formula Eq. 16 is valid at zero temperature and in the scaling regime, where $V,J^{\prime},U\ll W$141414From previous works bi.mi.17 , we know that at $U=2$ one has to restrict to $J^{\prime}\lesssim 0.5$ and $V\lesssim 2$ to be in the scaling regime.. In this way, $I/V_{c}$ becomes a universal function of the scaled voltage $V/V_{c}$ alone and in particular does not depend on the hybridization strength $J^{\prime}$.

V.1 AMEA for the IRLM

In the IRLM, the interaction lives on the contact links to the leads and, therefore, the interacting region comprises the sites $r=\{-1,0,1\}$ which corresponds to having

[TABLE]

as indicated in Fig. 3. Since $H_{\text{B}_{L/R}}$ describe semi-infinite tight-binding chains in equilibrium, $\Delta_{\text{ph,L/R}}$ represent baths with a semicircular density of states with a bandwidth of $W=4$ and an electronic distribution function given by the Fermi-function. Within AMEA, a given parameter set $E^{\alpha},\Gamma^{\alpha,(1,2)}$ fixes both the density of states as well as the distribution function of the corresponding bath. Since the Hamiltonian Eq. 15 is particle-hole symmetric, it suffices to perform the fit only for one of the two baths, e.g. the left ones, and obtain the parameters of the right bath by particle-hole transformation. Thus, also the cost function for the left and right bath will be equal for a given bias voltage, $\chi_{L}=\chi_{R}$. To illustrate the mapping, we show in Fig. 4 two examples for such a fit with $L=13$ ($N_{B}=6$) and $L=19$ ($N_{B}=9$). Notis that the same fit can be used for any set of parameters in this model.

V.2 Extrapolation of the steady state Current

Since there are no free parameters in Eq.16 we can use this as a benchmark for our numerical approach and test the extrapolation scheme of Sec. IV. However, it should be noted that our results are obtained for $T=0.025$ while Eq. 16 is the result for zero temperature.

Given an auxiliary system of size $L$ we can evaluate the current over a physical bond $i$ in the auxiliary system151515In practice, we measure at all physical bonds and average accordingly. There is a tiny breaking of current conservation due to the numerics.

[TABLE]

where the parameters $E_{i,i+s}$ represent the hopping along the chain in the interacting region. In the following, we consider results obtained with $7\leq L\leq 19$. In Fig. 5 we plot the universal steady state current together with the corresponding data points obtained with AMEA for $J^{\prime}=0.5$ and $J^{\prime}=0.2$. Shown are the AMEA results for individual system sizes as well as the extrapolated current. For $J^{\prime}=0.2$ the hybridization strength $\Delta(\omega=0)=J^{\prime 2}=0.04$ becomes comparable to the temperature used in our calculations, $T=0.025$. Thus, for $J^{\prime}=0.2$ we disregard data points corresponding to small voltages $V\leq 0.6$ because they are significantly altered by the finite temperature present in the auxiliary system.161616At small voltages, the current is carried by states around the chemical potentials which are most affected by the finite temperature.

We see that the current improves significantly towards the analytic solution thanks to the extrapolation scheme. As discussed above, the analytic solution is only valid for not too large bias voltagesbi.mi.17 . Indeed, we see a systematically growing deviation between the analytic solution and the current from AMEA171717Of course, in the present case it is not strictly possible to distinguish between deviations coming from the finite temperature and ones originating from leaving the scaling regime. for voltages $V\gtrsim 2$, see the markers in Fig. 5

The inset in Fig. 5 shows an example extrapolation. As one would expect, the data points with bigger cost functions (smaller system sizes) show a stronger scattering from the linear fit than the more accurate points. While the points with low cost functions make for more confidence in the results, the accuracy of the extrapolated current does not suffer when the biggest system size, $L=19$, is excluded from the analysis. This suggests that when utilizing the extrapolation to zero cost function, it is probably not necessary to simulate the biggest system sizes within reach. Rather, one can check for a small fraction of points whether or not the - usually very cpu-time intensive - bigger system size(s) are worth calculating181818If error estimates are used, points at lower cost-functions will reduce the uncertainty in the final result..

V.3 Spectral function of the IRLM

In this section, we evaluate the steady-state single-particle Green’s function $G$ at the central impurity site. The calculation is carried out in the real time domain and we use the approach discussed in Sec. III, see also Sec. A.2. We use a step size of $dt=0.05$ and $10^{5}$ time steps to first reach the steady state at $t_{0}=5*10^{3}$. We have verified that expectation values of static observables don’t change after this time. Then we sample the Green’s function $G(t-t_{0})$ for later times beyond $t_{0}$ up to $t_{\text{end}}=t_{0}+6000dt$. This is sufficient, since here $G(t_{\text{end}}-t_{0})<10^{-6}$

Finally, we average $G$ over $O(10^{5})$ realizations and determine the spectral function by direct Fourier transform 191919We have tried to take the statistical error into account for the Fourier transform within the framework of Bayesian probability theory but it did not lead to more satisfactory results. It turns out that an accurate error estimate for a given frequency is simply given by linear error propagation in the numerical Fourier integration. All results presented in this section are obtained with an auxiliary system of size $L=13$. The corresponding hybridization function is shown in the lower panels of Fig. 4.

Like any non-equilibrium approach, AMEA is also applicable in equilibrium situations which is just the special case when $\mu_{l}=\mu_{r}=0$ allowing us to compare our results against the literature. In Fig. 6, we compare our results to the equilibrium density of states obtained by Braun and Schmitteckert via MPS br.sc.14 . For interaction strengths $U<2$ that are small compared to the bandwidth, we observe a very good agreement with the reference over the whole frequency range. At the self dual point $U=2$ we start to see small quantitative deviations of peak heights but still obtain an satisfactory agreement. When the interaction becomes comparable to the bandwidth, $U=3$, the deviations become significant and continue to grow as the interaction is increased (not shown). The reason for the growing deviations is that in the present AMEA mapping the region outside the bandwidth is not well reproduced, see also Fig. 4. While these states do not play a role as long as all energy scales in the system are small compared to the bandwidth, i.e. in the scaling regime, the details of the leads at higher energies become important when the interaction becomes comparable to the bandwidth. The latter does, however, not mean, that AMEA is not at all applicable in this parameter regime, rather one has to make sure that the region outside the bandwidth is also faithfully reproduced by the auxiliary system. This can be achieved by using a differently distributed cost function in the fit or by going to larger auxiliary system sizes.

V.4 Performance:

From a numerical point of view, the stochastic wave function (SWF) method has two main advantages. First, since one evolves wave functions there is no need to square the Hilbert space as when one deals with the density matrix. This means that one can use a twice as large $L$, and thus, achieve a much better accuracy. 202020This does not hold for approaches in which the system size is not a limitation, such as tensor network states, i.e. MPS, where the entanglement entropy encoded into the state limits the simulation.. Second, individual realizations of possible time evolutions are independent which means that the method is easily parallelizable. This makes SWF very suitable for future cluster facilities which thrive on highly parallel algorithms.

However, the prize to pay is a cpu time that is about twenty times longer than solving an auxiliary system with the same value of the cost-function by MPS212121We note that the size $L$ of the auxiliary system solved by MPS must be larger, compared to the present approach, to reach the same value of the cost-function. This is because in MPS, the matrices $E,\Gamma^{(1)},\Gamma^{(2)}$ have to be tridiagonally restricted which lowers the number of fit parameters available for a given $L$.. On the other hand, thanks to parallelization, the wall-time222222The time it takes before the result is known can obviously be made almost arbitrarily small. For example, the GF’s for $L=13$ in Fig.6 where averaged over about half a million realizations where a single one takes around one second. For comparison, the solution with the super-fermion plus ED approach for $L=13$ would be in the order of minutes.

VI Summary, Conclusion and Outlook

We reported on technical developments within the auxiliary master equation approach and applied it to the Interacting Resonant Level Model (IRLM) in and out of equilibrium to benchmark the new techniques. We successfully applied the Stochastic Wave Function (SWF) algorithm to determine the steady state properties of the auxiliary Linblad system. On the one hand the SWF algorithm is highly parallelizable allowing to reach very low wall-times. On the other hand, we found that in the current implementation of SWF+ED the total cpu-time for a spectral function is twenty times higher than in available alternatives for the solution of the auxiliary system introduced by AMEA. Further, we saw that an auxiliary system size of $L=13$ is enough to obtain reliable spectral information of the IRLM for interactions $U\lesssim W/2$.

We obtained a further significant improvement by extrapolating physical quantities, most notably the current, to the $N_{B}\to\infty$ limit. In fact, it turns out to be more effective to extrapolate linearly in the cost function $\chi$, which then would correspond to an exponential extrapolation in $N_{B}$. Such an extrapolation is able to improve the results significantly and possibly circumvents the need to go to larger system sizes.

In addition, we introduced a variable rank parametrization of the auxiliary Lindblad matrices which typically reduces the number of fitting parameters in the AMEA mapping. Employing the new parametrization together with an optimization routine from machine learning, we were able to maintain an exponential decrease of the cost function also for larger system sizes where the previously used parallel tempering algorithm failed

In the current work, we calculated spectral functions only in equilibrium, where we can compare to the literature, while the present method also allows to calculate spectral functions in the non-equilibrium situation. Since the current through a system can also be expressed in terms of GF’s, we can investigate how the negative differential conductance in the IRLM arises from the spectral properties out of equilibrium. However, this is beyond the scope of the present work and will be presented elsewhere.

Further improvement in accuracy and computational time could be possibly be achieved by linear prediction, in order to extrapolate the Green’s function to large times, and by combining MPS with the present SWF approach.

Acknowledgements.

We would like to thank Irakli Titdvinidze, Daniel Bauernfeind and Gerhard Dorn for fruitful discussions. A special mention goes to Franz Scherr who introduced us to the machine learning environment tensorflow and provided a first implementation for the AMEA mapping. We are grateful to Peter Schmidteckert for providing us with the reference data for the spectral functions. This work was partially supported by the Austrian Science Fund (FWF) within Projects P26508 and F41 (SFB ViCoM), as well as NaWi Graz. The calculations were partly performed on the dCluster and lCluster Graz as well as the VSC-3 cluster Vienna

Appendix A Technical details of the SWF algorithm

In order to present the SWF algorithm, we consider a general Lindblad system for a generalized “density-matrix” $\tilde{\varrho}=f(\{c^{(\dagger)}\})\rho$ where $f(\{c^{(\dagger)}\})$ denotes some function of fermionic operators,

[TABLE]

It is composed of a central region with Hamiltonian $H$ and the corresponding Liouvillian $\mathcal{L}_{H}$,

[TABLE]

and a dissipative part described by $\mathcal{L}_{D}$,

[TABLE]

Here, $i$ and $j$ run over all $L$ sites of the system and $\Gamma^{(1)/(2)}$ are $L\times L$ matrices. The minus sign in $Eq.\thinspace{}(\ref{eq:L_D})$ is valid, if $\tilde{\varrho}$ is odd in the number of fermion operators, i.e $\tilde{\varrho}=c^{(\dagger)}_{i}\rho$. This is the case with Green’s functions, where we need to propagate $c_{i}^{(\dagger)}\rho$.

In order to obtain the jump operators one has to diagonalize the matrices $\Gamma^{(\beta)}$, $\beta=1,2$,

[TABLE]

and end up with the eigen-decomposition of the dissipator,

[TABLE]

The anti-commutators in Eq. (25) are included into the effective, non-hermitian Hamiltonian 232323For this generalize $[H,\rho]$ to $H_{\mathrm{eff}}\rho-\rho H_{\mathrm{eff}}^{\dagger}$.

[TABLE]

With this Hamiltonian and the jump operators $L_{k}^{(\beta)}$, Eq. (26), one formulates the SWF algorithms in Sec. III, Fig. 1 and Fig. 2.

A.1 Jump-time search and Arnoldi

As mentioned in Sec. III we use the so-called Arnoldi algorithm kn.ar.11.ec for the time evolution. Arnoldi is a Krylov space method analogue to Lanczos but for non-hermitian Hamiltonians. For a given initial state, $|\psi_{0}\rangle$ and time interval $dt$, a Krylov space, spanned by $Q$, is generated by iteratively applying $H_{\mathrm{eff}}$ to the starting vector until a satisfactory approximation for the time evolution operator $e^{-iH_{\mathrm{eff}}dt}\approx Q^{\dagger}e^{-iH_{\mathrm{K}}dt}Q$ is found. For any given time $t$ up to the maximal time $dt$, the state and the corresponding norm needed for the SWF algorithm are given by

[TABLE]

where we have used the property that $Q|\psi_{0}\rangle$ is nothing else than the first Krylov vector and $QQ^{\dagger}=\mathds{1}$ is the identity242424Note that for non-hermitian problems $Q^{\dagger}Q\neq\mathds{1}$. We want to point out that by virtue of Eq. 31 the norm can be calculated within the Krylov space representation itself, which is typically of size $\text{dim}_{K}=O(10)$, without the need to use the transformation matrices $Q$ which are of dimension $\text{dim}_{Q}=\text{dim}_{F}\cdot\text{dim}_{K}$ where $\text{dim}_{F}$ is the dimension of the Hilbert space (many-body Fock space). Differentiating Eq. 31, yields

[TABLE]

which allows to determine the jump time $t_{j}$ in the SWF algorithm, satisfying $\|\psi(t_{j})\|^{2}-r_{j}=0$, by applying Newtons method.

A.2 Practical implementation for the steady state situation

Here we want comment on the practical implementation for the special case of steady state quantities.

A.2.1 Steady state observables

We start with the simpler case of sampling a steady state observable. A steady state expectation value is obtained like in a Monte Carlo (MC) simulation. We start with a random starting state and time evolve the system until it reaches the steady state, where the system is time-translational invariant (like the thermalization in a MC simulation). Once we are in the steady state, we start measuring the observable generating an autocorrelated time series from which an estimator of the expectation value can be obtained. As usual the time-series needs to be long enough to have overcome autocorrelations, which can be checked for example by a Binning plot.

For the present case we typically recorded $N_{m}=2^{18}$ measurements separated by a time $\Delta t=N_{t_{\text{skip}}}dt$ with a time-step $dt=0.05$ and $N_{t_{\text{skip}}}dt=16dt\approx 10\bar{t}_{j}$, where $\bar{t}_{j}$ is the average jump time. For thermalization we performed additionally $10\%$ of the total time evolution leading to $O(10^{5})$ thermalization time-steps. Parallelization can be achieved by running several individual walkers on a single cluster node, where each walker is bound to one core for instance.

A.2.2 Steady state single particle GF’s

To obtain steady state GF’s of the Lindblad system we follow Ref do.ga.15 . In short, it is best to calculate the lesser and greater steady state GF, defined by

[TABLE]

where $\langle\cdot\rangle_{\infty}=T_{r}\{\cdot\rho_{\infty}\}$ denotes the expectation value in the steady state. We sample the GF by first time-evolving into the steady state like above. Next, we apply the operator $c_{r}^{(\dagger)}$, construct the doubled Hilbert space, continue to time-evolve in the doubled Hilbert space and measure according to the SWF algorithm in the doubled Hilbert space.

As stated in the main text the time steps needed for GF’s is of $O(10^{3})$ and to reach the accuracy needed for smooth spectral functions, we had to average over $O(10^{5})$ realizations. Further, we perform $O(10^{5})$ time-steps to get into the steady state. For the performance in terms of cpu-time, it is crucial that the time steps into the steady state are done only for a small fraction of the realizations; the corresponding final states are saved252525In the present case the time-steps into the steady state make for about ten per cent of the total run time.. Another realization starts from a state obtained by time evolving such a saved state for some time $\Delta t\approx 100\bar{t}_{\text{j}}$, to make sure that individual realizations are independent to a very good approximation262626Only early times will be correlated as the realizations gain in independence through the jumps in the time-evolution. One can test for autocorrelations when considering the different realizations for a specific (early-) time step as a time series and apply autocorrelation analysis..

Multistates:

One can sample multiple correlation functions, $G_{B_{i}A_{i}}(t,t^{\prime})$, together when generalizing the doubled Hilbert space to a multiple Hilbert space. For this, generalize

[TABLE]

with the excited states $\phi_{i}=A_{i}\mathinner{|{\psi}\rangle}$. For instance, this allows to sample the lesser and greater GF together in a tripled Hilbert space or multiple components of a cluster GF. The advantage is that $|\psi\rangle$ is only time evolved ones, where as in the individual approach, with only a doubled Hilbert space, $|\psi\rangle$ is time evolved $n$-times.

Destroyed states in the multiple Hilbert space

Here, we want to elaborate on the fact that part of the state may be destroyed when applying the SWF algorithm in the multiple Hilbert space. For simplicity, we consider in the following a doubled Hilbert space. Part of the state can get destroyed, when the system leaves the physical particle sectors through the application of a jump operator272727Naturally, this happens more often in smaller systems.. For instance, a state can get destroyed when the system is in the $N=L$ particle sector and a jump operator $L^{(2)}_{k}$ gets chosen that increases the particle number.

First, let us note that this cannot happen in the single Hilbert space since the corresponding weight $w_{\beta k}\propto\mathinner{\lvert\lvert{L_{k}^{(2)}\psi(t_{j})}\rvert\rvert^{2}}$ is zero and this jump operator will never be chosen.

The situation is different in the doubled Hilbert space when the two components of a state reside in different particle sectors. Too see this, let us consider the case of the greater GF. Here, if $|\psi\rangle$ is in sector $N$, $|\phi\rangle$ will always describe a state with $N+1$ particles, since the jump operator applied is the same for both components. If at some time $t_{\text{kill}}$, $|\phi\rangle$ is in the sector $L$, the weight for a jump operator that increases the particle number, $w_{\beta k}\propto\mathinner{\lvert\lvert{L_{k}^{(2)}\psi(t_{j})}\rvert\rvert^{2}}+\mathinner{\lvert\lvert{L_{k}^{(2)}\phi(t_{j})}\rvert\rvert^{2}}$, might be non-zero since the first part can be non-vanishing.

If part of the state is destroyed, all subsequent measurements in this specific realization of the time series for the GF will all be zero.

It is important to realize that this is the correct behavior. It exemplifies why the doubled Hilbert space is needed when calculating correlation functions and why it would be wrong to simply consider an independent time evolution for the excited state and the initial state separately. In fact, in the independent approach, any correlation between the initial state and the final state would be lost very quickly through the stochastic process and it is key that the two states always jump together, thereby mediating the correlation.