Quantum impurity models using superpositions of fermionic Gaussian states: Practical methods and applications

The coherent superposition of non-orthogonal fermionic Gaussian states has been shown to be an efficient approximation to the ground states of quantum impurity problems [Bravyi and Gosset,Comm. Math. Phys.,356 451 (2017)]. We present a practical approach for performing a variational calculation based on such states. Our method is based on approximate imaginary-time equations of motion that decouple the dynamics of each Gaussian state forming the ansatz. It is independent of the lattice connectivity of the model and the implementation is highly parallelizable. To benchmark our variational method, we calculate the spin-spin correlation function and R\'enyi entanglement entropy of an Anderson impurity, allowing us to identify the screening cloud and compare to density matrix renormalization group calculations. Secondly, we study the screening cloud of the two-channel Kondo model, a problem difficult to tackle using existing numerical tools.


I. INTRODUCTION
Quantum impurity models -systems of a few strongly interacting degrees of freedom coupled to a large bath of noninteracting fermions -constitute an important class of problems in condensed matter physics. Despite the small number of interacting modes involved, this class of problems can exhibit rich many-body physics phenomena. The archetypical phenomenon is the Kondo effect, where even weak interactions can lead to strong nonperturbative corrections to the ground state [1]. Such models also appear as effective models in many embedding methods, such as dynamical mean-field theory [2], that solve extended quantum many-body systems by approximately mapping them to quantum impurity problems.
Over the years, various numerical methods have been developed to tackle quantum impurity problems. A particularly successful approach is Wilson's numerical renormalization group (NRG) [1] and its extensions [3] which have allowed to study this class of problems in the thermodynamic limit. A related set of variational methods, based on the density matrix renormalization group (DMRG) [4] has also been used extensively and compared to NRG [5,6]. Finally, Quantum Monte Carlo methods have been successfully applied to systems where the sign-problem is mild, see e.g. Refs. 7 and 8. Despite these methods being very powerful, they each come with their limitations. For example, NRG is limited by an exponential scaling in the number of degrees of freedom of the impurity and the number of channels in the non-interacting bath; DMRG scales more favorably in the size of the impurity, but (since it does not exploit the non-interacting nature of the bath) retains an exponential scaling with the number of channels in the bath, rendering it very challenging to study, e.g., mesoscopic problems with several leads. Quantum Monte Carlo methods, on the other hand, typically suffer from severe sign problems for multi-orbital systems.
A natural question is whether a well-chosen class of variational states could exploit the structure of quan-tum impurity models to circumvent the limitations of these established approaches. Recently, this question was affirmatively answered by proving that the ground states of quantum impurity problems can be approximated by a superposition of non-orthogonal fermionic Gaussian states [9]. This is obviously the case when the number of states in the superposition -which we will refer to as the rank of the ansatz -grows exponentially with the full system size, as the states then form a complete many-body basis. More interestingly, the rigorous mathematical bounds of Ref. [9] demonstrate that the minimal rank to obtain a good approximation of the ground state scales only with the size of the impurity and the desired precision, while being independent of the size of the bath. This superposition of Gaussians (SGS) ansatz can be seen as a generalization of the generalized Hartree-Fock (GHF) method [10,11], which aims to find the approximate ground state of a system using a variational minimization over the field of fermionic Gaussian states. While exact for a noninteracting system, GHF corresponds to a mean-field approximation for interacting systems. In the context of quantum chemistry, approaches related to our work are known as multi-component Hartree-Fock-Bogoliubov methods [12,13].
Having chosen this ansatz, the challenge is to device practical algorithms to perform numerically efficient computations. Here, we focus on the problem of finding the lowest-energy state within the variational manifold. Multiple generic approaches exist for performing the energy minimization within a variational space, such as gradient descent and imaginary time evolution [14,15]. In particular, formal solutions for related ansatzes were originally developed under the name of resonating Hartree Fock [16][17][18]. However, due to the large number of variational parameters in the SGS ansatz and the presence of several non-linear constraints, their numerical implementation can become prohibitively costly for large systems.
In this work, we propose a simpler and numerically less costly path towards energy minimization within the arXiv:2105.01088v1 [cond-mat.str-el] 3 May 2021 variational manifold based on several key approximations to the imaginary time equations of motion. First, at each step we project the dynamics onto the subspace orthogonal to the one spanned by the current set of Gaussian states forming the SGS. Furthermore, we alternate the evolution of the coefficients of the coherent superposition of states and the (normalized) Gaussian states themselves. This allows us to decouple the equations of motion for each Gaussian state at each step in the evolution. While the projection of these equation of motions onto the variational manifold does not exactly correspond to imaginary time evolution, we show that under this evolution, the energy is non-increasing. The variational state therefore converges to a local energy minimum within the manifold and can thus be used to study ground state properties of quantum models (we note that a guarantee on convergence to a global minimum typically cannot be given for variational algorithms).
To illustrate the power of the method, we apply it to two canonical impurity models: the single-impurity Anderson model [19] and the two-channel Kondo effect [20]. The former has been studied using a variety of methods and is well-understood both analytically and numerically, thus allowing us to confirm the validity of our method. We find that using comparable computational resources, our method is able to achieve an error in the ground state energy that is about one order of magnitude better than DMRG. The two-channel Kondo model, on the other hand, is much more challenging to study numerically, and real-space correlation functions for fermionic leads had previously eluded numerical simulations. Instead, prior numerical studies have required either mapping to related problems in the same universality class [21] or focusing on quantities that can be calculated using the local dynamics of the impurity [22].
The remainder of this work is structured as follow. In Sec. II, we set the notation and describe the structure of the SGS ansatz. In Sec. III, we describe a generic minimization procedure for finding a variational approximation to the ground state and its numerical implementation. In Sec. IV, as a first demonstration of the method, we study the screening cloud of a single impurity Anderson model and benchmark our results using DMRG. Finally, in Sec. V we extend the calculations of the previous section by considering the two-channel Kondo model.

