Exact real-time dynamics of single-impurity Anderson model from a single-spin hybridization-expansion

In this work we introduce a modified real-time continuous-time hybridization-expansion quantum Monte Carlo solver for a time-dependent single-orbital Anderson impurity model: CT-1/2-HYB-QMC. In the proposed method the diagrammatic expansion is performed only for one out of the two spin channels, while the resulting effective single-particle problem for the other spin is solved semi-analytically for each expansion diagram. CT-1/2-HYB-QMC alleviates the dynamical sign problem by reducing the order of sampled diagrams and makes it possible to reach twice as long time scales in comparison to the standard CT-HYB method. We illustrate the new solver by calculating an electric current through impurity in paramagnetic and spin-polarized cases.

Here we propose a method that alleviates the dynamical sign problem in case of the singleorbital Anderson impurity model (AIM). It is based on the continuous-time hybridizationexpansion QMC (CT-HYB-QMC [24]) and can be viewed as a kind of a bold-line QMC. In our method, called CT-1/2-HYB, only the diagrams for one out of two spins need to be sampled (hence the name). The contribution from the other spin can be summed up semianalytically for each expansion term. The advantage of such an approach, which to our best knowledge has not yet been applied neither to equilibrium nor non-equilibrium calculations, is that the Monte Carlo average expansion order is decreased leading to the reduction of the sign problem. The problem is however that the semi-analytic evaluation of dressed diagrams for a continuous-bath is a non-trivial task. The solution we come up with in this work is to discretize the bath, which limits the spectral resolution and the timescale obtainable within our method.
The paper is organized as follows. In Section 2 we describe CT-1/2-HYB method in detail introducing two possible ways of semi-analytical evaluation of the bold diagrams. In Section 3 we present the strategy for discretizing the bath and results for a current through impurity in a paramagnetic and a spin-polarized case, as well as discuss the computational performance of the algorithm. In section 4 we summarize our findings and draw conclusions.

Single-spin hybridization-expansion
The system for which CT-1/2-HYB is designed is the single-orbital AIM with general timedependent parameters. Since the test case we consider in this paper is electric current calculation, it is convenient to introduce a setup with two baths, where each bath (lead) is labelled by index α. The Hamiltonian is composed of the local, bath and hybridization parts and The contour starts at 0 + at the upper real branch, goes to t max , then goes back to 0 − along the lower real branch and eventually reaches −iβ along the imaginary axis. Full and empty circles denote virtual creation and annihilation processes of impurity electrons. Bold blue line denotes the semi-analytic dressing of spin-down expansion diagrams by all the spin-up processes. reads The method of solution is based on the perturbative expansion of the dynamical partition function Z with respect to hybridization strength on the Keldysh contour C with an additional imaginary time branch, defined by inverse temperature β and maximum evolution time t max (Fig. 1). The imaginary time branch enables us to start the time-evolution from an equilibrium thermal state at t = 0. An accessible introduction to this formalism containing imaginary-and real-time evolution can be found e.g. in [25].
Contrary to CT-HYB, the expansion is performed only for one spin species. We will be referring to spin-down as the expansion channel, which is an arbitrary choice made for a convenience of the discussion. In the following we work in the Schrödinger picture: contourtime arguments of operators serve only as contour labels (they do not denote Heisenberg-like evolution). Thus, we expand Z in H hyb↓ obtaining The contribution of each diagram, specified by expansion order k along with sets of impurityelectron creation {t c m } and annihilation times {t a m }, can be represented as a product of the spin-down bath weight where C † ↓ (t) ≡ p V p↓ (t)c † p↓ , and the local weight dressed by spin-up hybridization The spin-down bath weight can be readily evaluated following the procedure known from CT-HYB method [1]. Defining the hybridization function where g αpσ (t, t ) is a free-electron Green function [25] we obtain where Z bath↓ = Tr c ↓ T C e −i C dtH bath↓ (t) is a diagram-independent constant. Namely, for each diagram one needs to evaluate a determinant of a k × k spin-down hybridization matrix.

