Entanglement entropy and computational complexity of the Anderson impurity model out of equilibrium II: driven dynamics

TL;DR

This paper investigates how entanglement entropy and computational complexity evolve in driven Anderson impurity models, proposing orbital reordering techniques to optimize matrix product state simulations of periodically driven quantum systems.

Contribution

It introduces a novel approach to reduce entanglement growth in driven impurity models by reordering bath orbitals based on Floquet quasi-energy, enhancing simulation efficiency.

Findings

Reordering bath orbitals reduces exponential growth of entanglement entropy.

Growth of entanglement is linked to bath orbital ordering and driving period.

Optimized orbital ordering improves matrix product state simulation performance.

Abstract

We study the growth of entanglement entropy and bond dimension with time in density matrix renormalization group simulations of the periodically driven single-impurity Anderson model. The growth of entanglement entropy is found to be related to the ordering of the bath orbitals in the matrix product states of the bath and to the relation of the driving period $T$ to the convergence radius of the Floquet-Magnus expansion. Reordering the bath orbitals in the matrix product state by their Floquet quasi-energy is found to reduce the exponential growth rate of the computation time at intermediate driving periods, suggesting new ways to optimize matrix product state calculations of driven systems.

Figures

Equations

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

$$H_d(t)=\sum _{\sigma }{ϵ}_d(t)n_{d\sigma }+U(t)(n_{d\uparrow }-\frac{1}{2})(n_{d\downarrow }-\frac{1}{2}),$$

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

$$H_{\mathrm{mix}}(t)=\sum _{k\sigma }{V}_k(t)d_{\sigma }^{†}{c}_{k\sigma }+h.c.,$$

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

$$=\sum _{k\sigma }{V}_k(t){\hat{d}}_{\sigma }^{†}(t){\hat{c}}_{k\sigma }(t)+h.c.,$$

$${\hat{c}}_{k\sigma }(t)=c_{k\sigma }{e}^{-i{\int }_0^t{ϵ}_k(t^{\prime })d{t}^{\prime }},$$

$${\hat{d}}_{\sigma }(t)=d_{\sigma }{e}^{-i{\int }_0^t[{ϵ}_d(t^{\prime })+U(t^{\prime })(n_{d\bar{\sigma }}-\frac{1}{2})]d{t}^{\prime }},$$

$$|\Psi (t)⟩=\sum _i{c}_i(t{)}_{}|i{⟩}_d\otimes |{\Psi }_i(t){⟩}_{\mathrm{bath}},$$

$$|\Psi (t+\Delta t)⟩$$

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

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

$${\tilde{V}}_{k\sigma }(t)$$

$$\times \mathrm{sinc}\left(\frac{{ϵ}_d(t)+U(t)(n_{d\bar{\sigma }}-1/2)-{ϵ}_k(t)}{2}\Delta t\right).$$

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