II. ANSATZ AND PROBLEM STRUCTURE
We start by describing the structure of fermionic quantum impurity models and of the SGS ansatz. We also introduce the covariance matrix formalism for Gaussian states which will be used throughout this work.

A. Generic quantum impurity model
We consider a lattice model of N fermionic degree of freedoms with Hamiltonian H = H 2 + H 4 . We choose to work within a formalism of Majorana operators, noting that any fermion problem (both with and without particle-number conservation) can be rewritten in this form. The noninteracting part of the Hamiltonian is given by where A is a real and skew-symmetric matrix, and c k = c † k is a set of 2N Majorana operators that obey standard anticommuation relations {c j c k } = 2δ j,k . The interacting part of the Hamiltonian reads where the rank-4 tensor U is skew-symmetric with respect to the exchange of any neighboring indices. The interaction involves at most M N distinct Majorana operators making U a sparse tensor. We remark that no assumption with regards to the lattice connectivity is made in our model definition. Although we focus here and below on quartic interaction terms, there is no fundamental limitation to including interaction terms involving a larger number of operators.
The ground state of the quadratic Hamiltonian H 2 will be a Slater determinant if H 2 conserves the number of particles or, more generally, a fermionic Gaussian state [23] in the case where the U (1) symmetry is broken and only the parity of the number of particles is preserved. As Slater determinants constitute a subset of the fermionic Gaussian states, we will focus our discussion on Gaussian states. Although there might be a slight numerical overhead associated with working in this enlarged class of states, it has the advantage of naturally allowing the treatment of (mean-field) superconductivity.

B. Covariance matrix formalism
Our description and usage of the covariance matrix formalism closely follows Ref. [9]. Any fermionic Gaussian state |φ obeys a Wick theorem and thus can be fully described by a covariance matrix (CM) where k, l ∈ 1 . . . 2N . This matrix is real and skewsymmetric by construction. The expectation value of any product of Majorana operators can then be calculating as the Pfaffian of a submatrix of Γ. For a normalized pure state the elements of the CM are subject to the constraint Γ 2 = −1.
Since the covariance matrices Γ are matrices of expectation values, they are invariant under a gauge transformation |φ → e iθ |φ , with θ a real number. For calculations involving multiple Gaussian states |φ µ (µ = 1, 2, . . . ), it is often necessary to fix this gauge freedom. Following Ref. [9], this can be achieved by choosing a reference state |φ 0 and taking φ 0 |φ µ to be real and positive for all µ. Overlaps of Gaussian states and matrix elements can be obtained using the respective covariance matrices of the states [9].
C. Sum of Gaussian states ansatz Following Ref. [9], the variational ansatz considered in this work is formulated as 1 with {|φ µ } a set of nonorthogonal Gaussian states, λ µ complex scalar amplitude, and R the rank of the ansatz. This variational state is characterized by the set of covariance matrices and amplitudes {Γ µ , λ µ |µ = 1 . . . R}. This corresponds to O(RN 2 ) variational parameters, subject to two normalization constraints. First, the normalization of the variational state requires where we introduce G µ,ν = φ µ |φ ν the Gram matrix characterizing the overlap between states. Second, we take each Gaussian state to be pure and normalized leading to the constraint on each covariance matrix (Γ µ ) 2 = −1. The sets of parameters obeying these constraints form the variational manifold. Our aim is to find the state |ψ 0 that minimize the energy within this manifold. Before turning to the full minimization problem, we first consider the following simpler problem. Given a set of Gaussian states {|φ µ } we wish to find the amplitudes {λ µ } which minimize the energy E = ψ| H |ψ . These optimal amplitudes, leading to a normalized state with the lowest energy within the subspace, can be obtained by diagonalizing the Hamiltonian projected onto the subspace spanned by the set of Gaussian states. This leads to the generalized eigenvalue problem [9,16] hλ = EGλ (6) where h is the R × R matrix with elements h α,β = φ α | H |φ β . In the case of orthogonal states G = 1 and this reduces to a regular eigenvalue problem. This standard result will be at the core of the minimization approach introduced in Sec. III as we will alternate between updating the amplitudes by solving Eq. (6) and updating the covariance matrices assuming fixed amplitudes.

