Optimized auxiliary representation of non-Markovian impurity problems by a Lindblad equation

In a different class of approaches, one tries to 'eliminate' the degrees of freedom of the reservoir and take into account its effects on the dynamics of the interacting central region [63–67]. One way to do this is by treating the coupling to the reservoir within a Lindblad equation [63–65]. In this way, the effect of the reservoir is to introduce nonunitary dynamics in the time dependence of the reduced density operator ${\rho }_{f}$ of the central region leading to the Lindblad equation, which is a linear, time-local equation for ${\rho }_{f}$ preserving its hermiticity, trace, and positivity. One important precondition for the validity of this mapping, however, is the Markovian assumption that the decay of correlations in the reservoir is much faster than typical time scales of the central region. As pointed out, e.g. in [63, 66], the approximations leading to the Markovian Lindblad master equation are justified provided the typical energy scale Ω of the reservoir is much larger than the reservoir–central region coupling. However, for a fermionic system, Ω can be estimated as $\min (W,\max (| \mu -\varepsilon | ,T))$, where W is the reservoir's bandwidth, and $\varepsilon$ is a typical single-particle energy of the central region. Therefore, even in the wide-band limit $W\to \infty$, the validity of the Markov approximation is limited either to high temperatures or to chemical potentials far away from the characteristic energies of the central region. As a matter of fact, the effect of a noninteracting reservoir with $W,| \mu | \to \infty$ (or $T\to \infty$ with finite $\mu /T$) can be exactly written in terms of a Lindblad equation. This can be easily deduced from the 'singular coupling' derivation of the Lindblad equation [63]. This is valid independently of the strength of the coupling between central region and reservoir. A nontrivial situation is obtained by introducing different reservoirs with different particle densities. The pleasant aspect of this limit is that the Lindblad parameters depend on the properties of the reservoir and of its coupling with the central region only, but not on the ones of the central region.

This is in contrast to the more standard weak-coupling Born–Markov version in which the Lindblad couplings (see, e.g. [63–65]) depend on the central region's properties. To illustrate this, consider a central region consisting of a single site with energy ${\varepsilon }_{f}$, i.e. with Hamiltonian

$$\begin{eqnarray}&&{H}_{f}={\varepsilon }_{f}\ {f}^{\dagger }f\end{eqnarray} \tag{ 1 }$$

(omitting spin) and reduced density matrix ${\rho }_{f}$. The part of the Lindblad operator ${{ \mathcal L }}_{b}$ describing the coupling to a noninteracting reservoir is given by

$$\begin{eqnarray}&&{{ \mathcal L }}_{b}\ {\rho }_{f}={{\rm{\Gamma }}}_{1}(2f{\rho }_{f}{f}^{\dagger }-\{{f}^{\dagger }f,{\rho }_{f}\})+{{\rm{\Gamma }}}_{2}(2{f}^{\dagger }{\rho }_{f}f-\{{{ff}}^{\dagger },{\rho }_{f}\})\end{eqnarray} \tag{ 2 }$$

with

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{1}={\rm{\Gamma }}(1-{f}_{{\rm{F}}}({\varepsilon }_{f}))\quad \quad {{\rm{\Gamma }}}_{2}={\rm{\Gamma }}\ {f}_{{\rm{F}}}({\varepsilon }_{f}).\end{eqnarray} \tag{ 3 }$$

Here, Γ is proportional to the reservoir's density of states at the energy ${\varepsilon }_{f}$, and f_F is the Fermi function which obviously contains the information on the chemical potential and temperature of the reservoir but also on the onsite energy ${\varepsilon }_{f}$ in the central region. This could be unsatisfactory since one would like to describe the effect of the reservoir in a form which is independent of the properties of the central region, especially when the latter consists of many coupled sites.

One possible way to eliminate the dependence of the Lindblad couplings on the parameters of the central region is to use an intermediate auxiliary buffer zone (mesoreservoir) between the Lindblad couplings and the central region (see, e.g. [68–71]). The buffer zone consists of isolated discrete sites (levels), each one coupled to a Markovian environment described by Lindblad operators with the same T and μ as given in equations (2) and (3). If the buffer zone is sufficiently large, i.e. if its levels are dense enough, then one can show that the buffer zone including Lindblad operators yields an accurate representation of the reservoir, which becomes exact in the limit of an infinite number of levels. Importantly, the parameters of this buffer zone do not depend on the central region's properties. The disadvantage of this approach is, however, that one needs a quite large number of buffer levels, especially at low temperatures where the Fermi function is sharp. Consequently, the many-body Hilbert space becomes very large, resulting in a challenging problem for the treatment of the correlated situation.

A very effective technique to solve low-dimensional correlated systems are matrix product states (MPS), and great progress has been made in recent years to apply MPS techniques to interacting Lindblad equations [72–85]. With this, rather large systems are within reach even though the numerical effort for Hubbard-type problems is rather high due to an extensive entanglement growth with system size [77, 85].

In this paper, rather than focusing on solution techniques for the interacting Lindblad equation, we investigate different strategies for the mapping of a general physical ('ph') reservoir onto an auxiliary ('aux') one, consisting of a small number N_B of noninteracting fermionic sites (auxiliary levels) and arbitrary Lindblad terms. This is important since the accuracy of the buffer-zone idea discussed above can be significantly improved by allowing for more general Lindblad couplings, which are determined through an optimization procedure. In this way, the associated interacting impurity solver becomes exponentially accurate with increasing N_B. Moreover, if the appropriate geometry is chosen (figure 1), already a modest value of N_B produces an accurate representation of the physical reservoir. An interesting aspect is that, as in the buffer-zone approach, the combination of intermediate bath levels and couplings to Markovian environments via Lindblad terms allows one to describe a non-Markovian bath seen on the impurity site.

Zoom In Zoom Out Reset image size

We have used these ideas in previous works [85–88] to address nonequilibrium impurity physics in the Kondo regime, with manybody solution techniques based on exact diagonalization (ED) and MPS. Especially in the latter case, we achieved very accurate results for the splitting of the Kondo peak with applied bias. Furthermore, in the equilibrium limit, we found a close agreement with numerical renormalization group (NRG) temperatures well below T_K. The sizes of the auxiliary systems were rather small, N_B = 16 for MPS and N_B = 6 for ED only,³ which demonstrates the significant improvement provided by using a Lindblad coupling in the appropriate geometry.

Here, we want to elaborate in much more detail the advantages of considering long-ranged Lindblad couplings in combination with the optimization strategy. In particular, we investigate five different geometries for the auxiliary Lindblad system, see figure 1.These setups, described in detail in section 2.6, feature different connections between the Markovian bath and the auxiliary levels. The important point is that the accuracy of our approach for fixed N_B crucially depends on the choice of the appropriate geometry. A systematic analysis of the performance of these geometries is, therefore, the main content of this paper. Besides determining the scaling of the accuracy as a function of the number of bath sites N_B, we discuss the importance of the different couplings and relate our results to the commonly used buffer-zone idea. With the applicability of manybody solution techniques in mind, we consider also the case of sparse geometries which are well-suited for MPS. Even with this limitation we find drastic differences between the different geometries. This highlights the huge potential for improvement and furthermore yields important insights into the underlying mechanisms. For ED approaches, any of the discussed geometries can be applied and one is generally interested in finding the best possible mapping for a fixed and low value of N_B. In this case the geometry with long-ranged and complex Lindblad terms outperforms the other choices, as shown below.

Finally, we also discuss results for the interacting spectral function of the Anderson impurity model in and out of equilibrium, obtained with the different geometries. In particular, we show here for the first time to our knowledge that allowing for complex Lindblad couplings produces a drastic improvement in the description of the Kondo resonance.