$$\displaystyle\epsilon_{d}(t)=\left\{\begin{array}[]{ll}-|\epsilon_{d}|,&\displaystyle 0<t<\frac{T}{2},\\ +|\epsilon_{d}|,&\displaystyle\frac{T}{2}<t<T,\end{array}\right.$$

$$e^{-i{H}_{F}T}\equiv e^{-i{H}_{+}T/2}{e}^{-i{H}_{-}T/2},$$

$$H_F=\bar{H}+i|{ϵ}_d|\tan \left(\frac{T}{4}{\mathrm{ad}}_{\bar{H}}\right)n_d+\mathcal{O}(|{ϵ}_d{|}^2),$$

$$S_{N/2}\sim c\ln T\ln t,$$

$$e^{-i{H}_{F}T}=\mathcal{T}{e}^{-i{\int }_0^{T}dt[H_0+ϵH_1(t)]}.$$

$$\mathcal{T}{e}^{-i{\int }_0^{T}dt[H_0+ϵH_1(t)]}=e^{-i{H}_{0}T}$$

$$-iϵ{\int }_0^{T}dt{e}^{-i{H}_0(T-t)}{H}_1(t)e^{-i{H}_{0}t}+\mathcal{O}({ϵ}^2).$$

$$H_F=H_0+ϵ\delta H_F^{(1)}+\mathcal{O}({ϵ}^2).$$

$$e^{-i{H}_{F}T}=e^{-i(H_0+ϵ\delta H_F^{(1)})T}+\mathcal{O}({ϵ}^2)=e^{-i{H}_{0}T}$$

$$-iϵ{\int }_0^{T}dt{e}^{-i{H}_0(T-t)}\delta H_F^{(1)}{e}^{-i{H}_{0}t}+\mathcal{O}({ϵ}^2).$$

$${\int }_0^{T}dt{e}^{i{H}_{0}t}{H}_1(t)e^{-i{H}_{0}t}={\int }_0^{T}dt{e}^{i{H}_{0}t}\delta H_F^{(1)}{e}^{-i{H}_{0}t},$$

$$e^{i{H}_{0}t}{H}_1(t)e^{-i{H}_{0}t}=\sum _{n=0}^{\infty }{\frac{(it{)}^n}{n!}}_{}{\mathrm{ad}}_{H_0}^n[H_1(t)],$$

$$H_1(t)=H_1\mathrm{sgn}\left(t-\frac{T}{2}\right),0\le t<T,$$

$$\sum _{n=0}^{\infty }\frac{(iT{)}^n}{(n+1)!}\left(1-\frac{1}{2^n}\right){\mathrm{ad}}_{H_0}^n(H_1)$$

$$=\sum _{n=0}^{\infty }{\frac{(iT{)}^n}{(n+1)!}}_{}{\mathrm{ad}}_{H_0}^n(\delta H_F^{(1)}),$$

$$\frac{(e^{i\frac{T}{2}{\mathrm{ad}}_{H_0}}-1{)}^2}{iT{\mathrm{ad}}_{H_0}}{H}_1=\frac{e^{iT{\mathrm{ad}}_{H_0}}-1}{iT{\mathrm{ad}}_{H_0}}\delta H_F^{(1)}.$$

$$\delta H_F^{(1)}=i\tan \left(\frac{T}{4}{\mathrm{ad}}_{H_0}\right)H_1.$$

$$⟨m|\delta H_F^{(1)}|n⟩=i⟨m|H_1|n⟩\tan \left(\frac{E_m-E_n}{4}{}_{}T\right).$$

Full text

I Introduction

The control of strongly correlated electron systems via laser-induced oscillating electric fields is becoming an active area of research Paoli et al. (2006); Banerjee and Shah (2018); Basak et al. (2018); Boschker et al. (2015); Shaw et al. (2014), making it important to study the dynamics of driven strongly correlated models, i.e., evolution of the system with time-dependent model parameters. The Anderson impurity model (AIM) Anderson (1961), a correlated orbital coupled to a noninteracting bath, is a model of fundamental interest and is important as an auxiliary model in the dynamical mean-field approach Georges et al. (1996); Kotliar et al. (2006); Gramsch et al. (2013) to correlated electron physics. The development of efficient real-time impurity solvers Hettler et al. (1998); Cohen et al. (2015); Balzer et al. (2015); Blümer (2007) to study the nonequilibrium properties of strongly correlated models such as the Anderson impurity model is a key challenge in this field of research.

We analyze the use of density matrix renormalization group (DMRG) White (1992); Schollwöck (2005); Kin-Lic Chan and Sharma (2010) as a nonequilibrium impurity solver Cazalilla and Marston (2002); White and Feiguin (2004); Wolf et al. (2014) for the periodically driven AIM, where the system’s wave function is represented by a matrix product state (MPS). The key issue in DMRG calculations is the growth of entanglement entropy of the MPS with simulation time. In our previous study He and Millis (2017) of “quench” physics (i.e., the evolution following an instantaneous change of interaction and/or hybridization parameters from one set of constant values to another), we found that different arrangements of bath orbitals could dramatically affect the growth of entanglement entropy, and that a particular arrangement (the “star geometry” Wolf et al. (2014), associated with a proper energy ordering of bath orbitals) led to a very slow (logarithmic) growth of entanglement entropy, enabling simulations of the long-time behavior at a computational cost that grew only polynomially with the simulation time.

In this paper, we study the periodically driven single-impurity Anderson model (SIAM) and find that driven systems are in general more expensive to simulate than quenched systems. There is a critical driving period $T_{c}$, namely $2\pi$ over the band width of the bath density of states, such that if the driving period $T\!<\!T_{c}$, the driven system is as easy to simulate as a quenched SIAM, while if $T\!>\!T_{c}$, the simulation becomes exponentially hard. For driving periods $T\!>\!T_{c}$, we find that an ordering of the bath orbitals exists, which we call quasi-energy ordering, such that the asymptotic entanglement entropy growth is slow (logarithmic in time, for the noninteracting model, and with a small linear coefficient for the interacting model). However, the initial transient growth of entropy for this bath ordering can be very rapid before the asymptotic limit is reached, which limits the maximum simulation times reachable in practice.

As in our previous work He and Millis (2017), we use the 4-MPS scheme developed previously to simulate the real-time dynamics of the driven SIAM. The wave function is represented by a Schmidt decomposition between 4 states of the impurity orbital ($|0\rangle,\left|\uparrow\right>,\left|\downarrow\right>,\left|\uparrow\downarrow\right>$) and 4 bath MPSs correspondingly. The Hamiltonian $H$ of the quenched model is time-independent, while the Hamiltonian $H(t)$ of the driven model oscillates in time periodically.

The rest of the paper is organized as follows. Section II describes the driven SIAM we solve and generalization of our 4-MPS method He and Millis (2017) to the driven model. In Sec. III, we present results obtained for the noninteracting SIAM to obtain the asymptotics of entropy growth behaviors in various parameter regimes to scan the complexity diagram. In Sec. IV, we simulate the interacting SIAM using our method. We will discuss entropy growth and show some physical results of the model. Section V is a conclusion and summary.

II Theory and method

We consider a single-impurity Anderson model (SIAM) with periodically oscillating model parameters. The most general form of the time-dependent model Hamiltonian $H(t)$ is given by

[TABLE]

[TABLE]

where $n_{d\sigma}=d_{\sigma}^{\dagger}d_{\sigma}$ and $\sigma=\,\uparrow,\downarrow$ is the spin label. We go to the interaction picture of $H_{0}(t)\equiv H_{d}(t)+H_{\mathrm{bath}}(t)$. The $H_{\mathrm{mix}}(t)$ part in the interaction picture becomes

[TABLE]

where $U_{0}(t,0)=\mathcal{T}e^{-i\int_{0}^{t}H_{0}(t^{\prime})dt^{\prime}}$ is the time-ordered unitary evolution from [math] to $t$ and $U_{0}(0,t)=[U_{0}(t,0)]^{\dagger}$. Since $H_{0}(t)$ does not couple the $d$ orbital to the bath, each bath orbital evolves independently in the interaction picture as given by

[TABLE]

with $\bar{\sigma}$ denoting the opposite spin of $\sigma$. Notice that $\hat{n}_{d\bar{\sigma}}(t)=n_{d\bar{\sigma}}$ does not evolve in the interaction picture of $H_{0}(t)$ and that $n_{d\bar{\sigma}}$ commutes with $d_{\sigma}$, which together lead to Eq. (6b).

As in the 4-MPS scheme developed in our previous work He and Millis (2017), the wave function $|\Psi(t)\rangle$ is represented by

[TABLE]

where $i$ sums over the $4$ impurity states $|0\rangle$, $\left|\uparrow\right>$, $\left|\downarrow\right>$, and $\left|\uparrow\downarrow\right>$, and every $|\Psi_{i}(t)\rangle_{\mathrm{bath}}$ is a matrix product state. The wave function $|\Psi_{i}(t)\rangle$ is evolved according to

[TABLE]

with the exponential Taylor expanded to 4th order of $\Delta t$ to ensure good unitarity. The time-averaged Hamiltonian $\tilde{H}_{\mathrm{mix}}(t)$ in a time step $\Delta t$ is now given by

[TABLE]

with effective hopping amplitude given by

[TABLE]

The Hamiltonian $\tilde{H}_{\mathrm{mix}}(t)$ is represented by a matrix product operator (MPO) with bond dimension 2 to act upon the wave function in 4 MPSs He and Millis (2017).

Up to now everything has been general for the single-impurity Anderson model. For concreteness, we consider the evolution starting from a product state

[TABLE]

where $|\mathrm{FS}\rangle_{\mathrm{bath}}$ is a half-filled Fermi-sea state of the bath with a semicircle density of states (DOS) as shown in Fig. 1. The $\mathcal{N}\rightarrow\infty$ bath orbitals have fixed energies $\epsilon_{k}(t)=\epsilon_{k}$ and fixed equal hopping amplitudes $V_{k}(t)=V/\sqrt{\mathcal{N}}$ to the impurity $d$ orbital. The Hubbard $U$ on the $d$ orbital is also fixed. The only time-dependent quantity is the $d$-orbital energy

[TABLE]

which oscillates in a square wave every half driving period $T/2$. A piecewise constant Hamiltonian is numerically easier to handle because Eq. (8) can be made exact by choosing $\Delta t$ such that $T/2$ is a multiple of $\Delta t$. We will provide physical results which show the local quantities on the $d$ orbital and complexity results which show the growth of entanglement entropy of the bath MPSs at the maximum entropy bond, which is often the one closest to the Fermi level.

III Noninteracting Results

We first do some noninteracting calculations of the driven SIAM using a standard Slater-determinant-based method. In this section, the Hubbard $U=0$ and the initial state is $|0\rangle_{d}\otimes|\mathrm{FS}\rangle_{\mathrm{bath}}$, an empty $d$-orbital and a half-filled Fermi sea in Fig. 1. The impurity-bath coupling $V/E=0.25$. Bath size $N=1000$ orbitals. The time-evolution of a Slater-determinant state by a noninteracting Hamiltonian is numerically cheap and is not limited by the growth of entanglement entropy. The 4-MPS method will be applied in the next section to an interacting SIAM with Hubbard $U>0$.

III.1 Energy-ordered bath

We use the entanglement entropy $S_{\mathrm{occ}}$ between the $N/2$ bath orbitals below the Fermi level and the rest of the system to estimate the maximum entanglement entropy encountered in an MPS-based simulation when the bath orbitals are energy-ordered. We find that for long driving periods $T>T_{c}=\pi/E$, since our bath DOS is gapless (see Fig. 1), an arbitrarily small amplitude $|\epsilon_{d}|$ changes the logarithmic growth of entropy to linear. This is shown in Fig. 2, where we did a simulation with $N=1000$ bath orbitals and driving period $ET=10$. The entanglement entropy $S_{\mathrm{occ}}$ between the $500$ bath orbitals below the Fermi level and the rest of the system is plotted in Fig. 2 over time. The quenched model exhibits a logarithmic growth of entanglement entropy and the entanglement entropy growth in the periodically driven model is linear in simulation time.

The linear growth of entanglement entropy in a driven SIAM can be intuitively understood in an entropy pumping picture. The up and down motion of the $d$ orbital acts as an elevator that transports some electrons from the occupied bath orbitals to the unoccupied bath orbitals (and holes in the opposite direction). So if the entanglement entropy $S_{\mathrm{occ}}$ increases by a constant $(\Delta S_{\mathrm{occ}})_{T}$ in every period, then the linear growth rate of $S_{\mathrm{occ}}$ would be $(\Delta S_{\mathrm{occ}})_{T}/T$.

The critical driving period $T_{c}=\pi/E$, or $2\pi$ over the band width, separates the logarithmic growth ($T<T_{c}$) and linear growth ($T>T_{c}$) of $S_{\mathrm{occ}}$. Fig. 3 plots the maximum growth rate of entropy $(\Delta S_{\mathrm{occ}})_{T}/ET$ v.s. the period $T$ as an envelope of the growth rate v.s. $T$ curves at fixed driving amplitudes $|\epsilon_{d}|$. Each of these curves is tangent to the envelope at some points and they all intersect with zero at the same critical driving period $T_{c}$. For periods $T<T_{c}$, the linear growth of entropy cannot be maintained. To understand this critical period, one needs to consider the Floquet Hamiltonian $H_{F}$ of the driven system defined by

[TABLE]

where $H_{\pm}$ corresponds to $\epsilon_{d}=\pm|\epsilon_{d}|$ respectively. The time-independent Floquet Hamiltonian $H_{F}$ reproduces the unitary evolution of the time-dependent system $H(t)$ over full periods. It turns out $ET=\pi$ is the convergence radius of the Floquet-Magnus expansion of $H_{F}$ in terms of $H_{+}$ and $H_{-}$. Within the convergence radius and for small driving amplitude $|\epsilon_{d}|$, we have

[TABLE]

where $\bar{H}=(H_{+}+H_{-})/2$ is the SIAM Hamiltonian with $\epsilon_{d}=0$, $\mathrm{ad}_{\bar{H}}=[\bar{H},\cdot\,]$ is the adjoint representation of $\bar{H}$, and $\tan(\cdot)$ is defined via its Taylor expansion. A detailed derivation of Eq. (16) is in Appendix A. The convergence radius is $ET=\pi$. We expect this to hold also for the interacting SIAM, because in the thermodynamic limit $N\rightarrow\infty$, the spectral radius $\|\bar{H}\|$ is mainly determined by the band width of the bath DOS (unless a bound state is formed on the impurity). As a result, $\|\bar{H}\|\approx E$ is equal to the half band width $E$ of the bath. Since $\tan(\cdot)$ is singular at $\pi/2$, the series expansion of Eq. (16) fails to converge if $\|\mathrm{ad}_{\bar{H}}\|_{\,}T/4=\|\bar{H}\|_{\,}T/2\approx ET/2>\pi/2$, i.e., for long periods $ET>\pi$.

Once the perturbative expansion does not converge, surprising new physics emerges. Outside the convergence radius ($T\!>\!T_{c}$), numerics shows that the bath orbital energies in $H_{F}$ are aliased to $[-\pi/T,\pi/T]\subset[-E,E]$, breaking the original ordering of the bath orbitals. This situation is shown in Fig. 4. The inter-bath-orbital hopping amplitudes remain very small. The ascending order of bath orbital energies is violated because of energy aliasing, which gives overlap of energy spectra between the occupied and unoccupied bath orbitals and thus a linear growth of entanglement entropy.

III.2 Quasi-energy-ordered bath

What happens then if one reorders the bath orbitals in the MPS in ascending order of quasi-energy rather than energy when the driving period $T\!>\!T_{c}$? We again use the Slater-determinant-based noninteracting simulation to make an estimation. The initial state $|0\rangle_{d}\otimes|\mathrm{FS}\rangle_{\mathrm{bath}}$ in the star geometry remains a product state (an MPS with bond dimension $=1$). We use the entanglement entropy $S_{N/2}$ between the $N/2$ bath orbitals with negative quasi-energies (within $[-\pi/T,0)$) and the rest of the system to estimate the maximum entanglement entropy encountered in an MPS-based simulation when the bath orbitals are quasi-energy-ordered. The $S_{N/2}$ defined here becomes equivalent to the entanglement entropy $S_{\mathrm{occ}}$ used in the previous subsection if $ET<\pi$, when the bath orbitals are energy-ordered.

We redo the same simulation as in Fig. 2 using $N=1000$ bath orbitals ordered by their quasi-energies of $ET=10$. The entropies of the energy-ordered simulation in Fig. 2 (blue and red lines) are compared with the new results (green and purple lines) in Fig. 5 and the time $t$ is put on log scale. It is found that the growth of $S_{N/2}$ is logarithmic for both the quenched and driven models. This is because the Floquet Hamiltonian $H_{F}$ is now energy-ordered, as opposed to the aliased situation in Fig. 4. But the driven model is still harder to simulate than the quenched model, because the slope of the $S_{N/2}$ v.s. $\ln t$ curve is greater for the driven model.

For the quenched model, the steady-state slope of $S_{N/2}$ v.s. $\ln t$ is unchanged when the bath orbitals are quasi-energy-ordered. Only the steady-state intercept is shifted up by a constant $\Delta S_{N/2}$, which is found to be approximately proportional to $\ln(T/T_{c})$ (see Fig. 6a). This is the price to pay for not ordering the quenched bath by energy, which is better than a randomly shuffled bath (see Fig. 7 in He and Millis (2017)), whose entropy $S_{N/2}$ would grow linearly with time $t$.

For the driven model, the driving period $T$ changes the slope of the $S_{N/2}$ v.s. $\ln t$ curve. Fig. 6b shows how the slope increases from that of the quenched model ($T\rightarrow 0$ at fixed $|\epsilon_{d}|$ is equivalent to quench) to unboundedly large values proportional to $\ln T$. This indicates that the leading-order term in the entropy $S_{N/2}$ is

[TABLE]

where $c$ depends on $|\epsilon_{d}|$ but is found to be bounded (see Fig. 7). At very large $|\epsilon_{d}|\gtrsim E$, the coefficient $c$ goes down, which is likely to come from the bound state formed on the impurity.

Eq. (17) means that the bond dimension in an MPS-based simulation using the quasi-energy-ordered bath arrangement is $D\sim e^{S_{N/2}}\sim t^{\,c\ln T}$. The time complexity of the singular value decomposition (SVD) step is then $\mathcal{O}(D^{3})=\mathcal{O}(t^{3c\ln T})$.

Since the power of $t$ for the quasi-energy-ordered method is unbounded for long driving periods $T$, the complexity is still beyond polynomial time. Another drawback of quasi-energy ordering is delocalization of maximum entanglement entropy throughout the MPS, while in energy-ordered MPSs, the maximum entanglement entropy tends to concentrate near the Fermi level. This gives the quasi-energy-ordered method a prefactor of the bath size $N$.

IV Interacting results

In the previous section, we have been estimating what would happen in an MPS-based simulation using a Slater-determinant-based code for the noninteracting SIAM. Now let us do some real MPS-based simulations of the interacting SIAM using the 4-MPS method developed in Sec. II. We choose a fixed Hubbard $U/E=1$ and the impurity-bath coupling $V/E=0.25$ is the same as in Sec. III. We use $N=30$ bath orbitals to fit the hybridization function of the continuum bath DOS in Fig. 1 with good accuracy up to $Et\leq 75$ following He and Millis (2017). The SVD truncation error tolerance was $10^{-5}$. Noninteracting $d$-occupancies are reproduced with $2\sim 3$ decimal places as a benchmark.

IV.1 Physical results

First let us show how the physical results of the interacting SIAM differ from the noninteracting SIAM. The results of short periods $ET<\pi$ are not significantly different from the quenched SIAM with no oscillation of $d$-orbital energy. So we plot both Figs. 8 & 9 in the long-period regime $ET\!>\!\pi$. Both the energy-ordered and quasi-energy-ordered algorithms as discussed in Sec. III give the same physical results.

The Hubbard $U$ suppresses the double occupancy of the $d$-orbital for both the quenched and driven SIAMs. In Fig. 8, the dashed grey lines indicate the level of double occupancy in a noninteracting SIAM (estimated from the $n_{d}^{2}/4$ of the quenched $n_{d}$). The interacting double occupancy is appreciably lower than $n_{d}^{2}/4$ when the driving amplitude $|\epsilon_{d}|/E=0.1$ is small. For period $ET=10$, both $n_{d}$ (red line in Fig. 8a) and the double occupancy $D$ (purple line) oscillate in sinusoidal waves, even though the driving signal $\epsilon_{d}(t)$ is a square wave. When the period increases to $ET=20$, the wave forms approach a relaxed oscillation (Fig. 8b). The overshoots in every period disappear in a noninteracting simulation ($U=0$, not plotted), which produces simple monotonic decays to the square wave levels.

When the driving amplitude $|\epsilon_{d}|$ is increased, the wave form of $n_{d}$ distorts, and the relaxation to steady-state oscillation slows down, as is shown in Fig. 9. Also, there is an increase of the average double occupancy $\overline{D}$. At $|\epsilon_{d}|/E=0.8$, the double occupancy $D$ in its oscillation steady state is above $n_{d}^{2}/4$ almost the entire period. A possible explanation might be that the oscillating $d$-orbital energy is like a phonon mode that induces an effective intra-$d$-orbital attraction, which becomes greater than $U$ when the oscillation amplitude $|\epsilon_{d}|$ is big enough ($|\epsilon_{d}|/E\gtrsim 0.6$, at which $D\approx 1/4$). Whether this attractive interaction can lead to superconductivity is interesting for further studies.

IV.2 Complexity results

Obtaining results in Fig. 9a at medium to large driving amplitudes was not easy, as they were in the $ET\!>\!\pi$ regime. The linear growth of maximum entanglement entropy makes the maximum bond dimensions in the MPSs increase exponentially with the number of periods simulated. We used some extrapolation techniques to estimate the steady-state quantities in Fig. 9b, especially for $|\epsilon_{d}|/E=0.8$ where the relaxation is slow. Here, in this section we mainly check whether this linear growth of entropy (i.e., exponential difficulty) can be helped by reordering the bath orbitals in the MPSs inquasi-energy order.

We find that even though the entropy growth in the noninteracting SIAM changes from linear to logarithmic by quasi-energy ordering the bath orbitals, as is shown in Sec. IIIB, the entropy growth for the interacting SIAM is slightly faster than logarithmic (see red curve in Fig. 10a). We increase the number of bath orbitals to $N=40$ to reach $Et=100$, and then make a comparison of the energy-ordered and quasi-energy-ordered simulations in Fig. 10 under $|\epsilon_{d}|/E=0.1$, $ET=6$. As is shown in Fig. 10b, the quasi-energy-ordered 4-MPS simulation is slower than the energy-ordered simulation in the short run. The short-term growth of entropy, e.g. in the first few periods, is faster if the energies of the bath orbitals are not ordered. In the long run, the quasi-energy ordering is more favorable. The entropy growth only slightly curves up in the $S_{\mathrm{max}}$ v.s. $\ln t$ plot. The long-term growth rate of entropy and $\ln t_{\mathrm{CPU}}$ v.s. $t$ in Fig. 10b are clearly reduced. The hardness in the $ET>\pi$ regime is beyond polynomial time using either method, but is significantly reduced by the quasi-energy ordering method.

Figure 11 shows the crossing time of the maximum entanglement entropies $S_{\mathrm{max}}$ of the energy ordered and quasi-energy ordered simulations. In a wide range of driving periods the crossing time $t_{\mathrm{cross}}$ of the entropies in Fig. 10 exists and is minimum at intermediate driving periods $T$ at which the linear growth rate of entropy $S_{\mathrm{max}}$ of the energy ordered method is fastest. After the entropies cross, the quasi-energy-ordered method still needs to overcome two more short-term drawbacks: a) its maximum entanglement entropy being more widespread over the MPS bonds than the energy-ordered method, and b) the bigger entropy at short times less than one period, before the actual CPU-times cross.