III. PROJECTED EQUATIONS OF MOTION FOR ENERGY MINIMIZATION
A generic approach for finding the ground state of a quantum system is imaginary time evolution. Starting from an initial state |ψ and evolving according to the imaginary time Schrödinger equation where E ψ = ψ| H |ψ , allows to reach the ground state in the τ → ∞ limit as long as the initial state has finite overlap with the ground state. In the case where |ψ is a variational state, the equations of motion must be projected back onto the part of the variational manifold orthogonal to |ψ in order to best approximate the dynamics of the system and preserve the norm of the state. Different methods were introduced to perform this projection [24]. In the case of the time-dependent variational principle (TDVP), this projection requires the inversion of the Gram matrix of the tangent states (obtained by taking the derivative of |ψ with respect to each variational parameter), which can be very large [14]. For the parametrization considered in this work, this does not appear to be a scalable approach as it would require repeated operations on matrices of dimension RN 2 ×RN 2 leading to an O(R 2 N 6 ) numerical complexity.
In this section, we instead derive simplified projections of the imaginary time equation of motion for the SGS states. Our approach can be understood as a parametrized energy descent. We derive equations of motion for the covariance matrices Γ µ as a function of an external parameter s such that the energy E(s) = ψ(s)| H |ψ(s) decreases monotically and converges to a local energy minima as s → ∞. As we do not pretend the followed approach to be sufficient to recover the system dynamics, we denote the evolution parameter as s to distinguish the resulting equations from imaginary time (denoted τ above).

A. Path of energy descent
As eluded to at the end of Sec. II, at any instant s, one can separate the Hilbert space in two instantaneous subspaces, where one (referred to as "parallel subspace" below) is spanned by the set {|φ µ (s) } of Gaussian states forming the ansatz, and the other is the orthogonal complement ("orthogonal subspace" below). The energy minimization in the parallel subspace is easily performed by choosing amplitudes satisfying Eq. (6) and we will thus focus here on the orthogonal subspace assuming fixed amplitudes.
We consider as a starting point the imaginary time Schrödinger equation. However, instead of projecting on a subspace orthogonal to the instantaneous state |ψ(s) leading to Eq. (7), we project onto a subspace orthogonal to the set {|φ µ (s) } of Gaussian states where we have introduced the projector The inverse Gram matrix in the above equation ensures that (Π ⊥ ) 2 = Π ⊥ and is necessary since the Gaussian states considered are generally nonorthogonal (but linearly independent, insuring that G is non-singular). This equation of motion preserves the norm of the state since Π ⊥ (s) |ψ(s) = 0 leading to ∂ s ψ(s)|ψ(s) = 0.
The projection in Eq. (8) constitutes a first approximation to the equations of motions and its justification is twofold. First, assuming optimal state amplitudes, any dynamics lowering the energy should be orthogonal to the instantaneous parallel subspace. However, the resulting equation of motion is only approximate due to the additional implicit projection onto the Gaussian state manifold. This projection will become explicit through the use of the Wick theorem in Sec. III B. Second, the projection to the orthogonal subspace ensures that the different Gaussian states do not collapse into to a single state. Indeed, without it, all states would collapse to a mean-field approximation of the ground state independently.
Inserting the definition of the ansatz into the projected equation of motion Eq. (8), we obtain In order to move forward, we decouple the equation of motions of the different Gaussian states by postulating that there exists an effective Hamiltonian B µ such that and which satisfies Eq. (10). Upon inspection, one can find that decouples the equations, with c µ a complex scalar. Taking c µ = [λ µλ * µ ] −1 , with the renormalized amplitudẽ λ µ = φ µ |ψ , satisfies Eq. (10) confirming the validity of the decoupling scheme. Although this is not a unique choice, this approach ensures a monotonic decrease of energy, i.e. ∂ s ψ| H |ψ ≤ 0, upon simultaneous integration of Eq.(11) for all states.
More generally, ensuring a path of energy descent upon the evolution of a given state |φ µ leads to the constraint Re c µλ * µ λ µ > 0. Hence, taking c µ = λ µ λ * µ /|λ µ λ µ |, possibly up to a real and positive multiplicative factor for each state, is sufficient to decrease the energy of the variational state.
This decoupling scheme is the second major approximation to the equations of motion used in this work. While exact for a generic many-body basis, the decoupling is approximate in the case where the evolution is projected on a constrained variational manifold. As is shown numerically in Secs. IV and V, the optimization nevertheless converges towards the ground state.