This paper is organized as follows: in section 2.1, we introduce the models under study and define the basic notation. In section 2.2, we illustrate the most important aspect of this work, namely the mapping of the physical Hamiltonian problem onto an auxiliary open quantum system described by a Lindblad equation. In section 2.3, we present the expressions for the non-interacting Green's function of the auxiliary system, and in section 2.4 we illustrate the fit procedure. In section 2.5, we briefly discuss the relation with the interacting case. In section 3, we present in detail the convergence of the fit as a function of N_B for the different geometries presented in section 2.6 and for different temperatures, as well as a discussion on the advantages and disadvantages of these setups. As an example, we present results for the spectral function of the Anderson impurity model. Finally, in section 4 we summarize our results and discuss possible improvements and open issues. In three appendices we present technical details of the minimization procedure (appendix A), show the explicit form of the matrices for the different geometries (appendix B), and discuss certain redundancies of the auxiliary system (appendix C).

We begin with a general discussion, which we eventually apply to the single-site Anderson impurity model. In the general case the central region may represent a small cluster or molecule. The Hamiltonian of the physical system at study is written as

$$\begin{eqnarray}&&H=\displaystyle \sum _{\alpha }({H}_{\alpha }+{H}_{\alpha f})+{H}_{f}\end{eqnarray} \tag{ 4 }$$

where ${H}_{f}$ is the Hamiltonian of the central region describing a small cluster of interacting fermions, H_α is the Hamiltonian of the reservoir $\alpha$ describing an infinite lattice of noninteracting particles, and ${H}_{\alpha f}$ is the coupling between central region and reservoirs.

$$\begin{eqnarray}&&{H}_{f}={H}_{0f}+{H}_{U}\end{eqnarray} \tag{ 5 }$$

consists of a noninteracting part

$$\begin{eqnarray}&&{H}_{0f}=\displaystyle \sum _{{ij}}{h}_{{ij}}{f}_{i}^{\dagger }{f}_{j}\end{eqnarray} \tag{ 6 }$$

and an interaction term H_U. The fermions in the reservoirs can be described by

$$\begin{eqnarray}&&{H}_{\alpha }=\displaystyle \sum _{p,p^{\prime} }{\varepsilon }_{\alpha p,p^{\prime} }{d}_{\alpha p}^{\dagger }{d}_{\alpha p^{\prime} }\end{eqnarray} \tag{ 7 }$$

in usual notation. For simplicity, spin indices are not explicitly mentioned here. Quite generally, a suitable single particle basis can be chosen such that ${\varepsilon }_{\alpha p,p^{\prime} }\sim {\delta }_{p,p^{\prime} }$. In the context of impurity problems, such a basis is often referred to as 'star representation', see also figure 1. ${H}_{\alpha f}$ is taken to be quadratic in the fermion operators:

$$\begin{eqnarray}&&{H}_{\alpha f}=\displaystyle \sum _{p,i}{v}_{\alpha {pi}}{d}_{\alpha p}^{\dagger }{f}_{i}+h.c.\,,\end{eqnarray} \tag{ 8 }$$

and d_i (f_i) are the fermionic destruction operators on the reservoir's (central region's) sites i.

We are interested in a steady-state situation, although the present approach can be easily extended to include time dependence, especially if this comes from a change of the central region's parameters only. In the steady state, the Green's functions depend only on the time difference and we can Fourier transform so that the Green's functions depend on a real frequency ω only, which here is kept implicit. We assume that initially the hybridization ${H}_{\alpha f}$ is zero and the reservoirs are separately in equilibrium with chemical potentials ${\mu }_{\alpha }$ and temperatures ${T}_{\alpha }$. Then the ${H}_{\alpha f}$ are switched on and after a certain time a steady state is reached. We use the non-equilibrium (Keldysh) formalism [89–94] whereby the Green's function can be represented as a 2 × 2 block matrix

$$\begin{eqnarray}\mathop{G}\limits_{\_}=\left(\begin{array}{cc}{G}^{{\rm{R}}} & {G}^{{\rm{K}}}\\ 0 & {G}^{{\rm{A}}}\end{array}\right)\,,\end{eqnarray} \tag{ 9 }$$

where the retarded G^R, advanced G^A, and Keldysh G^K components are matrices in the site indices (i, j) of the central region. We will adopt this underline notation in order to denote this 2 × 2 structure. We use lowercase $\underline{g}$ (${\underline{g}}_{\alpha }$) to denote Green's function of the decoupled central region (reservoir $\alpha$), while uppercase $\underline{G}$ represent the full noninteracting Green's function of the central region. For simplicity we omit the subscript ${}_{0}$, since in this paper we deal mainly with noninteracting Green's functions anyway. We use the subscript ${}_{\mathrm{int}}$ for interacting ones. ${\underline{G}}_{}$ is easily obtained from the Dyson equation as

$$\begin{eqnarray}&&{\underline{G}}_{}={({\underline{g}}^{-1}-\underline{{\rm{\Delta }}})}^{-1},\end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray}&&{\underline{{\rm{\Delta }}}}_{{ij}}=\displaystyle \sum _{\alpha ,p,p^{\prime} }{v}_{\alpha {pi}}\ {v}_{\alpha p^{\prime} j}^{* }\ {\underline{g}}_{\alpha p,p^{\prime} }\end{eqnarray} \tag{ 11 }$$

is the reservoir hybridization function (commonly called bath hybridization function) in the Keldysh representation. The retarded Green's functions ${g}_{\alpha }^{{\rm{R}}}$ for reservoirs with non-interacting fermions in equilibrium can be determined easily by standard tools, and the Keldysh components ${g}_{\alpha }^{{\rm{K}}}$ can be obtained from the retarded ones by exploiting the fluctuation-dissipation theorem:

$$\begin{eqnarray}&&{g}_{\alpha }^{{\rm{K}}}(\omega )=({g}_{\alpha }^{{\rm{R}}}(\omega )-{g}_{\alpha }^{{\rm{R}}}{(\omega )}^{\dagger })\ {s}_{\alpha }(\omega ),\end{eqnarray} \tag{ 12 }$$

which is valid since the uncoupled reservoirs are in equilibrium. Here,

$$\begin{eqnarray}&&{s}_{\alpha }(\omega )=1-2\ {f}_{{\rm{F}}}(\omega ,{\mu }_{\alpha },{T}_{\alpha })\end{eqnarray} \tag{ 13 }$$

and ${f}_{{\rm{F}}}(\omega ,{\mu }_{\alpha },{T}_{\alpha })$ is the Fermi function at chemical potential ${\mu }_{\alpha }$ and temperature ${T}_{\alpha }$.

From now on, for simplicity of presentation, we restrict to the Anderson impurity model (SIAM) in which the central region, described by equation (5), consists of a single site, i.e. there is only one value for the index i, which we drop, and

$$\begin{eqnarray}&&{H}_{U}={{Un}}_{f\uparrow }{n}_{f\downarrow }\quad \quad {n}_{f\sigma }={f}_{\sigma }^{\dagger }{f}_{\sigma }.\end{eqnarray} \tag{ 14 }$$

The idea we are going to present in section 2.2 can be immediately extended to the case of a central region consisting of many sites in which each site is connected to separate reservoirs. In the most direct fashion, this can be done with exactly the same approach as formulated here for the SIAM by just mapping each reservoir independently onto auxiliary Lindblad bath sites. An interesting application is, for example, the case of an interacting chain coupled on both sides to reservoirs with different chemical potentials [95]. Also, the extension to the case of arbitrary (quadratic) couplings with the reservoirs that intermix the central region sites relevant, e.g., for cluster dynamical mean-field theory, is conceptually straightforward although more involved.