Dressed local weight from discretized bath
The novel aspect of CT-1/2-HYB is the presence of the spin-up hybridization dressing in the local weight, Eq. (7), which incorporates both local and spin-up bath degrees of freedom. This auxiliary system is governed by the Hamiltonian H 0 = H loc + H bath↑ + H hyb↑ and is effectively non-interacting since H 0 conserves n ↓ = d † ↓ d ↓ . The contour-time dependence of n ↓ is specified by the the locations of operators d † ↓ : {t c m } and d ↓ : {t a n }. Because n ↓ can equal only 0 or 1, the only non-vanishing contributions to the dressed local trace come from the following orderings of the spin-down operators: In the following we assume that only those two orderings of contour-times {t c m }, {t a n } are allowed. Note that for k = 0 both contributions have to be summed, corresponding to constant n ↓ = 0 or 1 on the whole contour.
The dressed local weight can be represented as The factor ϕ takes into account the permutation sign and the phases resulting from tracing over impurity spin-down occupations. If we denote the effective contour-time-dependent spindown impurity occupation by n ↓ (t) ≡ n ↓ (t; {t c m }, {t a n }) then The determinant part represents the trace of the evolution operator of the spin-up impurity and the spin-up bath subject to an effective contour-time-dependent single-particle Hamilto- The single particle evolution operator U is given by and the relation between the trace of the many-body evolution operator and the singleparticle evolution operator is provided by the identity valid for arbitrary fermions f and a non-interacting Hamiltonian h(t) Note that this identity is also used in the determinant QMC [26], which is similarly based on the summation of determinants of matrices in the single particle basis. It is known that care must be taken while evaluating determinants of the form of Eq. (15) for large β due to exponentially large and small matrix elements. We followed the advice from [27] to stabilize the necessary matrix multiplications and evaluation of determinants. Because the contour-time dependence of n ↑ (t) is step-wise constant we can take advantage of a convenient segment picture [1] to obtain a more practical representation of φ and U. We call a segment a time interval during which n ↓ = 1, i.e. after a creation time and before an annihilation time. Then we have and where u n ↓ is a single-particle time-evolution matrix

Discretizing contour-time instead of bath
Another approach to evaluate the dressed local weight is to discretize the contour-time instead of the bath. Within the path integral-formalism one can integrate out the spin-up bath from Eq. (7) and work directly with the hybridization function ∆ ↑ defined by Eq. (8). After calculating traces over the bath and spin-down impurity degrees of freedom we obtain where ϕ is the same as in Eq. (12) and (16), Z bath↑ is a constant and Z d↑ is given by the path integral with and where n ↓ (t) is determined by the sets of contour-times {t c m }, {t a n }. When G ↑ is understood as a matrix in time domain, the path integral Z d↑ is given by the formula [28] Because it is undesirable to discretize the time derivative operator we can recast this formula using Eq. (20) where • denotes contour-time convolution and z d↑ is a generalized partition function for a system with time-local action which can be explicitly evaluated For each diagram a generalized Green function G 0↑ should be calculated as a solution to Eq. (21). However, this has to be done for several contour-time spacings ∆t as an extrapolation of the result down to ∆t = 0 is required. Provided the convolution in Eq. (23) is done carefully as described in [29] the error of the discrete approximation should scale as ∆t 2 .
We have decided not to follow this method due to its high computational cost. Nevertheless, this approach might be of a future interest since it allows one for an introduction of a retarded interaction U ↑↓ (t, t ) by adding an additional term −U ↑↓ (t, t )n ↓ (t) to the left hand side of Eq. (20), which is only a small numerical effort.