B. Equation of motion for the covariance matrices
In order for the above equations of motion to be useful, there must exist an efficient numerical implementation of them. We now derive the counterpart of Eq. (11) in the covariance matrix (CM) formalism.
Taking the derivative of the CM, as defined in Eq. (3), with respect to the evolution parameter s we obtain (taking k = l) where c.c. denotes the complex conjugate. Inserting Eq. (11) and the definition of Π ⊥ one obtains after some algebra In the special case of a rank 1 ansatz, Eq. (14) falls back onto the equation of motion for imaginary time evolution in the GHF approximation as derived for example in Ref. [11]. We now specialize to the quantum impurity model with quartic interacting defined in Sec. II A. In order to rewrite the differential equation purely as a matrix equation, we introduce the complex skew-symmetric matrices [9] which allows to compute easily matrix elements between different Gaussian states φ β | c k c l |φ α = iG β,α ∆ α,β k,l (k = l). More generally matrix elements involving n distinct Majorana operators are proportional to the Pfaffian of an n × n submatrix of ∆ α,β . Similarly to GHF, we introduce generalized Fock matrices which for α = β (∆ α,α = Γ α ) falls back on the standard Fock matrix for GHF [11].
With these definition, the matrix differential equation for each CM takes the form where the first and second line are respectively reminiscent of the equation of motion for real and imaginary time evolution of a Gaussian state [11]. The third line ensures normalization and cancels in the case R = 1.
We have verified numerically that Eq. (17) preserves the norm and the purity of the state which requires due to the constraint Γ 2 = −1 for a pure state. Equation (17) constitutes one of the main results of this work.

C. Numerical implementation
In order to find a good approximation of the ground state, we alternately evolve the covariance matrices for a small step s → s + δs and update the amplitudes by solving Eq. (6). The computational complexity of each of these steps is O(R 2 N 3 ) (assuming R N ). Since in each step, each Gaussian state can be evolved separately, one can easily parallelize over them and achieve a significant speedup. Appendix A presents additional details on the numerical implementation.
As an aside, we note that considering P interaction terms with m > 4 Majorana operators could, in the worst case, lead to additional operations of computational complexity O(P R 2 N 2 m 3 ). However, in the most relevant case where m, P N , this shouldn't affect the scaling of the overall computational complexity of the method. We note that the above estimate is an upper bound on the complexity of the right-and-side of Eq. (14) separately for each term of weight m and each matrix element. Taking advantage of the skewsymmetric matrix structure of the CM might allow further algorithmic improvement. For example, in the case of P = N 2 quadratic operators (m = 2) the naive es- In order to reduce the risk of converging to a local energy minima, we use in the remainder of this work the following strategy. Starting from the mean field (GHF) solution we gradually increase the rank of the ansatz by 1 after every n iterations. The additional Gaussian state is chosen by applying a random special orthogonal transformation Q to one of the current states of the ansatz with a high amplitude. For n sufficiently large and ||Q − 1|| F sufficiently small we have found numerically this approach to converge towards a good approximation of the ground state. This approach is needed due to our choice of optimal amplitudes at each step of the algorithm. If all states are added at once, there is a strong risk to converge to a solution where some of the amplitudes vanish, effectively leading to an approximate ground state of a lower rank.

IV. ANDERSON IMPURITY AND THE SCREENING CLOUD
As a benchmark of the method, we study realspace properties of the single impurity Anderson model (SIAM) [19]. This well-known model was previously studied numerically using a multitude of methods including NRG [25] and DMRG [26,27]. We will use the latter method to benchmark our results and confirm their validity.
We consider a 1D lattice where the SIAM Hamiltonian takes the form H = H 0 + H I , with describing a bath of L − 1 free fermions with hopping parameter t and chemical potential µ, and the impurity Hamiltonian Here, d † are fermionic creation operators, n r = d † r d r = σ d † r,σ d r,σ is the number operator, and when omitted, spin indices are summed. The hopping between the impurity (site 0) and the first site of the lead is t and U > 0 is a repulsive interaction. We focus on the particle-hole symmetric point at half-filling where d = −U/2 and µ = 0.
In the weak-coupling limit t < U , charge fluctuations are suppressed on the impurity site creating an effective spin 1/2 impurity. The Hamiltonian then maps to the Kondo model with effective coupling strength J = 8Γ/U , where Γ = (t 2 )/t is the broadening of the impurity energy level by the leads [28]. The ground state is then a singlet state formed by the impurity spin and the collective spin of a delocalized so-called Kondo cloud of electrons from the leads. One can associate to the cloud a length scale ξ K ∼ v F /T K where v F is the Fermi velocity and T K is the Kondo temperature [29].