A crucial point in the following considerations is the fact that, even in the interacting case, the influence of the reservoirs upon the central region is completely determined by the bath hybridization function $\underline{{\rm{\Delta }}}(\omega )$ only. In other words, any interacting correlation function of the central region solely depends on the central region Hamiltonian ${H}_{f}$ and on $\underline{{\rm{\Delta }}}$. This result is well known at least in equilibrium, and can be easily proven; for example, diagrammatically (see footnote⁴ ). The argument holds independently on whether one works with equilibrium or nonequilibrium Green's functions. Moreover, it crucially depends on the fact that the reservoir is noninteracting.

This can be exploited to choose different representations for the reservoir. A convenient discretization in equilibrium is to represent the reservoir by a finite number of bath sites, as commonly used in the context of NRG [25, 27] or exact-diagonalization-based dynamical-mean-field theory (ED-DMFT) [96–98]. Here, the desired physical hybridization function ${\underline{{\rm{\Delta }}}}_{\mathrm{ph}}$ is approximated by an auxiliary one ${\underline{{\rm{\Delta }}}}_{\mathrm{aux}}$,⁵ corresponding to a bath with a finite number of sites. In ED-DMFT, the parameters of these N_B bath sites are obtained through fitting the hybridization function in Matsubara space. As can be readily shown, only $2{N}_{{\rm{B}}}$ bath parameters are of relevance per spin degree of freedom, and it is common to choose a 'star' or 'chain' representation.

Clearly, the same fit strategy is inconvenient out of equilibrium for several reasons. First of all, the auxiliary system cannot dissipate since it is finite, and a steady state cannot be reached. In [86, 99] we have suggested an auxiliary master equation approach (AMEA), which adopts an auxiliary reservoir consisting of a certain number N_B of bath sites which are additionally coupled to Markovian environments described by a Lindblad equation

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}\rho ={ \mathcal L }\ \rho =({{ \mathcal L }}_{H}+{{ \mathcal L }}_{D})\rho .\end{eqnarray} \tag{ 15 }$$

Here, the Hamiltonian for the auxiliary system is given by (we reintroduce spin)

$$\begin{eqnarray}&&{H}_{\mathrm{aux}}=\displaystyle \sum _{{ij}\sigma }{E}_{{ij}}{c}_{i\sigma }^{\dagger }{c}_{j\sigma }+{{Un}}_{f\uparrow }{n}_{f\downarrow },\end{eqnarray} \tag{ 16 }$$

and enters the unitary part of the Lindblad operator

$$\begin{eqnarray}&&{{ \mathcal L }}_{H}\rho =-i[{H}_{\mathrm{aux}},\rho ].\end{eqnarray} \tag{ 17 }$$

The dissipator ${{ \mathcal L }}_{D}$ describes the coupling of the auxiliary sites to the Markovian environment and is given by

$$\begin{eqnarray}&&{{ \mathcal L }}_{D}\rho =2\displaystyle \sum _{{ij}\sigma }{{\rm{\Gamma }}}_{{ij}}^{(1)}\left({c}_{j\sigma }\rho {c}_{i\sigma }^{\dagger }-\displaystyle \frac{1}{2}\{\rho ,{c}_{i\sigma }^{\dagger }{c}_{j\sigma }\}\right)+2\displaystyle \sum _{{ij}\sigma }{{\rm{\Gamma }}}_{{ij}}^{(2)}\left({c}_{i\sigma }^{\dagger }\rho {c}_{j\sigma }-\displaystyle \frac{1}{2}\{\rho ,{c}_{j\sigma }{c}_{i\sigma }^{\dagger }\}\right).\end{eqnarray} \tag{ 18 }$$

in terms of matrices ${{\boldsymbol{\Gamma }}}^{(1/2)}$ with matrix elements ${{\rm{\Gamma }}}_{{ij}}^{(1/2)}$ to be determined by an optimization procedure, as discussed below. The indices $i,j$ in equations (17) and (18) run over the impurity i = f (we identify ${c}_{f\sigma }={f}_{\sigma }$) and over the N_B bath sites⁶ . Similar to the case of the ED-DMFT impurity solver previously mentioned, the idea is to optimize the parameters of the auxiliary reservoir in order to achieve a best fit to the physical bath hybridization function equation (11), i.e., for a given N_B, ${\underline{{\rm{\Delta }}}}_{\mathrm{aux}}(\omega )$ should be as close as possible to ${\underline{{\rm{\Delta }}}}_{\mathrm{ph}}(\omega )$:

$$\begin{eqnarray}&&{\underline{{\rm{\Delta }}}}_{\mathrm{aux}}(\omega )\approx {\underline{{\rm{\Delta }}}}_{\mathrm{ph}}(\omega ).\end{eqnarray} \tag{ 19 }$$

As for the ED-DMFT case, one can choose a single-particle basis for the auxiliary bath such that the matrix ${\boldsymbol{E}}$ is sparse,⁷ , i.e. it has a 'star' or a 'chain' form and is real valued. However, there is no reason why the matrices ${{\boldsymbol{\Gamma }}}^{(1)}$ and ${{\boldsymbol{\Gamma }}}^{(2)}$ should be sparse and real in the same basis as well, and, in fact, as discussed below, for an ED treatment of the Lindblad problem it is convenient to allow for a general form in order to optimize the fit. This larger number of parameters allows one to fulfill equation (19) to a very good approximation. The introduction of dissipators equation (18) additionally allows us to carry out the fit directly for real ω (see section 2.4) since ${\underline{{\rm{\Delta }}}}_{\mathrm{aux}}(\omega )$ is a continuous function. This makes this approach competitive with ED-DMFT for the equilibrium case as well.

Notice that equation (18) is not the most general form of the dissipator, and one could think of including Lindblad terms that contain four or more fermionic operators or also anomalous and spin-flip terms. This would increase the number of parameters available for the fit. However, the latter would violate conserved quantities and the former would describe an interacting bath, so that the argument of section 2.2 (see footnote 3) does not apply. As a matter of fact, the exact equivalence to a noninteracting bath (see footnote 6) only holds for a quadratic form of the Lindblad operator as in equation (18).

Once the optimal values of the matrices ${\boldsymbol{E}}$, ${{\boldsymbol{\Gamma }}}^{(1)}$ and ${{\boldsymbol{\Gamma }}}^{(2)}$ for a given physical model are determined for the non-interacting system, one could solve for the dynamics of the correlated auxiliary system defined by equation (15), which amounts to a linear equation for the reduced many-body density matrix. If the number of sites of this system is small, one can solve exactly for the steady state and the dynamics of this interacting system by methods such as Lanczos exact diagonalization or matrix product states (MPS) [81–83, 85].

In order to carry out the fit equation (19), we need to compute the auxiliary reservoir hybridization function ${\underline{{\rm{\Delta }}}}_{\mathrm{aux}}(\omega )$ for many values of the bath and Lindblad parameters. This can be done in an efficient manner since only noninteracting Green's functions are needed, see also equation (10) and the preceding discussion. Computing the single-particle Green's function matrix $\underline{{\bf{Z}}}$, obtained from the steady-state solution of equation (15), amounts to solving a noninteracting fermion problem, which scales polynomially with respect to the single-particle Hilbert space . A method to deal with quadratic fermions with linear dissipation based on a so-called 'third quantization' has been introduced in [100]. We adopt the approach of [69] in which the authors recast an open quantum problem like equation (15) into a standard operator problem in an augmented fermion Fock space with twice as many sites and with a non-Hermitian Hamiltonian [69, 101, 102]. This so-called super-fermionic representation is convenient for our purposes, not only to solve for the noninteracting Green's functions but also to treat the many-body problem in an analogous framework to Hermitian problems. An analytic expression for the noninteracting steady-state retarded and Keldysh auxiliary Green's functions was derived in [86]. An alternative derivation, which does not rely on super-fermions, is given in [71]. For the retarded component, we get (see footnote 7)