Measurement of observables
The partition function Z from Eq. (5) can be viewed as a sum of the weights w(C) of diagrams C = {k; {t c m }, {t a n }} for all expansion orders k and contour-times {t c m }, {t a n } The space of diagrams C is sampled within a Markov-chain Monte Carlo simulation. We have implemented the segment removal and addition moves, which guarantee the ergodicity of the sampling algorithm, and the shift moves (shifting the beginning or the end of a segment), which decrease the autocorrelation time. A detailed description of those moves along with Metropolis-Hastings acceptance rates can be found in [1]. We use fast updates for determinants of the hybridization matrices (bath weight) but calculate the dressed local weight explicitly for each Monte Carlo diagram. The diagram weights discussed above can be in general any complex numbers. To work around this issue, we sample the diagrams with probabilities given by the absolute value of weights while simultaneously measuring their average (complex) sign. The expectation value of an observable A is given by where A(C) is the expectation value of A in diagram C, the complex sign is defined as and the Monte Carlo average follows the prescription It is known that CT-HYB-QMC for a single-impurity Anderson model does not suffer from the fermionic sign problem, i.e. the average sign on the imaginary axis is always one. However, the presence of real-time branches causes the dynamical sign problem. For each diagram having operators on the real branch there exists another one which has exactly the same weight but the opposite sign. This results from the fact that operator closest to t max on the real axis can be moved freely from the upper to the lower branch (or conversely) and this operation changes only the sign of the diagram. This means that the average complex sign sgn is a real quantity which equals the probability of having no operators on real-time branches. Practically, the average sign decreases exponentially with t max which effectively prohibits an accurate evaluation of observables using Eq. (26) for large t max . To overcome the inaccuracy of the results we have to scale the total number of Monte Carlo samples exponentially with t max until it is feasible.
In the following we present explicit measurement formulas for certain observables of interest.
The time-dependent impurity spin-down occupation is the average of the diagram-dependent function n ↓ (t; {t c m }, {t a n }) where t can lie either on the upper or the lower real-time branch The measurement of spin-up observables requires the evaluation of the elements of occupation matrix for the system of coupled spin-up impurity and bath. Namely, it can be shown This result follows from the general formula analogous to the one found in [26] Tr The spin σ electron current from lead α to the impurity is given by which can be obtained by solving the Heisenberg equation of motion for the impurity particle number operator d † σ d σ and employing the continuity equation While evaluation of the expectation value of the current operator requires a reference solution within the standard CT-HYB algorithm (see e.g. Eq. (60) in [5]), it is possible to extract the full spin-up single-particle occupation matrix directly in CT-1/2-HYB and by that all the relevant matrix elements for the spin-up current. Similarly to (30) where αp is the index of the column corresponding to pth electronic level from bath α. In practice, one measures the whole spin-up single-particle occupation matrix as defined in (30) and (35), from which the observables of interest can be constructed.

Discrete representation of the bath
Before presenting the results, we show in detail the parameters of the bath we have chosen for our simulation. We want to model a bath defined by the following time-and spin-independent coupling function which represents a flat band of half-width D = 4 with a soft cutoff controlled by parameter ν = 3, coupled to the impurity with the coupling strength Γ = 1. However, within the discrete bath approach, only a peak-structured coupling function can be constructed In order to approximate the continuous bath-impurity coupling function, Eq. (36), by a discrete one, Eq. (37), we apply the equal weight method from Ref. [30], obtaining equal hybridization parameters V p The resulting discrete representations of the bath for several values of bath size N are presented in Fig. 2.

Current through impurity
We measure the electric current in a simulation of a voltage quench from φ(t = 0) = 0 to φ(t > 0) = 2. In order to analyze the convergence of result with increasing bath size, we perform QMC calculations for several numbers N of bath sites. We set inverse temperature β = 5 and Coulomb interaction on the impurity U = 4. We consider the half-filled case, such that the (time-independent) local energy level In Fig. 3 we plot the spin-up electric current through impurity as a function of time. Due to the particle-hole symmetry the left hand side of Eq. (34) vanishes, thus currents from the left and right are equal up to the sign. Also, in the paramagnetic case the total current can be recovered by multiplication by factor of 2. For reference we also provide a result of an exact-diagonalization calculation for N = 6 which confirms the validity of our method.
As also observed in a similar study [18], the current relatively quickly (around 1 unit of Γ −1 ) reaches its steady state value. We note a quick convergence of the result with N . The difference between results for N = 18 and N = 24 is only minimal. It is instructive to look at the comparison of the average expansion order and the average Monte Carlo sign between CT-1/2-HYB and CT-HYB, which we present in Fig. 4 for N = 12 (the dependence of those quantities on N is negligible). In the paramagnetic case CT-1/2-HYB leads to the reduction of the average expansion order by a factor of two since only spin-down diagrams have to be sampled while spin-down and spin-up diagrams are contributing to Z with equal weights. Noteworthy, the expansion order grows linearly with time and the rate of growth is smaller by a factor of 2. According to the established exponential decrease of the average sign with the growing expansion order [4,5] twice as long time scales can be reached within CT-1/2-HYB method as compared to CT-HYB, when the smallest average sign feasible for obtaining the result is fixed.