A. Convergence and comparison to DMRG
As a first test of the convergence of the method, we compute the variance of the Hamiltonian in the variational ground state, where H = ψ| H |ψ . Figure 1(a) presents δ 2 H as a function of the ansatz rank and for different number of iterations of the variational minimization (see Sec. III C) for a system size L = 500. The variance δ 2 H should converge to zero as the variational state converges to an eigenstate of the system (δ 2 H = 0 for any eigenstate of H). We estimate in Fig. 1(b) the converged ground state energy E * using a linear extrapolation of the variational energy as δ 2 H → 0 for δ 2 H 2 × 10 −5 . As a second test of the validity of our variational SGS state results, Fig. 2 compares the energies obtained using the SGS ansatz to energies obtained using stateof-the-art DMRG simulations [30]. In both cases, we plot the results as a function of the refinement parameter of the ansatz, which is R for the SGS ansatz and the bond dimension m for DMRG. For the parameters considered, which were chosen to span similar computation times, the SGS-based method reaches lower energies than DMRG. along the lead where the spin operator is S(r) = 1 2 s,s d † r,s σ s,s d r,s . From a low-energy expansion of the bath operators [31] away from the impurity (k F r 1, with k F the Fermi wavevector), the correlation function is expected to be the sum of uniformly decaying and oscillating functions where, at zero temperature, C U (r) and C 2k F (r) are smoothly decaying functions. As we focus on a halffilled lattice where k F = π/2a (with a the lattice constant), the correlation function is the sum of uniform and staggered (cos(2k F r) → (−1) r ) contributions. Figure 3 compares the correlation functions computed using the approximate ground states obtained with SGS ansatz and DMRG. To plot more easily the highly oscillating function C(r), we introduce the correlation function on the even (e) and odd (o) sites, denoted as C e (r) and C o (r), respectively. For the odd sites ( Fig. 3(a)), the uniform and staggered parts of the correlation function are both negative, leading to a larger amplitude than the even sites ( Fig. 3(b)) where the contributions have opposite signs. Focusing first on the odd sites, the mean-field solution (R = 1, solid blue curve) differ qualitatively from the higher-precision DMRG results as expected. Modestly increasing the rank of the SGS ansatz, we recover the same behavior as high-precision DMRG. Small discrepancies between the two methods are observed far from the impurity where the amplitude of the correlation function is small. These differences are consistent with the expected precision of both methods and we expect that further increasing the values of m and R would reduce these differences. Similarly, for the even sites, both methods agree close to the impurity. Notably, the DMRG results shows a change of sign of C e (r) far from the impurity. This effect has been previously observed in other DMRG studies of the SIAM [26,27] and appears to disappear as the bond dimension of the MPS is increased. No such effect is observed for the SGS state, thus suggesting that this feature in the DMRG results is an artifact of the truncation of the MPS bond dimension.
As a final comparison between methods, we consider in Fig. 3

(c) the integrated correlation function
which allows characterization of the screening cloud of the impurity [26]. In the ground state, the total spin S 2 tot is expected to be zero for even L, leading to the sum rule Σ(L − 1) = 0. While this sum rule is not explicitly enforced, it is approximately respected and the violation converges towards zero as the rank of the ansatz is increased (Σ(L − 1) ≈ 2 × 10 −6 for R = 16).
To conclude this section, we use the SGS ansatz to study the screening cloud of the Anderson impurity for a large system of L = 1000 sites and different coupling strengths Γ = (t ) 2 /t. Figure 4(a) shows −C o (r) for three different couplings between the impurity and the bath. In the Kondo regime, the correlation function C o (r, Γ) is expected to collapse to a universal function C(r/ξ K ) through the relationC(r/ξ K ) = ξ K (Γ)C o (r, Γ) with ξ K (Γ) a coupling-dependent length scale. Figure 4(b) present the scaling collapse of ten different curves for the odd distances from the impurity C o (r). 2 As expected from previous numerical studies [25] and analytical calculations [31] the scaled correlation function decays following a power lawC(r/ξ K ) ∝ (r/ξ K ) −ν with a crossover in the exponent ν near r/ξ K = 1. From power-law fits, we extract the exponents ν ∼ 0.76 for r ξ K and ν ∼ 1.84 for r ξ K . Given the sensitivity of the results of the fit on the range of parameters considered, this results are consistent with the previously established scenarios of a crossover from ν = 1 to ν = 2 near r ∼ ξ K .
Finally, Fig. 4(c) presents the impurity screening length scale as extracted using three different methods. First, the parameters ξ K (Γ) were obtained from the scaling collapse of the correlation function C o (r, Γ) (blue disks). Second, as a comparison, we also plot the length scale ξ c (Γ) (orange crosses) at which the integrated correlation function falls below a given threshold: for a threshold parameter c ∈ (0, 1) [26]. Although simpler, this second method has the disadvantage of being sensitive to finite-size effects and convergence, as illustrated by the saturation of ξ 0.1 (orange crosses) for smaller couplings. These effects are reduced in the case of the scaling collapse approach as it takes into account the correlation function calculated at all odd sites. As a third method, we show the length scale ξ K,S (Γ) (green square) obtained from the scaling collapse of the impurity contribution to the entanglement entropy (near the impurity). This method will be described in the following section. As scaling collapse methods determine the screening length only up to a global prefactor, we scale the data sets such that all three methods result in the same length scale for 1/Γ = 5.
In the Kondo regime (Γ/U 1), we verify that the length scale follows the expected scaling where the prefactor is obtained through a fit (dashed black line) to the blue disks. In the regime of intermediate coupling strength all methods follow the expected exponential scaling. Away from this regime, in the case of weak coupling different level of sensitivity to finite size effect and convergence leads to underestimate the screening length compared to the expected exponential scaling. Similarly, in the strong coupling regime the mapping from the SIAM to the Kondo model is no longer valid and deviations are expected.