$$\begin{eqnarray}&&{{\bf{Z}}}_{}^{{\rm{R}}}(\omega )={(\omega -{\boldsymbol{E}}+{i}({{\boldsymbol{\Gamma }}}^{(2)}+{{\boldsymbol{\Gamma }}}^{(1)}))}^{-1},\end{eqnarray} \tag{ 20 }$$

and the Keldysh component of the inverse Green's function reads

$$\begin{eqnarray}&&{({\mathop{{\bf{Z}}}\limits_{\_}}^{-1})}^{{\rm{K}}}=-2{i}({{\boldsymbol{\Gamma }}}^{(2)}-{{\boldsymbol{\Gamma }}}^{(1)}),\end{eqnarray} \tag{ 21 }$$

yielding the Keldysh Green's function

$$\begin{eqnarray}{{\bf{Z}}}_{}^{{\rm{K}}} & = & -\,{{\bf{Z}}}_{}^{{\rm{R}}}{({\mathop{{\bf{Z}}}\limits_{\_}}_{}^{-1})}^{{\rm{K}}}{{\bf{Z}}}_{}^{{\rm{A}}}\\ & = & 2{i}{{\bf{Z}}}_{}^{{\rm{R}}}({{\boldsymbol{\Gamma }}}^{(2)}-{{\boldsymbol{\Gamma }}}^{(1)}){{\bf{Z}}}_{}^{{\rm{A}}}.\end{eqnarray} \tag{ 22 }$$

The ff component of $\underline{{\bf{Z}}}$ is the auxiliary impurity Green's function

$$\begin{eqnarray}&&{\underline{G}}_{\mathrm{aux}}={(\underline{{\bf{Z}}})}_{{ff}}.\end{eqnarray} \tag{ 23 }$$

From this, one can determine the retarded component of ${\underline{{\rm{\Delta }}}}_{\mathrm{aux}}(\omega )$

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\mathrm{aux}}^{{\rm{R}}}(\omega )=1/{g}^{{\rm{R}}}-1/{G}_{\mathrm{aux}}^{{\rm{R}}}.\end{eqnarray} \tag{ 24 }$$

For the Keldysh component, one has to carry out two inversions of Keldysh matrices (see, e.g. [92]), yielding

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\mathrm{aux}}^{{\rm{K}}}(\omega )=-{({\mathop{G}\limits_{\_}}_{\mathrm{aux}}^{-1})}^{{\rm{K}}}=1/| {G}_{\mathrm{aux}}^{{\rm{R}}}{| }^{2}\ {G}_{\mathrm{aux}}^{{\rm{K}}},\end{eqnarray} \tag{ 25 }$$

where the contribution from g^K is infinitesimal.

From the preceding equations, we can efficiently compute ${\underline{{\rm{\Delta }}}}_{\mathrm{aux}}(\omega )$ for a given set of parameters of the auxiliary reservoir. The numerical effort for a single evaluation is low and scales only at most as ${ \mathcal O }({N}_{{\rm{B}}}^{3})$. We introduce a vector of parameters ${\boldsymbol{x}}$ which yields a unique set of matrices ${\boldsymbol{E}}$, ${{\boldsymbol{\Gamma }}}^{(1)}$ and ${{\boldsymbol{\Gamma }}}^{(2)}$ within a chosen subset (see, e.g. figure 1 and appendix B) and quantify the deviation from equation (19) through a cost function

$$\begin{eqnarray}\chi {({\boldsymbol{x}})}^{2} & = & \displaystyle \frac{1}{{\chi }_{0}^{2}}{\displaystyle \int }_{-{\omega }_{{\rm{c}}}}^{{\omega }_{{\rm{c}}}}| {\mathop{{\rm{\Delta }}}\limits_{\_}}_{\mathrm{ph}}-{\mathop{{\rm{\Delta }}}\limits_{\_}}_{\mathrm{aux}}| W(\omega ){\rm{d}}\omega ,\\ | {\mathop{{\rm{\Delta }}}\limits_{\_}}_{\mathrm{ph}}-{\mathop{{\rm{\Delta }}}\limits_{\_}}_{\mathrm{aux}}| & = & \displaystyle \sum _{\xi \in \{{\rm{R}},{\rm{K}}\}}{\mathfrak{I}}m\,\{{{\rm{\Delta }}}_{\mathrm{ph}}^{\xi }(\omega )-{{\rm{\Delta }}}_{\mathrm{aux}}^{\xi }(\omega ;{\boldsymbol{x}})\}{}^{2},\end{eqnarray} \tag{ 26 }$$

and minimize $\chi ({\boldsymbol{x}})$ with respect to ${\boldsymbol{x}}$. The normalization ${\chi }_{0}$ is hereby chosen such that $\chi ({\boldsymbol{x}})=1$ when ${\underline{{\rm{\Delta }}}}_{\mathrm{aux}}(\omega )\equiv 0$. It is important to note that both the retarded and the Keldysh component must be fitted. Due to Kramers–Kronig relations, the real part of ${{\rm{\Delta }}}_{\mathrm{ph}}^{{\rm{R}}}(\omega )$ is fully determined by its imaginary part, provided the asymptotic behavior is fixed. Therefore, we can restrict to fit its imaginary part, while ${{\rm{\Delta }}}_{\mathrm{ph}}^{{\rm{K}}}(\omega )$ is purely imaginary. Furthermore, in equation (26) we introduced a cut-off frequency ${\omega }_{{\rm{c}}}$ and a weighting function $W(\omega )$. In this work, we take $W(\omega )=1$ and ${\omega }_{{\rm{c}}}=1.5\,D$, with D the half-bandwidth of ${\underline{{\rm{\Delta }}}}_{\mathrm{ph}}(\omega )$. Different forms of $W(\omega )$ can be used, for instance in order to increase the accuracy of the fit near the chemical potentials. In general, other definitions for equation (26) are possible such as a piecewise definition with varying interval lengths. By this, one may combine the present approach with NRG ideas such as the logarithmic discretization, and work along these lines is in progress (see also [71]). On the whole, the minimization of equation (26) constitutes a multi-dimensional optimization problem and appropriate numerical methods for it are discussed in appendix A.

As asymptotic limit, we require here ${\underline{{\rm{\Delta }}}}_{\mathrm{aux}}(\omega )\to 0$ for $\omega \to \pm \infty$, which is obtained when setting ${{\rm{\Gamma }}}_{{ff}}^{(1/2)}=0$. Semipositivity further requires ${{\rm{\Gamma }}}_{{if}}^{(1/2)}={{\rm{\Gamma }}}_{{fi}}^{(1/2)}=0$. For simplicity, we restrict here to the particle-hole symmetric case. This reduces the number of free parameters in ${\boldsymbol{E}}$, ${{\boldsymbol{\Gamma }}}^{(1)}$ and ${{\boldsymbol{\Gamma }}}^{(2)}$. For the case in which the impurity site f is located in the center and that one has an even number of bath sites N_B, particle-hole symmetry in the auxiliary system is obtained when

$$\begin{eqnarray}{E}_{{ij}} & = & {(-1)}^{i+j+1}{E}_{{N}_{{\rm{B}}}+2-j,{N}_{{\rm{B}}}+2-i},\\ {{\rm{\Gamma }}}_{{ij}}^{(1)} & = & {(-1)}^{i+j}{{\rm{\Gamma }}}_{{N}_{{\rm{B}}}+2-j,{N}_{{\rm{B}}}+2-i}^{(2)},\end{eqnarray} \tag{ 27 }$$

for $i,j\in \{1,\ldots ,\,{N}_{B}+1\}$. More details for the particular form of ${\boldsymbol{E}}$, ${{\boldsymbol{\Gamma }}}^{(1)}$, and ${{\boldsymbol{\Gamma }}}^{(2)}$ are given in appendix B.