V Conclusion

We have generalized our previously developed 4-MPS method to time-dependent Hamiltonians to study the driven SIAM and have analyzed the computational complexity in the short-period $ET<\pi$ and long-period $ET>\pi$ regimes for both the noninteracting ($U=0$) and interacting ($U>0$) models. The model behavior in the $ET<\pi$ regime is not significantly different from the quenched model. This is the regime in which the Floquet-Magnus expansion converges. Both the interacting and noninteracting models are as easy to simulate as the quenched models (polynomial time). In the $ET>\pi$ regime, the entropy grows linearly in the energy-ordered algorithm, which is therefore exponentially hard to reach long times (many periods). Using quasi-energy ordering reduces the entropy growth of the noninteracting model from linear to logarithmic with a coefficient of the logarithm that grows unboundedly with the driving period $T$ (proportional to $\ln T$). For the interacting model, quasi-energy ordering significantly reduces the linear growth rate of entropy and thus the exponential hardness grows with time more slowly in the long run.

The main challenge of future research is the lack of analytic results for the Floquet Hamiltonian of the interacting model beyond the convergence radius of the Magnus expansion ($ET\!>\!\pi$). Therefore the “good basis” of the MPSs has to be guessed. In our work, we made a guess to use quasi-energy ordering based on the noninteracting model. To make better guesses that further reduce entanglement entropy growth, one would need to know the Floquet Hamiltonian of the driven interacting model in greater detail.