C. Impurity entropy
As a second probe of the SGS variational ground state, we consider the contribution of the Anderson impurity to the entanglement entropy. This quantity offers a different approach to study the screening cloud of an impurity [32][33][34].
We consider a bipartition of the sites in subsystems A and B and compute the second Rényi entropy where ρ A = Tr B |ψ ψ| is the reduced density matrix for subsystem A. Throughout, subsystem A is formed of the l first sites of the lead and the impurity. While in the case of a single Gaussian state the entropy can easily be computed from the decomposition in normal modes of the covariance matrix [35,36], the presence of coherences between the different Gaussian states forming the ansatz in the density matrix modify this calculation. In Appendix B we use the fermionic coherent state formalism [37] to derive expressions of S 2 (A, |ψ ) for an SGS state of arbitrary rank. Even though S 2 differs from the more common von Neumann entropy S 1 , many of the same universal properties can be extracted. In particular, using conformal field theory (CFT) calculations, the constant contribution of the boundary to the entropy was shown to be independent of the order of the Rényi entropy at criticality (see e.g. Refs. [38,39]). In addition the form of the leading corrections to the scaling of the entropy with subsystem size l due to irrelevant boundary operators was shown to be independent of the order n of the Rényi entropy [40,41].
As for the correlation function in Sec. IV B, for a lattice model at half-filling, the entanglement entropy is the sum of a uniform and a staggered contribution The uniform and staggered contributions to the entropy can be extracted using a local polynomial interpolation [33]. The impurity entropy S imp (l, Γ, L) is then obtained by the subtraction where S 0,U 2 (l, L) is the uniform part of the entropy in the absence of the impurity.
To make the calculation of S imp more explicit, Fig. 5(a) presents the raw data for a system of length L = 400 and Γ = 1/12. The corresponding impurity entropy S imp is the dashed red curve in Fig. 5(b). As expected, for small couplings (Γ 1) and close to the impurity (l ξ K ), S imp ∼ ln(2) indicating the entanglement of the impurity with the lead. Consistent with (c)  Fig. 4(c), the characteristic length scale over which the S imp decays increases as Γ is reduced.
To extract more quantitatively the screening length ξ K , we again consider a scaling collapse of the data. The impurity entropy was previously found to be amenable to scaling using [32] S imp (l, Γ,  Fig. 4(c) (green squares ξ K,S ). In the regime where the mapping from the SIAM to the Kondo model is valid (U/Γ 1) and ξ K L, we obtain a good agreement between the different method to extract the Kondo screening length.

V. TWO-CHANNEL KONDO EFFECT
As a second application of the variational method, we consider the two-channel Kondo (2CK) impurity model. We focus on the symmetric regime where the impurity is coupled with the same strength J to both channels. In this regime, the presence of a second channel of free electrons lead to very different ground state properties from the single channel case and the physics of the model includes an intermediate coupling fixed point with non-Fermi liquid behavior [20].
We focus on the ground state of the lattice version of the symmetric two-channel Kondo Hamiltonian where α = 1, 2 is the channel index, H J is the coupling between the impurity spin and the channels of free fermions and H BC accounts for the choice of boundary condition away from the impurity. We consider directly an antiferromagnetic spin-spin interaction of coupling strength J > 0 between the impurity spin and the first site of each channel: where S α (r) = 1 2 s,s d † r,α,s σ s,s d r,α,s is the fermionic spin operator. As the SGS ansatz is fermionic, we take the impurity to be a fermionic site with fixed singleoccupancy. 3 We also introduce the unitless coupling parameter g ≡ ρ 0 J with ρ 0 = 1/4t the density of state at the Fermi Level.
For finite-size systems, the two-channel Kondo model is known to exhibit important differences between the case where the total number of sites is even or odd [21]. In order to preserve the symmetry between the two channels and have an even number of sites (including the impurity), we introduce an additional site at the opposite end of the chain coupling the two channels leading to a total of 2L sites in the model. The 2CK model exhibits an intermediate coupling fixed point, as well as a duality between the weak and strong coupling regimes [42]. In order to explore this physics in the SGS variational ground state, we consider the spin-spin correlation function where r is the distance to the impurity and α is the channel index. We separate again the correlation function in a uniform and a staggered part K(r) = −K U (r) − cos(2k F r)K 2k F (r). Figure 6 shows the uniform part of the correlation function K U (r) for various unitless coupling strengths g = 0.15 − 1.5 on a symmetric logarithmic axis. Considering first K U close to the impurity, we observe a reduction of correlations as the coupling strength is increased. In absolute value, |K U (r)|, reaches a minimum for g * ∼ 0.7 followed by a sign change of K U (r) near the impurity. We associate this point where correlations are minimal (smallest screening cloud) as the intermediate coupling strength fixed point of the 2CK model. This result is in agreement with previous NRG results for the Kondo temperature T 2CK [42].
In the weak coupling regime (g g * ), the Kondo screening length of the 2CK model is to leading order ξ 2CK ∝ e 1/g [31,42]. In order to extract ξ 2CK (g), we use   a scaling collapse of K(r). We focus on the staggered part of the correlation function in the weak coupling regime where correlations are larger. To limit finite-size effects, we consider only correlations up to a distance L/2 from the impurity. Figure 7(a) presents a scaling collapse following the same procedure as in Sec. IV. For g = 0.15 − 0.75, all correlation functions can be collapsed on a universal scaling function exhibiting a crossover between two different power laws near r ∼ ξ 2CK . Using power law fits, we obtain the exponent ν = 1.04 near the impurity and ν = 1.58 far from the impurity (dotted and black dashed lines). An exponent ν < 2 is coherent with the expected non-Fermi liquid behavior of an overscreened multichannel Kondo model. In particular, for r ξ 2CK the scaling K 2k F (r) ∝ (r/ξ 2CK ) 1.5 was calculated in Ref. [31].
The extracted screening length ξ 2CK (g) is shown in Fig. 7(b) and compared to the expected exponential scaling (orange curve). Small deviations from the expected exponential behavior is expected for smaller values of g where finite-size effects become important (ξ 2CK L). This result is also in agreement with previous DRMG calculations for an effective spin chain representation of the symmetric 2CK model where ξ 2CK was extracted from a scaling collapse of the impurity entanglement entropy [21].

VI. CONCLUSIONS
In this work, we have developed a practical method for finding the variational ground state of an impurity problem using a coherent superposition of fermionic Gaussian states ansatz. The approach has a computational complexity O(R 2 N 3 ), thus scaling polynomially with the rank R of the ansatz and the number of fermionic modes N . Combined with the favorable scaling of the accuracy with rank R, as guaranteed by the results of Ref. [9], this gives a powerful new approach to study the ground state of quantum impurity models. In particular, the approach is independent of spatial locality or lattice connectivity, allowing for more flexibility than methods such as DMRG. In addition, its implementation is highly parallelizable allowing for further speed improvements.
In order to verify the method we have first studied the single impurity Anderson model. Comparing the ground state energy obtained using DMRG to the one obtained using the SGS ansatz, we have shown that an SGS state of rank R ∼ 6 can rival with high precision DMRG calculations with a large bond dimension m = 3000. Further increasing the ansatz rank up to R = 16 allows to improve the precision of the ground state energy estimate by an order of magnitude. To demonstrate the quality of the variational ground state, we have also carefully examined the spin-spin correlation functions and the impurity entanglement entropy, and found excellent agreement with analytical and previous numerical results.
To highlight the potential of the method, we have also studied the two-channel Kondo model. Studying again a spin-spin correlation function, we found signatures of the intermediate coupling fixed point in the form a sign change of the uniform part of the correlation function, and were able to confirm the expected exponential scaling of the Kondo screening length. We find a power law decay of the correlation function far from the impurity with exponent ν ≈ 1.58, which is close to the expected non-Fermi liquid behavior with exponent ν = 1.5 predicted by CFT calculations [31]. These results showcase the power of the SGS ansatz for the study of real-space properties of multichannel impurity models, a space of applications previously ill-covered by current standard numerical methods.
An open problem is whether there are computationally efficient ways of extracting the impurity Green's function (in real or imaginary time or frequency) within this class of ansatz states. Such a method would allow this approach to be integrated as impurity solver in embedding methods, such as the dynamical mean-field theory, which requires the impurity Green's function to achieve self-consistency. Furthermore, it appears possible to further reduce the scaling of the method by replacing the Gaussian covariance matrices by Gaus-sian fermionic matrix product states [43], which could in principle further reduce the scaling with the number of fermionic modes N . ACKNOWLEDGMENTS DMRG calculations were performed using the ITensor Library [30]. Numerical calculations of pfaffians were performed using the PFAPACK package [44]. The authors thank A.W.W. Ludwig for insightful discussions.
Appendix A: Numerical solution of Eq. (17) In this Appendix, we provide additional details on the numerical implementation of the variational method. There are two main approaches to solving Eq. (17) numerically. The first is to solve the equation using a generic differential equation solver, for example taking to first order with δs the step size. The drawback of this approach is that the accumulation of small numerical errors due to the finite step size will rapidly lead to covariance matrices which do not represent normalized and pure states. One must then frequently correct the normalization by decomposing the covariance matrix in the canonical form with R a real orthogonal matrix and rescaling the coefficients λ j to unity. A second approach is to rewrite Eq. (A1) as an orthogonal transformation [11] with the orthogonal matrix where we used (Γ µ ) 2 = −1 and neglected terms quadratic in δs. This second approach is numerically more expensive than the former due to the required matrix exponentiation. However, it preserves normalization up to machine precision. Our numerical experiments showed this second approach to perform better in some cases as it allows for taking larger time steps and one need not frequently correct the normalization of the covariance matrices using costly canonical form decomposition. In this Appendix, we derive the expression used to compute the order 2 Rényi entropy S 2 (|ψ , A) in Sec. IV C. We consider a bipartition of an SGS state |ψ in parts A and B. Defining the density matrices ρ α = |φ α φ α |, the reduced density matrix for subsystem A is since the states in the SGS ansatz are non-orthogonal and thus G α,β = φ α |φ β = 0. In order to evaluate Eq. (B1), we expand the density matrix in a fermionic coherent state basis [37]. The derivation is sketched below with the final result given by Eq. (B16). Using Eq. (B1), the Rényi entropy of order 2 of an SGS state is where we introduced the rank-4 tensor which obeys the relations T α,β γ,δ = (T β,α δ,γ ) * and T α,β γ,δ = T γ,δ α,β and where I α,β γ,δ = Tr A (Tr B (ρ α ρ β ) Tr B (ρ γ ρ δ )) . The final result, as a function of the covariance matrices of the Gaussian states, is given by inserting Eq. (B28) in the above equation.

Coherent state operator expansion
Setting first the required notation, we introduce for each fermionic mode the coherent states where η i , η i (i = 1, . . . N ) are Grassman variables obeying the usual algebra To lighten the notation we introduce the states |η = i |η i and η| = i η| i , as well as the shorthands η · η = i η i η i for products and d N η = N j=1 dη j dη j for differentials. In addition, we make the dependency on barred variables implicit such that F (η, η) → F (η).
Following the results of Ref. [37], any operator O can be represented by the integral where χ(O, η) is the characteristic function of the operator and the operator F (η) is Following the language of quantum optics D(η) = exp d † · η − η · d is the fermionic analogue of the bosonic displacement operator. In this work, the operator O in Eq. (B6) is the density matrix of a fermionic Gaussian state [23,35] with the canonical modes b j = i R i,j c i with R ∈ SO(2N ). The rotation matrix R and the eigenvalues λ j are defined through the normal form decomposition of the covariance matrix introduced in Eq. (A2). Using this decomposition and rewriting the displacement operator in this canonical basis, the trace in Eq. (B7) can be evaluated mode per mode leading to the characteristic function where the 2N × 2N matrix V N makes the rotation of the Grassman operators from a complex fermion basis to a real Majorana basis: The characteristic function Eq. (B10) allows to relates the Grassman variable representation used e.g. in Refs. [9,23] to expressions involving the density matrix through Eq. (B6).

Reduced density matrix
We now turn to deriving the necessary expressions for evaluating numerically Eq. (B1). To this end, we need to compute partial traces of the form Tr B (ρ α ρ β ), where Tr B is the partical trace over subsystem B which is composed of N B modes such that N = N A + N B .
Denoting states and Grassman variables associated with the subsystem by a subscript, the partial trace in the coherent state basis takes the form Using Eq. (B6) and noting that the F operator is separable such that F (η) = F (η A )F (η B ), we obtain where we introduced the kernel Using Eq. (B8) and evaluating the resulting Gaussian integral we obtain where we used for example −ψ|η = exp(−ψη). Inserting Eq. (B15) in Eq. (B13), we now isolate the integrals over subsystem B with the effective characteristic function for the subsystem A In order to take advantage of the structure of Eq. (B11), we formalize the basis change of Eq. (B11) by introducing the new Grassman variables θ i = j (V N ) i,j (η, η) j with the differentials transforming as dη j dη j = −i 2 dθ 2j−1 dθ 2j . Using(·) to indicate functions with arguments in the rotated basis, Eq. (B18) takes the form where Dθ B = 2N j=2N A +1 dθ j . In the case N B = N , i.e. when tracing over the whole system, we recover the result of Ref. [23].
Introducing the block structure of the covariance matrix where Γ A , Γ B are skewsymetric and Γ AB = −Γ T BA , one can evaluate the Gaussian integrals resulting iñ where we introduced the 4N J × 4N J skewsymmetric matrices (J ∈ {A, B}) (B21) Finally, one can transform back to the original basis to obtain Ξ(η A , η A , λ A , λ A ), by taking θ A = V N (η A , η A ) and φ A = V N (λ A , λ A ).