We now briefly discuss here some relevant issues in connection with the evaluation of particular observables of the physical system from results of the auxiliary system. More details can be found in [85, 86].

As already discussed, by mapping onto an auxiliary interacting open quantum system of finite size described by the Lindblad equation (equation (15)), we obtain a manybody problem which can be solved exactly or at least with high numerical precision, provided N_B is not too large. In [86], we presented a solution strategy based on exact diagonalization (ED) with Krylov space methods, and in [85] one based on MPS. In the end, both techniques allow us to determine the interacting impurity Green's function ${\underline{G}}_{\mathrm{aux},\mathrm{int}}(\omega )$ of the auxiliary system. As discussed previously, in the limit ${\underline{{\rm{\Delta }}}}_{\mathrm{aux}}(\omega )\to {\underline{{\rm{\Delta }}}}_{\mathrm{ph}}(\omega )$ (i.e. for large N_B) this becomes equivalent to the physical one ${\underline{G}}_{\mathrm{ph},\mathrm{int}}(\omega )$. However, this equivalence only holds for impurity correlation functions, and, for example, it does not apply for the current flowing from a left ($\alpha =l$) to a right ($\alpha =r$) reservoir across the impurity. Therefore, the current evaluated within the auxiliary Lindblad system does not necessarily correspond to the physical current even for large N_B,⁸ unless one fits the bath hybridization functions ${\underline{{\rm{\Delta }}}}_{\mathrm{ph},\alpha }(\omega )$ for the left and right reservoirs separately. Such a separate fit, however, is not necessary and would simply worsen the overall accuracy for a given N_B. Once the approximate ${\underline{G}}_{\mathrm{ph},\mathrm{int}}(\omega )\approx {\underline{G}}_{\mathrm{aux},\mathrm{int}}(\omega )$ is known, the current of the physical system can be evaluated by means of the well-known Meir–Wingreen expression [92, 103, 104], albeit by using the Fermi functions and density of states (hybridization functions) of the two physical reservoirs separately. Therefore, the knowledge of ${\underline{G}}_{\mathrm{aux},\mathrm{int}}(\omega )$ enables one to compute most quantities of interest.

An additional step consists in extracting just the self-energy from the solution of the auxiliary impurity system

$$\begin{eqnarray*}&&{\underline{{\rm{\Sigma }}}}_{\mathrm{aux}}(\omega )={\underline{G}}_{\mathrm{aux}}^{-1}(\omega )-{\underline{G}}_{\mathrm{aux},\mathrm{int}}^{-1}(\omega ).\end{eqnarray*}$$

and inserting it into the Dyson equation for the physical system with the exact physical noninteracting Green's function

$$\begin{eqnarray}&&{\underline{G}}_{\mathrm{ph},\mathrm{int}}(\omega )\approx {({\underline{G}}_{\mathrm{ph}}{(\omega )}^{-1}-{\underline{{\rm{\Sigma }}}}_{\mathrm{aux}}(\omega ))}^{-1}.\end{eqnarray} \tag{ 28 }$$

Clearly, this step is only useful when the relation equation (19) is approximate, since for ${\underline{{\rm{\Delta }}}}_{\mathrm{aux}}(\omega )\to {\underline{{\rm{\Delta }}}}_{\mathrm{ph}}(\omega )$ also the noninteracting Green's functions ${\underline{G}}_{\mathrm{ph}}(\omega )$ and ${\underline{G}}_{\mathrm{aux}}(\omega )$ would coincide, i.e. in the hypothetical ${N}_{{\rm{B}}}\to \infty$ case, and one could just set ${\underline{G}}_{\mathrm{ph},\mathrm{int}}(\omega )\to {\underline{G}}_{\mathrm{aux},\mathrm{int}}(\omega )$. For finite N_B, this substitution has the advantage that in equation (28) the noninteracting part ${\underline{G}}_{\mathrm{ph}}(\omega )$ is exact, and the approximation equation (19) only affects the self energy.

With the goal in mind of providing the best approximation to the full interacting impurity problem described by the Hamiltonian equation (4), we would like to approximate ${\underline{{\rm{\Delta }}}}_{\mathrm{ph}}(\omega )$ by ${\underline{{\rm{\Delta }}}}_{\mathrm{aux}}(\omega )$ as accurately as possible for a given number of bath sites N_B. In principle, one has the freedom to choose different geometries for the auxiliary system, and a generic set of five different setups is depicted in figure 1. (An explicit form of the corresponding matrices for N_B = 4 is given in appendix B.) For large N_B all geometries converge to the exact solution ${\underline{{\rm{\Delta }}}}_{\mathrm{aux}}(\omega )\to {\underline{{\rm{\Delta }}}}_{\mathrm{ph}}(\omega )$, and the question is how fast. In section 3, we want to elaborate on this point in detail and present results obtained with those geometries, which we briefly discuss and motivate here.

In all cases one can restrict the geometries to a sparse (e.g. tridiagonal) and real-valued matrix ${\boldsymbol{E}}$. As is commonly true for impurity problems, the physics on the impurity site is invariant under unitary transformations among bath sites only. For an arbitrary unitary transformation ${\boldsymbol{U}}$ with ${U}_{{if}}={U}_{{fi}}={\delta }_{{if}}$ to new fermionic operators, one obtains an analogous auxiliary system with modified bath parameters ${\boldsymbol{E}}^{\prime} ={{\boldsymbol{U}}}^{\dagger }{\boldsymbol{E}}{\boldsymbol{U}}$, ${{\boldsymbol{\Gamma }}}^{{({\bf{1}})}^{\prime }}={{\boldsymbol{U}}}^{\dagger }{{\boldsymbol{\Gamma }}}^{({\bf{1}})}{\boldsymbol{U}}$ and ${{\boldsymbol{\Gamma }}}^{{({\bf{2}})}^{\prime }}={{\boldsymbol{U}}}^{\dagger }{{\boldsymbol{\Gamma }}}^{({\bf{2}})}{\boldsymbol{U}}$. It is easy to check that the ff-component of the Green's functions equations (20) and (22) is not affected by this transformation. Therefore, we choose without loss of generality ${\boldsymbol{E}}$ to be sparse as well as real, and for ${{\boldsymbol{\Gamma }}}^{(1/2)}$ in the most general case dense matrices with ${ \mathcal O }({N}_{{\rm{B}}}^{2})$ parameters. The particular form of ${\boldsymbol{E}}$ is irrelevant, i.e. whether it is diagonal for bath sites (star) or tridiagonal (chain), as long as the ${{\boldsymbol{\Gamma }}}^{(1/2)}$ matrices are transformed accordingly.

Such a general geometry with sparse ${\boldsymbol{E}}$ and dense ${{\boldsymbol{\Gamma }}}^{(1/2)}$ is referred to as 'full' setup in the following. Here, we will further distinguish between the case in which the ${{\boldsymbol{\Gamma }}}^{(1/2)}$ are real or they have complex elements ('full complex'). In addition, we consider the four sparse cases '2 chains n.n.', '2 chains onsite', 'star', and '1 chain n.n.', in which also the ${{\boldsymbol{\Gamma }}}^{(1/2)}$ are sparse. The meaning of these abbreviations is given in figure 1, see also appendix B. These sparse geometries are, however, not linked to each other by unitary transformations and represent nonequivalent subsets of the 'full' setup. Which one of these is advantageous in practice is not obvious a priori, and is discussed in the next section⁹ .

