Entanglement entropy and computational complexity of the Anderson impurity model out of equilibrium I: quench dynamics

TL;DR

This paper demonstrates that in the Anderson impurity model, entanglement entropy grows logarithmically during quench dynamics, facilitating long-time simulations, with the behavior explained via a noninteracting chain model and basis optimization.

Contribution

It shows that logarithmic entanglement growth is generic in quenched impurity models when using energy-ordered bath orbitals, enabling extended real-time studies.

Findings

Logarithmic entropy growth observed in both interacting and noninteracting models.

Entropy growth behavior explained by a critical transition in the energy spectrum.

Basis choice significantly affects the entropy growth pattern.

Abstract

We study the growth of entanglement entropy in density matrix renormalization group calculations of the real-time quench dynamics of the Anderson impurity model. We find that with appropriate choice of basis, the entropy growth is logarithmic in both the interacting and noninteracting single-impurity models. The logarithmic entropy growth is understood from a noninteracting chain model as a critical behavior separating regimes of linear growth and saturation of entropy, corresponding respectively to an overlapping and gapped energy spectra of the set of bath states. We find that with an appropriate choices of basis (energy-ordered bath orbitals), logarithmic entropy growth is the generic behavior of quenched impurity models. A noninteracting calculation of a double-impurity Anderson model supports the conclusion in the multi-impurity case. The logarithmic growth of entanglement entropy enables studies of quench dynamics to very long times.

Figures

Equations

$$H=H_d+H_{\mathrm{bath}}+H_{\mathrm{mix}},$$

$$H_d=\sum _{\sigma }{ϵ}_d{d}_{\sigma }^{†}{d}_{\sigma }+U{d}_{\uparrow }^{†}{d}_{\uparrow }{d}_{\downarrow }^{†}{d}_{\downarrow },$$

$$H_{\mathrm{bath}}=\sum _{k\sigma }{ϵ}_k{c}_{k\sigma }^{†}{c}_{k\sigma },$$

$$H_{\mathrm{mix}}=\sum _{k\sigma }{V}_k{d}_{\sigma }^{†}{c}_{k\sigma }+h.c..$$

$$|{\Psi }_{t=0}⟩=|{\Psi }_0{⟩}_d\otimes |\mathrm{FS}{⟩}_{\mathrm{bath}},$$