Effect of magnetic field
A valid research question is whether the CT-1/2-HYB method can lead to an even better performance improvement over CT-HYB in presence of a magnetic field breaking the SU(2) symmetry of the model. In such a case, the average expansion orders in the two spin channels can differ, and hence one can choose to expand Z in the channel with the lower expansion order. In the following, we modify the definition of the local energies E dσ such that it takes into account non-zero chemical potential µ and magnetic field B  The maximum spin-down average order occurs when the spin-down impurity occupation is close to 0.5, which is evident from Fig.  5(b). For further analysis, it helps to realize that the following relation holds where k ↑ is a hypothetical average expansion order in case the hybridization expansion was performed in the spin-up channel. Then, we define the improvement factor η as where the sum k ↓ + k ↑ is a hypothetical average expansion order in the standard CT-HYB algorithm. We plot the equilibrium improvement factor as a function of |B| in Fig. 5(c).
In the paramagnetic case the improvement factor η is always 2. For µ = 0 the improvement factor also equals 2 and does not depend on B, which results from the particle-hole symmetry.
Only by simultaneous introduction of non-zero µ and B the average expansion orders can be made different in the two spin channels, leading to η > 2.  Fig. (a) can be understood as spin-down currents for B = −2, just as Fig. (b) can be viewed as spin-down currents for B = 2 While Fig. 5(c) proves that CT-1/2-HYB can more than halve the average expansion order on the imaginary axis, it is not yet clear if the same applies to the real-time evolution. Ultimately, it is the growth of the expansion order with the evolution time that is causing the dynamical sign problem. That is why we decided to estimate this quantity for the calculation of time-dependent current in the presence of chemical potential and magnetic field. We choose µ = 2 and B = 2, −2, which is a case where an improvement factor in equilibrium of around 2.5 can be achieved.
In the polarized case currents from the right and from the left bath are different. Fig. 6 presents the left and the right spin-up currents in the presence of magnetic fields B = −2, 2.
The spin-up current is lower for B = 2 since then the spin-up impurity level lies almost at the edge of the bath which hinders transport dynamics. Fig. 7 shows the growth of the expansion order with time for the paramagnetic case compared with the two polarized cases. Contrary to the imaginary-time calculation, the expansion order grows at almost the same rate irrespective of chemical potential and magnetic field. Noteworthy, for B = −2 the growth rate of the average expansion order is slightly higher than for other cases, even though the initial average order (for t = 0) is the lowest in this case.
The different behaviors of the average expansion orders in imaginary and real-time expansions were noted previously in [5]. The average expansion order on the Keldysh contour has no interpretation as a quantity proportional to the kinetic energy and is thus only weakly influenced by the position of the impurity level relative to the conduction bath. Although the initial (equilibrium) average expansion order depends on the magnetic field B and chemical potential µ, the growth rate of the expansion order is hardly influenced

Computational performance
The biggest computational burden of the method comes from the need to overcome the dynamical sign problem. For that reason we have scaled the number of Monte Carlo steps exponentially with t max in order to keep the error approximately constant for all times. We have been able to generate enough samples such that the statistical errors are smaller than size of markers on all the plots presented in this work. The dependence of the average sign on the physical parameters is the same as in the standard CT-HYB method on a Keldysh contour. Namely, keeping the hybridization strength as energy unit, the sign is largely unaffected by interaction strength, chemical potential or magnetic field but decreases with bath bandwidth [4]. Now we discuss the computational complexity of diagram evaluations. For the discretized bath approach the number of bath sites N plays the most important role. Generally, the expected scaling of this approach is O( k ↑ N 3 ) and results from multiplications single-particle evolution matrices of size (N + 1) × (N + 1), whose average number scales linearly with k ↑ . The actual scaling can deviate from this prediction as certain implementation-dependent optimizations might be possible.
For the discretized time approach the number of contour time steps N t determines the computational scaling. Both the solution of the Dyson equation [25] and the evaluation of the determinant in time-space will scale as O(N 3 t ). Moreover, it takes on average k ↑ Monte Carlo moves to decorrelate subsequent measurements. This implies an additional scaling factor of k ↑ [13]. The average order k ↓ is linear in inverse temperature β and t max : k ↓ = c β β + c t t max . Leaving the dynamical sign problem aside and neglecting the cost of bath weight evaluation, the overall theoretical complexity of the discretized bath approach is O( k ↓ 2 N 3 ) and that of the discretized contour-time approach is O( k ↓ N 3 t ) (where N t is the number of time-steps, including those on the imaginary branch).
The computational scaling of our implementation of diagram evaluations can be learned from Fig. 8. The benchmark problem is the evaluation of the paramagnetic current from It is also constructive to understand the scaling of the necessary computational resources with growing N . In Fig. 8(b) we present the CPU time per MC step as a function of N for a fixed t = 1. The fitted computational complexity is O(aN 2 + bN 3 ) with a b. The most general version of our implementation of CT-1/2-HYB involves a precomputation of all u 0 (t, t ) and u 1 (t, t ) on a two-dimensional time-grid and subsequent interpolations during a Monte Carlo run. Unfortunately, this approach can quickly lead to very high memory demands, especially for large N . However, when the time-dependence of Hamiltonian parameters is step-wise, which is the case for the voltage quench, the time-evolution generated by the operators u n ↓ can be performed in their eigenbases. For this purpose, one needs to store only a small set of eigenenergies and rotation matrices, which not only solves the problem of high memory requirement but also avoids the interpolation-related numerical inaccuracy and presumably improves performance. We have applied this strategy for the calculations presented here, which possibly explains the O(N 2 ) contribution to the scaling due to presence of diagonal matrices.