The 'full' geometry comprises all other ones and thus, obviously, gives the best possible fit for a given N_B. In addition, one can allow for the off-diagonal matrix elements of the ${{\boldsymbol{\Gamma }}}^{(1/2)}$ to be complex, thus extending the set of fit parameters. Nevertheless, the sparse setups may be of great value for sophisticated manybody solution strategies for the interacting Lindblad equation, such as MPS. We made use of the 'full' setup (with real parameters) in the ED treatment [86], which is applicable to dense ${{\boldsymbol{\Gamma }}}^{({\bf{1}})}$ and ${{\boldsymbol{\Gamma }}}^{({\bf{2}})}$ matrices and could consider up to N_B = 6. Larger systems are prohibitive due to the exponentially increasing Hilbert space (see footnote [2]). On the other hand, within MPS, we can currently consider up to N_B = 20 bath sites. However, in favor of the applicability of MPS methods, one should avoid long-ranged hoppings and we thus employed the '2 chains n.n.' geometry. As becomes evident also from the results below, the gain in N_B hereby outweighs the restriction of the fit setup, so that the MPS approach is clearly more accurate. Also, the other sparse setups investigated below are possible candidates for MPS, see also [105]. Besides this, approaches such as the previously mentioned buffer zone scheme and variations of it [68–70], which are often applied concepts in Lindblad-type representation of noninteracting environments, are related to the 'star' geometry; see also our later discussion.

We now present and discuss results for the interacting spectral function of the SIAM,

$$\begin{eqnarray}&&A(\omega )=-\displaystyle \frac{1}{\pi }{\mathfrak{I}}m\,{G}_{\mathrm{ph},\mathrm{int}}^{R}(\omega )\end{eqnarray} \tag{ 29 }$$

(see equation (28)) in and out of equilibrium and in particular investigate the accuracy obtained by each of the geometries of figure 1 for a fixed N_B. We also discuss the relation with the fidelity in reproducing the physical hybridization function, i.e. the fit just discussed. In figures 8 and 9, we confront the results of the fit on the left with the interacting spectral function on the right. Results for N_B = 6 are presented for two values of the bias voltage, $\phi =0$ (equilibrium case) and $\phi =3{\rm{\Gamma }}$. In particular, we focus on the region around the Kondo peak, which gets split at finite ϕ, since this is the more sensitive to approximations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Starting with the equilibrium case, we see that all but the two worse geometries are able to resolve the Hubbard side bands. Further improvements of the fit show that, the better an auxiliary system reproduces the Fermi jump in the Keldysh and the plateau in the retarded hybridization function, the sharper is the corresponding Kondo resonance. Only the '2 chains n.n.' setup seems to break this trend as, for $\phi =0$, its peak is clearly sharper than for the 'full' geometry although the fit is less accurate. However, this peculiarity can be readily explained by looking at the corresponding fit and in particular at the oscillations in the retarded component at low energies. Here, the 'full' geometry shows a substantial dip which in turn leads to a suppression of the Kondo peak. A way to improve this is to introduce a weight $W(\omega )$ (cf equation (26)) to the fit that emphasizes regions around the chemical potentials. One can also notice the drastic improvement in the description of the Kondo resonance obtained by allowing for complex values of the parameters ${{\rm{\Gamma }}}_{{ij}}^{(1/2)}$. The same accuracy is expected to be obtained for the real 'full' case for N_B = 8, but this is currently beyond the reach of a Krylov-space solution method.

Turning to the nonequilibrium case, we first note that the 'star' and '1 chain' setup are not able to follow the double Fermi step in the Keldysh hybridization function but rather produce a single jump at an elevated temperature. Thus, assuming for the moment a flat retarded component, one would expect only a single temperature-broadened Kondo resonance. The fact that the spectral function instead shows a two-peak structure is not a genuine effect but rather is connected to corresponding oscillations in the retarded component of the hybridization function. Notice that the two-peak structures occur here at the wrong frequency. In fact, all but the 'full complex' and 'full' setups fail to correctly reproduce the positions of the peaks, despite resolving the double Fermi step, as oscillations in the retarded component interfere with the physical effect. In general, this example shows that one has to be careful with interpretations of structures in the interacting spectral function when the fit does not resolve important features and/or displays unphysical oscillations. On the other hand, if deviations from the physical hybridization function are small, one can be confident that the results are correct. Results suggest this to be the case already with 6 bath sites for the 'full' setups and about 12 bath sites for the '2 chains' geometries.

The PT algorithm is outlined in detail in the following. For completeness, let us first briefly recap the basic ideas of the underlying Markov chain Monte Carlo (MCMC), and of the related simulated annealing algorithm.

MCMC techniques were originally developed to evaluate thermodynamic properties of classical systems which exhibit a very large phase space where simple sampling strategies fail. For our purposes here, we are interested in minimizing the cost function $\chi ({\boldsymbol{x}})$ as defined in equation (26) with respect to the parameter vector ${\boldsymbol{x}}$. For such high-dimensional minimization problems, one can adapt MCMC schemes by viewing $\chi ({\boldsymbol{x}})$ as an artificial energy and by introducing an artificial inverse temperature β. In the so-called simulated annealing, one samples from the Boltzmann distribution $p({\boldsymbol{x}})=1/Z\exp (-\chi ({\boldsymbol{x}})\beta )$ at a certain β, and then successively cools down the artificial temperature. Motivated by the behavior of true physical systems, one expects to end up in the low-energy state when letting the system equilibrate and when cooling sufficiently slowly. Analogous to thermodynamics, one can calculate the specific heat ${C}_{{\rm{H}}}={\beta }^{2}\langle {\rm{\Delta }}\chi {({\boldsymbol{x}})}^{2}\rangle$ and by this locate regions with large changes, i.e. phase transitions, where a slow cooling is critical. However, in practice it may be time consuming to realize the equilibration and sufficiently slow cooling, and for tests within AMEA we often ended up in local minima. In order to obtain a robust algorithm which can also start from previous solutions as needed, for instance, within DMFT, we sought a method which is able to efficiently overcome local minima and still systematically target the low-energy states. For this, a multicanonical and PT algorithm were tested and the latter turned out to be more convenient. In the following, we briefly outline the PT scheme used within AMEA and refer to [106–110] for a thorough introduction to MCMC, simulated annealing, multicanonical sampling, and PT.

As just stated, in a MCMC scheme, one typically samples from the Boltzmann distribution $p({\boldsymbol{x}})=1/Z\exp (-\chi ({\boldsymbol{x}})\beta )$ at some chosen inverse temperature β. This is done through an iteratively created chain of states $\{{{\boldsymbol{x}}}_{l}\}$, whereby one avoids the explicit calculation of the partition function Z. An effective and well-known scheme for this is the Metropolis–Hastings algorithm [106, 107]. One starts out with some state ${{\boldsymbol{x}}}_{l}$ and proposes a new configuration ${{\boldsymbol{x}}}_{k}$, whereby it has to be ensured that every state of the system can be reached in order to achieve ergodicity. The proposed state ${{\boldsymbol{x}}}_{k}$ is accepted with probability [106, 107]¹⁴

$$\begin{eqnarray}&&{p}_{\mathrm{pacc}.}^{l,k}=\min \left\{1,\displaystyle \frac{p({{\boldsymbol{x}}}_{k})}{p({{\boldsymbol{x}}}_{l})}\right\}=\min \{1,{{\rm{e}}}^{-(\chi ({{\boldsymbol{x}}}_{k})-\chi ({{\boldsymbol{x}}}_{l}))\beta }\}.\end{eqnarray} \tag{ A.1 }$$