$$\displaystyle\Delta_{\sigma}(t,t^{\prime})=\left\{\begin{array}[]{lr}\displaystyle-\frac{V^{2}}{E\tau}[H_{1}(E\tau)+iJ_{1}(E\tau)],&t\succ_{\mathcal{C}}t^{\prime},\\ \vspace{-1ex}&\\ \displaystyle-\frac{V^{2}}{E\tau}[H_{1}(E\tau)-iJ_{1}(E\tau)],&t\prec_{\mathcal{C}}t^{\prime},\end{array}\right.$$

$${\Delta }_{\sigma }(t,t^{\prime })$$

$$=-i\sum _{j=1}^N{V}_j^2\left[{\Theta }_{\mathcal{C}}(t,t^{\prime })-n_{j\sigma }^0\right]e^{-i{ϵ}_j(t-t^{\prime })}.$$

$$|\Psi (t)⟩=\sum _i$$

$$|\Psi (t)⟩=\mathcal{T}{e}^{-i{\int }_0^{t}d{t}^{\prime }{\hat{H}}_{\mathrm{mix}}(t^{\prime })}|{\Psi }_{t=0}⟩,$$

$${\hat{H}}_{\mathrm{mix}}(t)=e^{i{H}_{0}t}{H}_{\mathrm{mix}}{e}^{-i{H}_{0}t}$$

$$=\sum _{j\sigma }{V}_j{e}^{i(U{n}_{d\bar{\sigma }}+{ϵ}_d-{ϵ}_j)t}{d}_{\sigma }^{†}{c}_{j\sigma }+h.c.,$$

$${\tilde{H}}_{\mathrm{mix}}(t)$$

$$=\sum _{j\sigma }{\tilde{V}}_{j\sigma }(t)d_{\sigma }^{†}{c}_{j\sigma }+h.c.,$$

$${\tilde{V}}_{j\sigma }(t)=V_j{e}^{i(U{n}_{d\bar{\sigma }}+{ϵ}_d-{ϵ}_j)t}\mathrm{sinc}(\frac{U{n}_{d\bar{\sigma }}+{ϵ}_d-{ϵ}_j}{2}\Delta t).$$

$${\tilde{H}}_{\mathrm{mix}}(t)=\sum _{\sigma }(-1{)}^{n_{d\bar{\sigma }}}{d}_{\sigma }^{†}{\tilde{c}}_{\sigma }(t)+h.c.,$$

$${\tilde{c}}_{\sigma }(t)=\sum _j{\tilde{V}}_{j\sigma }(t)(-1{)}^{n_1+\cdots +n_{j-1}}{\tilde{c}}_{j\sigma },$$

$$c_{j\sigma }=(-1{)}^{n_d+n_1+\cdots +n_{j-1}}{\tilde{c}}_{j\sigma },$$

$${\tilde{c}}_{\sigma }(t)=\left[0,1\right]\prod _{j=1}^N\left[\begin{array}{cc}I& 0\\ {\tilde{V}}_{j\sigma }(t){\tilde{c}}_{j\sigma }& (-1{)}^{n_j}\end{array}\right]\left[\begin{array}{c}1\\ 0\end{array}\right],$$

$$|\Psi (t+\Delta t)⟩\approx e^{-i{\tilde{H}}_{\mathrm{mix}}\left(t+\frac{\Delta t}{2}\right)\Delta t}|\Psi (t)⟩$$

Full text

I Introduction

It has been a long-standing challenge to develop efficient real-time impurity solvers to study the dynamics of a system out of equilibrium. The Anderson impurity model Anderson (1961) (AIM), a single spin-degenerate orbital with an intra-orbital Hubbard interaction $U$ coupled to a bath of noninteracting orbitals, is of fundamental importance in its own right as a nontrivial but solvable Wiegmann and Tsvelick (1983); Tsvelick and Wiegmann (1983) interacting electron model and as an auxiliary problem for dynamical mean field theory Georges et al. (1996); Kotliar et al. (2006). The nonequilibrium properties of this model Aoki et al. (2014); Gramsch et al. (2013) provide an important laboratory for the development of real-time methods Hettler et al. (1998); Cohen et al. (2015); Balzer et al. (2015); Kreula et al. (2016).

The density matrix renormalization group (DMRG) Schollwöck (2005) is a powerful numerical method for solving low-dimensional electron problems including the Anderson impurity model. In DMRG, the wave function of the system is represented by a matrix product state (MPS). Every matrix in the MPS corresponds to a local degree of freedom. If these degrees of freedom are labeled by an index $n$, then the logarithm of the bond dimension $D_{n}$ between the $n$-th and $(n+1)$-th matrices of the MPS is greater than or of the order of the entanglement entropy $S_{n}$ between orbitals $1,2\ldots n$ and $n+1\ldots N$. The rate of growth with time of the entanglement entropy across the maximum entropy cut controls the maximum time achievable in a real-time DMRG calculation.

In this paper we study the growth of entanglement entropy in DMRG calculations, focusing primarily on the single-impurity Anderson model. Our work is motivated by the recent results of Wolf et al. Wolf et al. (2014), who found that a particular “star geometry” arrangement of bath orbitals provided a favorable growth of entanglement entropy. We consider a variety of bath Hamiltonians and find that the growth of entanglement entropy $S_{n}$ across the cut at $n$ is controlled by whether the spectra of MPS orbitals $1,2\ldots n$ overlaps that of the $n+1\ldots N$ orbitals. The star geometry with energy-ordered bath orbital arrangement in the MPS ensures no energy overlap at every bond, leading to a logarithmic growth of maximum entanglement entropy with time. This method of controlling entropy growth using energy separation is hypothesized to work for multi-impurity models as well; the hypothesis is supported by a calculation of a non-interacting double-impurity Anderson model.

The rest of the paper is organized as follows. Section II describes the single-impurity Anderson model (SIAM) and Sec. III describes our implementation of DMRG, which differs slightly from the standard implementation. In Secs. IV–VI, we show results obtained using our method for the SIAM, discuss the logarithmic growth of entropy and also present a few results for the quenched double-impurity Anderson model. Section VII is a conclusion and summary.

II Theory

We focus on the single-impurity Anderson model (SIAM) defined by the Hamiltonian

[TABLE]

Here $k$ labels the $\mathcal{N}\rightarrow\infty$ bath orbitals and $\sigma$ labels the spin. For simplicity we take the impurity-bath coupling amplitudes $V_{k}=V/\sqrt{\mathcal{N}}$ to be $k$-independent. We define the bath density of states $\mathrm{DOS}(\epsilon)=\frac{1}{\mathcal{N}}\sum_{k}\delta(\epsilon-\epsilon_{k})$ and consider a semicircle with a half band width $E$. The initial state that we consider is a product state

[TABLE]

where the Fermi-sea state $|\mathrm{FS}\rangle_{\mathrm{bath}}$ of the bath is initially half-filled. In the studies we present here, the $d$-orbital energies $\epsilon_{d}$ and $\epsilon_{d}+U$ are chosen to be symmetric about the Fermi level at [math]. This is not an essential requirement of our theory. The situation is depicted in Fig. 1. The formalism is easily generalized to a mixed state with an initial density matrix $\rho_{t=0}=(\rho_{0})_{d}\otimes(\rho_{0})_{\mathrm{bath}}$, but here we will focus on a pure initial state.

The numerical methods we use require a truncation of the bath to a finite number $N$ of bath orbitals. To choose the best truncation, we first calculate the hybridization function $\Delta_{\sigma}(t,t^{\prime})=-i\sum_{k}|V_{k}|^{2}\langle\mathcal{T_{C}}\,c_{k\sigma}(t)c_{k\sigma}^{\dagger}(t^{\prime})\rangle_{\mathrm{bath}}$. Here we work with the Keldysh-contour Green’s function $-i\langle\mathcal{T_{C}}\ldots\rangle_{\mathrm{bath}}$, where operators $c_{k\sigma}(t)$ and $c_{k\sigma}^{\dagger}(t^{\prime})$ are in the Heisenberg picture evolving via $H_{\mathrm{bath}}$ and the mean value $\langle\ldots\rangle_{\mathrm{bath}}$ is taken with respect to $|\mathrm{FS}\rangle_{\mathrm{bath}}$. With our choice of the semicircular DOS and constant $V_{k}$, the hybridization function can be evaluated at $\mathcal{N}\rightarrow\infty$ as

[TABLE]

where $\tau=t-t^{\prime}$, $H_{1}$ is the 1st-order Struve function and $J_{1}$ is the 1st-order Bessel function. The symbols $\succ_{\mathcal{C}}$ and $\prec_{\mathcal{C}}$ refer to Keldysh-contour ordering. Then we fit the hybridization function to that of a finite bath with only $N$ orbitals, i.e.,

[TABLE]

In fitting Eq. (9) with Eq. (10), all $2N$ real parameters $\epsilon_{j}$ and $V_{j}$ are varied to minimize the least-square error. The occupancies $n_{j\sigma}^{0}$ are chosen to be either [math] or $1$ to fit the $t\succ_{\mathcal{C}}t^{\prime}$ and $t\prec_{\mathrm{C}}t^{\prime}$ parts independently and to make the initial state of the finite bath a Slater determinant. This is possible even if the original bath was at nonzero temperature. Since our bath is particle-hole symmetric, we choose $N$ to be even to preserve this symmetry.

The number $N$ of bath orbitals controls the maximum time $t_{N}\lesssim 2\pi N/|E_{\mathrm{max}}-E_{\mathrm{min}}|=\pi N/E$ up to which the exact hybridization function is reproduced with good accuracy. For example, $N=40$ bath orbitals are enough to reach $Et\lesssim 100$ and $N=170$ orbitals can reach $Et\lesssim 500$. Adding more orbitals increases the maximal time that can be reached, but does not significantly improve the accuracy of the fit at shorter times.

III Method

We use DMRG/MPS methods to carry out the time evolution. We represent the SIAM wave function $|\Psi(t)\rangle$ as an entangled state between the impurity $d$ orbital and the bath, i.e.,

[TABLE]

where $i$ sums over the 4 impurity states $|0\rangle$, $\left|\uparrow\right>$, $\left|\downarrow\right>$ and $\left|\uparrow\downarrow\right>$. Every bath state $|\Psi_{i}(t)\rangle_{\mathrm{bath}}$ is a normalized matrix product state (MPS), with the $c_{i}(t)$’s being the normalizing coefficients. Eq. (11) is a Schmidt decomposition of $|\Psi(t)\rangle$ between the $d$ orbital and the bath if $|\Psi(t)\rangle$ is a simultaneous eigenstate of $N_{\uparrow}$ and $N_{\downarrow}$, the total numbers of spin-up and spin-down electrons. This representation differs from the conventional DMRG in that it removes the $d$ orbital from the MPS, enabling analysis of the entanglement among the bath orbitals.

We evolve the wave function $|\Psi(t)\rangle$ using the interaction picture of $H_{0}=H_{d}+H_{\mathrm{bath}}$. The wave function evolves according to

[TABLE]

where $\mathcal{T}$ is the time-ordering symbol and

[TABLE]

where $\bar{\sigma}$ is the opposite spin of $\sigma$. The main advantage of the interaction picture is that $\hat{H}_{\mathrm{mix}}(t)$ typically has a narrower spectral radius than $H_{0}$ (bath bandwidth $\sim E$ large compared with impurity level width $\sim V^{2}/E$).

We evaluate Eq. (12) by discretizing the time interval into time steps $\Delta t$. The operator at the time step centered on time $t$ is

[TABLE]

with the coupling amplitudes

[TABLE]

The errors of both the mid-point Hamiltonian $\hat{H}_{\mathrm{mix}}(t)$ and the time-averaged Hamiltonian $\tilde{H}_{\mathrm{mix}}(t)$ are $\mathcal{O}(\Delta t^{2})$. The latter choice is preferred if the bath bandwidth is large compared with the level width, because the very high and very low-energy bath orbitals are suppressed by the $\mathrm{sinc}$ function.

To apply the Hamiltonian $\tilde{H}_{\mathrm{mix}}(t)$ to the wave function $|\Psi(t)\rangle$ in Eq. (11), we work in the Jordan-Wigner transformed representation with the $d$ orbital being the first orbital ($d$ and $d^{\dagger}$ having no Jordan-Wigner signs). The Hamiltonian in Eq. (14) is rewritten as

[TABLE]

where the $\tilde{c}_{j\sigma}$ is the Jordan-Wigner transformed $c_{j\sigma}$. The two operators are related by

[TABLE]

so that the operators $\tilde{c}_{j\sigma}$ and $\tilde{c}_{j^{\prime}\sigma^{\prime}}$ for $j\neq j^{\prime}$ commute. We do the same Jordan-Wigner transform the two spins of the same orbital, so that $\tilde{c}_{j\uparrow}$ and $\tilde{c}_{j\downarrow}$ still anticommute. But this is easy to handle with a local $4\times 4$ matrix. The bath operator $\tilde{c}_{\sigma}(t)$ in Eq. (17) is then represented by a matrix-product operator (MPO)

[TABLE]

where the $j=1$ matrix is left-multiplied by $[0,1]$ to pick the second row, and the $j=N$ matrix is right-multiplied by $[1,0]^{T}$ to pick the first column. The MPO has a bond dimension of 2. We can similarly express $\tilde{c}_{\sigma}^{\dagger}(t)$ in terms of $\tilde{c}_{j\sigma}^{\dagger}$. The Hamiltonian $\tilde{H}_{\mathrm{mix}}(t)$ can then act on $|\Psi(t)\rangle$ following DMRG routines Schollwöck (2005).

The final evolution scheme is given by

[TABLE]

with the exponential factor Taylor expanded into a 4th-order polynomial of $\tilde{H}_{\mathrm{mix}}(t\!+\!\Delta t/2)$. The narrow spectral radius of $\tilde{H}_{\mathrm{mix}}\Delta t$ ensures good unitarity of the 4th-order truncation. Since the bath operators $\tilde{c}_{\sigma}$ and $\tilde{c}_{\sigma}^{\dagger}$ are long-range, we cannot locally exponentiate the Hamiltonian following the time-evolving block decimation (TEBD) Vidal (2003) method. We adjust the truncation error tolerance of the singular value decomposition (SVD) in DMRG routines according to the MPS norm so that higher-order terms do not take much time to calculate. We also parallelize the calculations of the 4 MPSs in Eq. (11) on 4 cores and use total numbers of spin-up and spin-down electrons as symmetries to speed up the calculation.

IV Results

In this section we show some results obtained for the interacting SIAM with $U/E=1$ using the method and other model parameters described in Secs. II and III. The impurity-bath coupling $V/E=0.1\sim 0.5$. This is the parameter range of interest. The impurity level width $\sim V^{2}/E$ remains smaller than the band width $\sim E$ while the Kondo temperature $T_{K}\approx 0.4Ve^{-\pi E^{2}/16V^{2}}$ (from Wang et al. (2008)) can change by orders of magnitudes.

Fig. 2a shows the charge relaxation dynamics, obtained by starting from an initially empty $d$ orbital $|\Psi_{0}\rangle_{d}=|0\rangle_{d}$ and a half-filled Fermi-sea state $|\mathrm{FS}\rangle_{\mathrm{bath}}$ for the bath. Our choice of particle-hole symmetric parameters ensures that $n_{d}=\langle n_{d\uparrow}\rangle+\langle n_{d\downarrow}\rangle$ always equilibrates to $1$ so long as the impurity-bath coupling $V$ is not big enough to form a bound state on the impurity. We see in agreement with previous work Eckstein et al. (2009); Wolf et al. (2014) that the charge equilibration proceeds relatively rapidly. The reciprocal of the time $t_{0.5}$ it takes to reach $n_{d}=0.5$ is plotted in Fig. 3a. At small $V/E\lesssim 0.1$, $t_{0.5}\sim V^{-2}$ is inversely proportional to the $d$-level width $\sim V^{2}/E$. For $V/E\gtrsim 0.15$, the rate $1/t_{0.5}$ crosses over to approximately linear in $V$ and the equilibration process in Fig. 2a becomes more oscillatory as we are approaching the formation of a bound state on the impurity. The variation of charge equilibration rates with $V$ can be seen in calculations performed for a noninteracting SIAM, suggesting that the charge relaxation physics is essentially due to hybridization. The Hubbard $U$ does not qualitatively change the behavior of the model.

Fig. 2b shows the spin relaxation dynamics obtained by starting from $\left|\uparrow\right>_{d}\otimes|\mathrm{FS}\rangle_{\mathrm{bath}}$, a fully spin-polarized $d$ orbital and the same half-filled Fermi-sea state $|\mathrm{FS}\rangle_{\mathrm{bath}}$ of the bath. The magnetization $m_{d}=\langle n_{d\uparrow}\rangle-\langle n_{d\downarrow}\rangle$ relaxes much more slowly than the charge, again in agreement with previous results Cohen et al. (2013); Rostami et al. (2016). The asymptotic behavior of $m_{d}$ v.s. $t$ shows approximately an exponential tail, with the relaxation rate $\Gamma_{\infty}\equiv d\ln m_{d}/dt|_{t\rightarrow\infty}$ plotted in Fig. 3b. $\Gamma_{\infty}$ is estimated by fitting $\ln m_{d}(t)$ v.s. $t$ to a straight line for $t_{\mathrm{max}}/2<t<t_{\mathrm{max}}$, where $t_{\mathrm{max}}$ is the maximum time reached in the simulation for the slope. The solid red line is a trend line. We also show as the dashed green line the analytical result — the Kondo temperature $T_{K}$ calculated using the formula in Wang et al. (2008) and interpreted as a relaxation rate.

The Kondo result has a similar magnitude and $V$ dependence to the calculated results. The numerical differences at large $V$ arise from relaxation processes associated with valence fluctuations not included in the Kondo limit, while the more pronounced differences at small $V$ are an intermediate asymptotics effect. For small $V$, even at the very long times ($Et\leq 600$) accessible to our method, the magnetization $m_{d}$ is still substantial, so the Kondo-limit expression, which gives the linear response relaxation for small magnetization ($m_{d}\rightarrow 0$), is not applicable. Evidently, the nonlinear response (relaxation of a finite $m_{d}$) is stronger than the linear response. Developing a theory of the relaxation in the small $V$ and intermediate $m_{d}$ regime is an interesting open question. For intermediate $V/E\simeq 0.25$, the theoretical result is within a factor of $2$ of the numerical one with the differences likely arising from the convention used for the Kondo temperature $T_{K}$.

V Logarithmic growth of entropy

A remarkable feature of the simulations reported here is the long time scales that can be reached; these time scales are necessary to reveal, for example, the magnetization decay. As we show in this section, this is possible because the maximum entanglement entropy of the 4 bath MPSs in Eq. (11) grows only logarithmically during the simulation, which means the long times are not exponentially hard to reach, but are of only polynomial time complexity.

V.1 Entanglement entropy growth in SIAM

In this section, we compare the maximum entanglement entropy of the interacting SIAM ($U/E=1$) with a noninteracting SIAM ($U=0$) with $\epsilon_{d}=0$ at the Fermi level. Both models start from the same initial condition $|0\rangle_{d}\otimes|\mathrm{FS}\rangle_{\mathrm{bath}}$ with an empty $d$ orbital and a half-filled bath in Fig. 1. Results of the entanglement entropy are shown in Fig. 4. The entropy growth starting from a spin-polarized impurity $\left|\uparrow\right>_{d}\otimes|\mathrm{FS}\rangle_{\mathrm{bath}}$ is also logarithmic but takes smaller values.

The curves in Figs. 4a and 4b are obtained in slightly different ways. Fig. 4a shows results obtained for a noninteracting simulation of $N=1000$ bath orbitals all coupled to one empty $d$ orbital at the Fermi level. We then plot the entanglement entropy between the 500 bath orbitals below the Fermi level with the other 501 orbitals up to $Et=1000$ and the data shows a logarithmic growth of entanglement entropy at all values of the impurity-bath coupling $V$. At long times, the slopes of the curves are the same; only the offset and the transient entropy growth depend on $V$.

Fig. 4b shows the interacting model ($U/E=1$, $\epsilon_{d}=-U/2$). At each time, the hybridization function was fitted with the minimal number of bath orbitals needed to obtain a fit with a root-mean-square error (RMSE) of $3\times 10^{-4}$. Notice that the hybridization fit is independent of $U$ and can be done before simulation. We then plot the maximum entanglement entropy $S_{\mathrm{max}}$ seen on all bonds of the 4 bath MPSs encountered during the simulation from $t=0$ to the maximum simulation time $t_{\mathrm{max}}$ allowed by the hybridization fit v.s. $t_{\mathrm{max}}$. $S_{\mathrm{max}}$ may be encountered before $t_{\mathrm{max}}$ due to the finite bath effect. So Fig. 4b takes into account the possibility of using the finite bath effect to limit entropy growth. But still the logarithmic growth of entropy and the independence of the steady-state slope of $S-\log t$ on the impurity-bath coupling $V$ is the same as the noninteracting SIAM simulation in Fig. 4a. These two properties mean $S\leq c\ln t$, and therefore the bond dimension $D\sim e^{S}\leq t^{c}$, which means the interacting SIAM can be simulated in polynomial time $\mathcal{O}(D^{3})=\mathcal{O}(t^{3c})$ of $t$.

V.2 Analysis of entropy growth

To understand the logarithmic growth of entropy, we consider a noninteracting chain model, as is shown in Fig. 5. In this model, the impurity is coupled to two semi-infinite chains. We choose a constant hopping amplitude between the bath sites in each chain. By adjusting the on-site energy difference of the two chains, we can vary the densities of states as shown in the lower panel of Fig. 5, obtaining either overlapping, gapped or just touching spectra. Our computation of the entanglement entropy $S_{\mathrm{occ}}$ across the impurity site shows that we have linear growth, logarithmic growth, and saturation, respectively. In the numerical test we did, chain $a$ was initially empty and chain $b$ was initially full. But the conclusion is found to hold for randomized initial occupancies, too.

The logarithmic growth of $S_{\mathrm{occ}}$ is seen at a critical $E=V_{a}+V_{b}$, at which the density of states (DOS) of the two semi-infinite chains touch at only one energy point. When $E>V_{a}+V_{b}$, the system is gapped and entropy growth saturates. This can be explained by the lack of energy eigenstates that are extended in both regions $a$ and $b$, which then means that particles (or holes) that are originally in $a$ cannot go into $b$ and vice versa beyond a penetration depth determined by the gap, which then puts an upper bound on the entanglement entropy between $a$ and $b$. This energy barrier works for a general initial occupancy. Starting from any product state, so long as the semi-chains $a$ and $b$ are gapped, the entropy must saturate.

When $E<V_{a}+V_{b}$, there is a finite overlap of the DOS of the two semi-infinite chains and we see a linear growth of entanglement entropy in Fig. 5. In rare cases this does not happen. For example, for a uniform chain $V_{a}=V_{b}=V_{c}$ and $E=\epsilon_{d}=0$, the entropy growth is logarithmic rather than linear. But this behavior depends on the initial occupancy. If the occupied sites are randomized, or if the model parameters are slightly modified to deviate from a uniform chain, the expected behavior of a linear growth of the entanglement entropy is seen between $a$ and $b$. The energy criterion guarantees that particles do not enter the forbidden regions of a noninteracting bath. But once the energy barrier is not at work, it is difficult in general, though not impossible, to organize the migrated particles into a low entanglement entropy state to make the MPS matrices small.

The logarithmic growth of entropy in Fig. 4b can be understood as the result of arranging the bath orbitals in the MPS in energy order, so that at any bond of the MPS, the left and right parts of the bath degrees of freedom always have touching energy spectra. This argument applies to an interacting model, too, because the bath is still noninteracting, and the Hubbard $U$ only reduces the chance for the impurity – the only bridge via which the bath orbitals can indirectly hop to one another – to be doubly occupied, thus reducing its bridging efficiency. The bath entanglement entropy of an interacting SIAM is therefore upper bounded by that of a noninteracting SIAM from this picture.

V.3 Bath in chain geometry

So far we have been working in the star geometry of the bath. Bath orbitals do not hop to each other directly. They only do so via the impurity. The diagonalization of bath orbitals in energy space leads to a logarithmic growth of entanglement entropy, according to the energy criterion in the previous section. In this section, we would like to emphasize again that the energy criterion is a sufficient but not necessary condition for the entropy to grow slowly. The example to give here is the evolution of the quenched SIAM in the chain geometry of the bath. The impurity is the head of the chain, which is directly connected to only one bath orbital, which in turn is connected to another bath orbital, and so on so forth. One can go from the star geometry to the chain geometry via Lanczos tridiagonalization starting from the impurity orbital, and from the chain back to the star by diagonalizing the bath. More details of the two geometries can be found in Wolf et al. (2014).

Starting from the initial state $|0\rangle_{d}\otimes|\mathrm{FS}\rangle_{\mathrm{bath}}$ with $|\mathrm{FS}\rangle_{\mathrm{bath}}$ being the same filled Fermi-sea state as in Fig. 1 transformed to the chain geometry, the maximum entropy on the chain (the entanglement entropy between the left and right parts of the chain at the maximum entropy cut) is still found to grow logarithmically. Fig. 6 shows the result of a noninteracting calculation. The initial occupancies on the chain are spatially uniform. Every site has an occupancy of 0.5 per spin except the empty impurity. The entanglement entropy $S_{n}$ between sites $1,2,\ldots n$ and $n+1,\ldots N$ on the chain are then plotted in Fig. 6 as a function of $n$ at equal intervals of time. On top of the logarithmic background of $S_{n}$ of the equilibrium state $|\mathrm{FS}\rangle_{\mathrm{bath}}$, an entropy peak propagates like a soliton from the impurity down the chain at a speed $\sim E$. The maximum entanglement entropy (height of the peak) therefore increases with time logarithmically, even though there is no separation of energy spectrum on the chain, i.e., partition of the bath into different regions with different energies like in the star geometry.

Starting from an inverted half-filled Fermi-sea state with the highest energies initially occupied, the same logarithmic growth of entropy is seen due to particle-hole symmetry. But starting from a state with random initial occupancies in the star geometry, the entanglement profile in the chain geometry becomes prohibitively high ($\mathrm{max}(S_{n})\propto N$) even at $t=0$. Also, a linear growth of entropy is seen starting from a product state in the chain geometry with [math]-$1$ random initial occupancies (see Fig. 7), while in the chain geometry, the entropy growth is still logarithmic starting from a random product state with [math]-$1$ occupancies. These results demonstrate that the logarithmic entropy growth in Fig. 6 is not guaranteed by the MPS basis of the Hamiltonian, but is due to the initial filled Fermi-sea state. For such a special initial state, the star geometry does not have a big advantage over the chain geometry, as they both give a logarithmic growth of maximum entanglement entropy. The benefit of the star geometry is its good behavior for more general initial states.

It is important to point out, as is shown in Fig. 7, that the star geometry alone does not guarantee a logarithmic entropy growth. The order of the bath orbitals in the MPS matters. The initial occupancies affect the transient growth of entropy, while the asymptotic entropy growth is determined by the ordering of the bath orbital energies.

VI Double-impurity model

In this section, we show that the logarithmic growth of entanglement entropy is not limited to the single-impurity Anderson model by doing a noninteracting simulation of a double-impurity Anderson model. The most general noninteracting double-impurity Anderson model can be pictorially represented in Fig. 8. Fig. 8a is in the basis in which the 2 impurity orbitals and all bath orbitals are diagonal, which is the double-impurity version of the star geometry. Fig. 8b shows the double-impurity version of the chain geometry by Lanczos tridiagonalizing the star geometry in Fig. 8a starting from the two impurities. One can also tridiagonalize the bath orbitals above and below the Fermi level separately (Fig. 8c) to obtain the double-impurity generalization of the two semi-infinite chain model in Fig. 5. Since the left and right semi-chains have touching energy spectra, the logarithmic growth of entropy is expected as a critical behavior between linear growth and saturation of entropy, as is discussed in Sec. V.2.

Figure 9 shows a sample result. We chose a half-filled bath with a semicircle DOS the same as Fig. 1, and put two $d$ orbitals at $\pm 0.2E$ ($E$ is the half band width) with $d-d$ hopping $0.15E$ to mimic typical crystal field splitting. The two $d$ orbitals are equally coupled to all bath orbitals. In the basis in which the two $d$ orbitals are diagonalized, their orbital energies are $\pm 0.25E$ and the original $d-d$ hopping makes the two $d$ orbitals now couple to the bath differently, which is more realistic. Then we plot the entanglement entropy $S_{\mathrm{occ}}$ between the initially occupied bath orbitals and the rest of the system.

The double-impurity model has a richer dynamics than SIAM. Since both impurities are initially empty, the one below the Fermi level leaks a hole into the bath, leading to a short-term entropy peak. The steady-state growth of $S_{\mathrm{occ}}$ is still logarithmic, but the slope of $S_{\mathrm{occ}}$ v.s. $\log t$ is not constant. This is because the two $d$ orbitals are not at the Fermi level (one is above and one is below). Their distances in energy to the Fermi level $|\epsilon_{d1,2}|$ relative to the impurity-bath coupling $V$ determine the slope, which approaches a maximum for the case of a $d$-orbital at the Fermi level ($\epsilon_{d}=0$) as $V$ gets large.

The logarithmic growth of entropy again shows that the quenched multi-impurity model is not exponentially hard in DMRG simulations, but is of only polynomial-time complexity. Whether the conclusion still holds for interacting models needs further investigation in DMRG, especially for those multi-impurity models with non-density-density (spin flipping and pair-hopping) terms, whose entanglement entropies need not be bounded by the corresponding noninteracting models.

VII Conclusion

We have studied the growth of entanglement entropy in quenched Anderson impurity models. It is found that the growth of entropy is determined by the representation of the bath orbitals in the matrix product state (MPS). The Hubbard $U$ on the impurity orbital does not change the qualitative behavior of the steady-state growth of entanglement entropy of the bath MPS. The crucial feature controlling the growth of entanglement entropy is the overlap in energy of the density of states of the two parts of the maximum entanglement partition. In the star geometry of energy-ordered bath orbitals, the touching-spectra condition is satisfied at every bond, so the maximum bond dimension is power law in $t$. The power is upper bounded by the case of a half-filled $d$-orbital at the Fermi level and does not grow with the impurity-bath coupling, which allows a simulation of the long-time dynamics of the quenched impurity models in polynomial time. The conclusion is likely to generalize to multi-impurity models.

The growth of entanglement entropy of an interacting quantum system and the associated computational cost has been studied previously Prosen and Žnidarič (2007); Pižorn et al. (2014) in terms of the integrability of the quantum model. Our study looks at the problem from a different perspective. We focus on a special class of quantum models — the impurity models — and think of the growth of entanglement entropy among the bath orbitals. Because of the sparsity of interactions in the model, the entropy growth in the noninteracting bath is controlled by the energy partitioning of the bath and the localization of bath electrons to the energies they belong to. Since the new criterion of energy-partitioning the bath is not related in obvious ways to the integrability of the whole model (bath $\!+\!$ impurity), hopefully this new view of entropy growth of complexity can help us find new polynomial-time solvable models, parameter ranges, and/or special initial conditions that are not covered by the previously established integrability criterion.

Acknowledgments: We thank Dr. Dante Kennes for helpful discussions. This research is supported by the Department of Energy under grant DE-SC0012375.