Conclusion
CT-1/2-HYB-QMC allows one to solve the single-orbital Anderson impurity model exactly on twice as longer time scales than standard CT-HYB-QMC due to a reduction of the average expansion order. The growth rate of the average expansion order with t max is independent of the location of impurity spin-up and spin-down levels in the considered time-evolution after the voltage quench, which rules out an even better performance gain in cases with non-zero chemical potential and magnetic field.
It is possible to extend this approach to multi-orbital impurity models with density-density interactions. If all the orbitals are interacting, the improvement factor will diminish with the number of fermionic flavors (combined spin and orbital numbers) n as n n−1 . This results from the fact that the semi-analytical evaluation of the diagrams will be possible for only one out of n flavors. Moreover, in this case the method can only be applied to Hamiltonians conserving orbital and spin quantum numbers, which poses a serious limitation to the modelling of realistic materials. However, if the interaction matrix is sparse, also some single-particle mixing terms can be allowed between non-interacting flavors, while the improvement factor can be better than n n−1 as more flavors can be treated semi-analytically. We note that the cost of the CT-HYB inchworm algorithm [17] increases exponentially with number of orbitals due to an exponential increase in the size of many-body Hilbert space, which does not depend on the sparsity of interactions. On the other hand, in CT-1/2-HYB the increase in computational cost with the number of orbitals is in principle only polynomial. The main challenge of a real-time QMC is however the dynamical sign problem, which still leads to an exponential CPU-time scaling with t max in CT-1/2-HYB and only to a polynomial scaling in case of inchworm algorithm. Even though CT-1/2-HYB-QMC may significantly extend accessible timescales for multi-orbital models with sparse interactions, it is still subject for future analysis to understand whether it will be enough to study any phenomena of interest.
It is worthwhile to consider whether the CT-1/2-HYB could be combined with bold or inchworm QMC. Their strength is that dressed many-body propagators can be used to evaluate Monte Carlo diagrams. Those propagators are actually of no use in CT-1/2-HYB as an attempt to include many-body corrections to CT-1/2-HYB would undermine the basic assumption about the effective single-particle dynamics. Nevertheless, the inchworm expansion could be done around propagators dressed by spin-up processes in the spirit of CT-1/2-HYB. In contrast to pure CT-1/2-HYB, some spin-up diagrams would be missed in such an inchworm procedure and would still have to be sampled separately. It is unclear at the moment how the efficiency of the inchworm algorithm would be affected by use of dressed propagators since no such attempts have been reported.
An interesting single-orbital application of CT-1/2-HYB method as an impurity solver for DMFT could be investigation of the cross-over behavior between Falicov-Kimball and Hubbard models, by tuning the spin-down hopping from zero to the value of the spin-up hopping. Our method is especially suitable for studying this problem since the zeroth order term of 1/2-HYB expansion corresponds to Falicov-Kimball-like impurity.