If the proposed configuration is accepted, then the next element ${{\boldsymbol{x}}}_{l+1}$ in the chain is ${{\boldsymbol{x}}}_{k}$, otherwise ${{\boldsymbol{x}}}_{l}$ again. From equation (A.1) it is obvious that ${p}_{\mathrm{pacc}.}^{l,k}=1$ when $p({{\boldsymbol{x}}}_{k})\gt p({{\boldsymbol{x}}}_{l})$, so that an importance of sampling towards regions where $p({\boldsymbol{x}})$ is large is achieved. One can show that the algorithm fulfills detailed balance and draws a set of samples $\{{{\boldsymbol{x}}}_{l}\}$ that follow the desired distribution $p({\boldsymbol{x}})$. However, stemming from the iterative construction, correlations in the chain are present which require a careful analysis for the purpose of statistical physics [106, 107]. For optimization problems, on the other hand, the situation is much simpler and one is just interested in the element in $\{{{\boldsymbol{x}}}_{l}\}$ which minimizes $\chi ({\boldsymbol{x}})$. Since a proposed step with $\chi ({{\boldsymbol{x}}}_{k})\lt \chi ({{\boldsymbol{x}}}_{l})$ is always accepted, the algorithm targets minima; however, uphill moves in configuration space are also allowed with a probability depending exponentially on the barrier height ${\rm{\Delta }}{\chi }_{k,l}=\chi ({{\boldsymbol{x}}}_{k})-\chi ({{\boldsymbol{x}}}_{l})$ and β. Effectively, uphill moves only take place when ${\rm{\Delta }}{\chi }_{k,l}\beta \lesssim { \mathcal O }(1)$. For small values of β, large moves in configuration space with large ${\rm{\Delta }}{\chi }_{k,l}$ are likely to be accepted, whereas for large β the distribution $p({\boldsymbol{x}})$ is peaked at minima in $\chi ({\boldsymbol{x}})$ so that those regions are especially sampled. For the latter case, configurations in the chain $\{{{\boldsymbol{x}}}_{l}\}$ are generally more correlated and once a ${{\boldsymbol{x}}}_{l}$ corresponds to a local minimum, the algorithm may stay there for a very long time.

One has great freedom in defining a proposal distribution from which the new state ${{\boldsymbol{x}}}_{k}$ is drawn given the current configuration ${{\boldsymbol{x}}}_{l}$ ¹⁵ . Common choices are, for instance, a Gaussian or a Lorentzian distribution with the vector difference ${{\boldsymbol{x}}}_{k}-{{\boldsymbol{x}}}_{l}$ as argument. We favored the former and updated each component i with a probability according to [106]

$$\begin{eqnarray}&&{q}_{l,k}^{i}=\displaystyle \frac{1}{\sqrt{2\pi }{\sigma }_{i}}{e}^{-\tfrac{{({{\boldsymbol{x}}}_{k}-{{\boldsymbol{x}}}_{l})}_{i}^{2}}{2{\sigma }_{i}^{2}}}.\end{eqnarray} \tag{ A.2 }$$

Hereby, a different step size ${\sigma }_{i}$ for each component is expedient since the potential landscape $\chi ({\boldsymbol{x}})$ around ${{\boldsymbol{x}}}_{l}$ is typically highly anisotropic. Ideally, one should make use of the covariance matrix ${{\rm{\Sigma }}}_{l}$ of $\chi ({{\boldsymbol{x}}}_{l})$ and consider as argument for the Gaussian instead ${({{\boldsymbol{x}}}_{k}-{{\boldsymbol{x}}}_{l})}^{T}{{\rm{\Sigma }}}_{l}^{-1}({{\boldsymbol{x}}}_{k}-{{\boldsymbol{x}}}_{l})$ [106]. However, we encountered a problem in that the estimation of the covariance matrix at run time was strongly affected by noise and thus not feasible. The adjustment of the step sizes ${\sigma }_{i}$, on the contrary, can be done after a short number of updates by demanding that a value of ${p}_{\mathrm{pacc}.}^{l,k}\approx 0.5$ is reached on average when modifying the component i. For this, we implemented a check at every single proposal that increases ${\sigma }_{i}\to 1.1{\sigma }_{i}$ when ${p}_{\mathrm{pacc}.}^{l,k}\gt 0.6$ and decreases ${\sigma }_{i}\to 0.9\,{\sigma }_{i}$ when ${p}_{\mathrm{pacc}.}^{l,k}\lt 0.4$. Analogous to the treatment of spin systems, we define one sweep as a single update of all the components of ${\boldsymbol{x}}$ ¹⁶ .

In a PT algorithm one considers, instead of sampling at one certain temperature, a set of different temperatures ${\beta }_{m}^{-1}$ and corresponding replicas ${{\boldsymbol{x}}}_{l}^{m}$, each of which is evolved through a Markov chain. The largest ${\beta }_{m}$ thereby target local minima whereas low ${\beta }_{m}$ values allow for large moves in configuration space. The key idea of the PT approach is to let the individual replicas evolve dynamically in the set of ${\beta }_{m}$. By this, one achieves a situation whereby a replica at high ${\beta }_{m}$ values systematically targets local minima but can overcome potential barriers again when its inverse temperature is changed to lower values. As a result, the time scales to reach an absolute minimum are drastically reduced and an efficient sampling of the low-energy states is achieved. For the purpose of calculating thermodynamic properties, one usually chooses a Metropolis–Hastings probability to swap two replicas with adjacent temperatures [109, 110]

$$\begin{eqnarray}&&{p}_{\mathrm{swap},l}^{m,m+1}=\min \left\{1,\displaystyle \frac{{p}^{m}({{\boldsymbol{x}}}_{l}^{m+1}){p}^{m+1}({{\boldsymbol{x}}}_{l}^{m})}{{p}^{m+1}({{\boldsymbol{x}}}_{l}^{m+1}){p}^{m}({{\boldsymbol{x}}}_{l}^{m})}\right\}=\min \{1,{{\rm{e}}}^{({\beta }_{m+1}-{\beta }_{m})(\chi ({{\boldsymbol{x}}}_{l}^{m+1})-\chi ({{\boldsymbol{x}}}_{l}^{m}))}\},\end{eqnarray} \tag{ A.3 }$$

with the Boltzmann distribution for each ${\beta }_{m}$ given by ${p}^{m}({\boldsymbol{x}})=1/{Z}_{m}\exp (-\chi ({\boldsymbol{x}}){\beta }_{m})$. Such swap moves are conveniently proposed after a certain number of sweeps, which satisfies the sufficient condition of balance for thermodynamics [110]. In practice, we chose 10 sweeps before swapping replicas. For the exchange to effectively take place, the underlying requirement is that the adjacent ${\beta }_{m}$ and ${\beta }_{m+1}$ values be close enough to each other, so that the two energy distributions ${\rm{\Omega }}[\chi ({\boldsymbol{x}})]{p}^{m}({\boldsymbol{x}})$ and ${\rm{\Omega }}[\chi ({\boldsymbol{x}})]{p}^{m+1}({\boldsymbol{x}})$ overlap, with ${\rm{\Omega }}[{\chi }_{0}]=\int {\rm{d}}{\boldsymbol{x}}\ \delta ({\chi }_{0}-\chi ({\boldsymbol{x}}))$ being the density of states of the cost function. This means that a replica at one temperature must represent a likely configuration for the neighboring temperature [110, 111]. In order to achieve this, a crucial point in the PT algorithm is to adjust the distribution of the inverse temperatures properly to the considered situation. Various criteria for this have been devised, see e.g. [110]. A common choice is to demand that the swapping probability equation (A.3) become constant as a function of temperature [112, 113], and in [114] a feedback strategy was presented which optimizes the round trip times of replicas. We tested the latter within AMEA but favored the simpler former criterion in the end, since it allows for a rapid feedback and quick adjustment to large changes in $\chi ({{\boldsymbol{x}}}_{l}^{m})$. In the simple situation of a constant specific heat C_H with respect to energy χ for instance, an optimal strategy is known since a geometric progression ${\beta }_{m}/{\beta }_{m+1}=\mathrm{const}.$ of temperatures yields a constant swapping probability [110, 111]. For interesting cases, in practice, this is rarely fulfilled, but within AMEA it served as a good starting point. The set of inverse temperatures is then optimized by averaging ${p}_{\mathrm{swap},l}^{m,m+1}$ over a couple of swappings to obtain the mean probability ${\bar{p}}_{\mathrm{swap}}^{m,m+1}$ and adjusting the ${\beta }_{m}$ thereafter. For this we chose a fixed lowest and highest ${\beta }_{m}$ value and changed the spacings in between according to