Acknowledgments: This research is supported by the Department of Energy under grant DE-SC0012375.

Appendix A Floquet Hamiltonian

The Floquet Hamiltonian $H_{F}$ of a periodically driven system $H(t)=H_{0}+\epsilon H_{1}(t)$ is given by

[TABLE]

For small amplitudes we have $\epsilon\rightarrow 0$. We can take derivative with respect to $\epsilon$ at $\epsilon=0$ to obtain

[TABLE]

We define an expansion

[TABLE]

Then we have

[TABLE]

Comparing Eqs. (19) and (21), we have from the first-order terms of $\epsilon$ that

[TABLE]

where we have multiplied on both sides by $e^{iH_{0}T}$ from the left. Then we use the nested commutator expansion

[TABLE]

where $\mathrm{ad}_{H_{0}}\equiv[H_{0},\cdot]$ is the adjoint representation of $H_{0}$, and $\mathrm{ad}_{H_{0}}^{n}[H_{1}(t)]=[H_{0},\mathrm{ad}_{H_{0}}^{n-1}[H_{1}(t)]]$ is the $n$-fold nested commutator of $H_{0}$ with $H_{1}(t)$. Using this formula on both sides of Eq. (22), and from the square wave model

[TABLE]

we have

[TABLE]

or in functional form

[TABLE]

All functions of $\mathrm{ad}_{H_{0}}$ are defined using their power series in Eq. (25). We now apply the multiplicative inverse of the power series of $\mathrm{ad}_{H_{0}}$ on the right-hand side to both sides and after some algebra obtain

[TABLE]

The formal solution in Eq. (27) can be evaluated in the eigenbasis of $H_{0}$ as

[TABLE]

where $|m\rangle$ and $|n\rangle$ are eigenstates of $H_{0}$ with eigen-energies $E_{m}$ and $E_{n}$. In case $H_{1}=|\epsilon_{d}|n_{d}$ with $n_{d}\equiv\sum_{\sigma}d_{\sigma}^{\dagger}d_{\sigma}$, Eq. (16) in the main text is derived. Since no assumption is made on $H_{0}$ except it is time independent, Eqs.(27), (28) & (16) hold for both the interacting and the noninteracting SIAMs.