$$\begin{eqnarray}&&{\rm{\Delta }}{\beta }_{m}^{\prime }=c\displaystyle \frac{{\rm{\Delta }}{\beta }_{m}}{\mathrm{log}({\bar{p}}_{\mathrm{swap}}^{m,m+1})},\end{eqnarray} \tag{ A.4 }$$

with ${\rm{\Delta }}{\beta }_{m}={\beta }_{m+1}-{\beta }_{m}$ and c adjusted properly so that $\max ({\beta }_{m}^{\prime })-\min ({\beta }_{m}^{\prime })=\max ({\beta }_{m})-\min ({\beta }_{m})$. In [112, 113], it was shown that a constant swapping probability of 20%–23% seems to be optimal. We determined the highest and lowest ${\beta }_{m}$ values by the changes in $\chi ({\boldsymbol{x}})$ we want to resolve or allow for, and the number of inverse temperatures ${\beta }_{m}$ was then set accordingly in order to roughly obtain ${\bar{p}}_{\mathrm{swap}}^{m,m+1}\approx 0.25$. Fixing the smallest and largest ${\beta }_{m}$ is, for our purposes, the most convenient choice among the many possibilities.

However, despite the feedback optimization of temperatures as just described, we encountered unwanted behavior in practice, whereby the set of parallel replicas effectively decoupled into several clusters. In suppressing this, we found it advantageous to introduce the following simple modification to equation (A.3)¹⁷ :

$$\begin{eqnarray}&&{p}_{\mathrm{swap},l}^{m,m+1}=\max \{{p}_{\mathrm{swap},l}^{m,m+1},{p}_{\mathrm{swap}}^{\mathrm{th}.}\},\end{eqnarray} \tag{ A.5 }$$

with a certain threshold probability ${p}_{\mathrm{swap}}^{\mathrm{th}.}$, e.g. ${p}_{\mathrm{swap}}^{\mathrm{th}.}\,=\,0.1$ or 0.05. In this way, one avoids the ${\beta }_{m}$ shifting unnecessarily close to each other and prevents very long time scales in which replicas oscillate only between two neighboring inverse temperatures.

$$\begin{eqnarray}{\boldsymbol{E}}={{\boldsymbol{E}}}_{t}\equiv \left(\begin{array}{ccccc}{x}_{1} & {x}_{3} & 0 & 0 & 0\\ {x}_{3} & {x}_{2} & {x}_{4} & 0 & 0\\ 0 & {x}_{4} & 0 & {x}_{4} & 0\\ 0 & 0 & {x}_{4} & -\,{x}_{2} & {x}_{3}\\ 0 & 0 & 0 & {x}_{3} & -\,{x}_{1}\end{array}\right)\end{eqnarray} \tag{ B.1 }$$

$$\begin{eqnarray*}{{\boldsymbol{\Gamma }}}^{(1)}=\left(\begin{array}{ccccc}{x}_{5} & {x}_{9} & 0 & {x}_{11} & {x}_{12}\\ {x}_{9}^{* } & {x}_{6} & 0 & {x}_{13} & {x}_{14}\\ 0 & 0 & 0 & 0 & 0\\ {x}_{11}^{* } & {x}_{13}^{* } & 0 & {x}_{7} & {x}_{10}\\ {x}_{12}^{* } & {x}_{14}^{* } & 0 & {x}_{10}^{* } & {x}_{8}\end{array}\right)\end{eqnarray*}$$

The parameters ${x}_{9}$ to ${x}_{14}$ can be complex. Therefore, it is straightforward to see that, for general N_B the number of independent (real) parameters is $C({N}_{{\rm{B}}})=\tfrac{{N}_{{\rm{B}}}}{2}({N}_{{\rm{B}}}+3)$ for the real case and $C({N}_{{\rm{B}}})={N}_{{\rm{B}}}({N}_{{\rm{B}}}+1)$ for the complex case.

$$\begin{eqnarray*}&&{\boldsymbol{E}}={{\boldsymbol{E}}}_{t}\end{eqnarray*}$$

$$\begin{eqnarray*}{{\boldsymbol{\Gamma }}}^{(1)}=\left(\begin{array}{ccccc}{x}_{5} & {x}_{9} & 0 & 0 & 0\\ {x}_{9} & {x}_{6} & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & {x}_{7} & {x}_{10}\\ 0 & 0 & 0 & {x}_{10} & {x}_{8}\end{array}\right)\end{eqnarray*}$$

so in general $C({N}_{{\rm{B}}})=3{N}_{{\rm{B}}}-2$.

$$\begin{eqnarray*}&&{\boldsymbol{E}}={{\boldsymbol{E}}}_{t}\end{eqnarray*}$$

$$\begin{eqnarray*}{{\boldsymbol{\Gamma }}}^{(1)}=\left(\begin{array}{ccccc}{x}_{5} & 0 & 0 & 0 & 0\\ 0 & {x}_{6} & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & {x}_{7} & 0\\ 0 & 0 & 0 & 0 & {x}_{8}\end{array}\right)\end{eqnarray*}$$

Here, $C({N}_{{\rm{B}}})=2{N}_{{\rm{B}}}$.

$$\begin{eqnarray*}{\boldsymbol{E}}=\left(\begin{array}{ccccc}{x}_{1} & 0 & {x}_{3} & 0 & 0\\ 0 & {x}_{2} & {x}_{4} & 0 & 0\\ {x}_{3} & {x}_{4} & 0 & {x}_{4} & {x}_{3}\\ 0 & 0 & {x}_{4} & -\,{x}_{2} & 0\\ 0 & 0 & {x}_{3} & 0 & -\,{x}_{1}\end{array}\right)\end{eqnarray*}$$

$$\begin{eqnarray*}{{\boldsymbol{\Gamma }}}^{(1)}=\left(\begin{array}{ccccc}{x}_{5} & 0 & 0 & 0 & 0\\ 0 & {x}_{6} & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & {x}_{7} & 0\\ 0 & 0 & 0 & 0 & {x}_{8}\end{array}\right)\end{eqnarray*}$$

Also here $C({N}_{{\rm{B}}})=2{N}_{{\rm{B}}}$.

Remember, here the impurity is on i = 1.

$$\begin{eqnarray*}{\boldsymbol{E}}=\left(\begin{array}{ccccc}0 & {x}_{1} & 0 & 0 & 0\\ {x}_{1} & 0 & {x}_{2} & 0 & 0\\ 0 & {x}_{2} & 0 & {x}_{3} & 0\\ 0 & 0 & {x}_{3} & 0 & {x}_{4}\\ 0 & 0 & 0 & {x}_{4} & 0\end{array}\right)\end{eqnarray*}$$

$$\begin{eqnarray*}{{\boldsymbol{\Gamma }}}^{(1)}=\left(\begin{array}{ccccc}0 & 0 & 0 & 0 & 0\\ 0 & {x}_{5} & {x}_{9} & 0 & 0\\ 0 & {x}_{9} & {x}_{6} & {x}_{10} & 0\\ 0 & 0 & {x}_{10} & {x}_{7} & {x}_{11}\\ 0 & 0 & 0 & {x}_{11} & {x}_{8}\end{array}\right)\end{eqnarray*}$$

In this case, $C({N}_{{\rm{B}}})=3{N}_{{\rm{B}}}-1$.