Breaking the entanglement barrier: Tensor network simulation of quantum transport

TL;DR

This paper introduces a basis transformation approach that significantly improves tensor network simulations of quantum transport by exploiting the natural frequency structure of particle scattering, enabling longer and larger simulations.

Contribution

The authors demonstrate that choosing the frequency basis for quantum transport problems reduces entanglement growth, greatly enhancing tensor network simulation capabilities for non-equilibrium systems.

Findings

Simulation efficiency can be exponentially improved using the frequency basis.

Tensor network methods can be extended to larger spatial and temporal scales.

The approach broadens the class of non-equilibrium quantum systems accessible to tensor network simulation.

Abstract

The recognition that large classes of quantum many-body systems have limited entanglement in the ground and low-lying excited states led to dramatic advances in their numerical simulation via so-called tensor networks. However, global dynamics elevates many particles into excited states, and can lead to macroscopic entanglement and the failure of tensor networks. Here, we show that for quantum transport -- one of the most important cases of this failure -- the fundamental issue is the canonical basis in which the scenario is cast: When particles flow through an interface, they scatter, generating a "bit" of entanglement between spatial regions with each event. The frequency basis naturally captures that -- in the long-time limit and in the absence of inelastic scattering -- particles tend to flow from a state with one frequency to a state of identical frequency. Recognizing this natural structure yields a striking -- potentially exponential in some cases -- increase in simulation efficiency, greatly extending the attainable spatial- and time-scales, and broadening the scope of tensor network simulation to hitherto inaccessible classes of non-equilibrium many-body problems.

Equations

$$\sqrt{T(0)}|0_{\mathcal{L}}{1}_{\mathcal{R}}⟩+\sqrt{1-T(0)}|1_{\mathcal{L}}{0}_{\mathcal{R}}⟩$$

$$S\approx H[T(0)]\frac{\mu t}{2\pi },$$

$$H=H_{\mathcal{S}}+H_{\mathcal{I}}+H_{\mathcal{L}}+H_{\mathcal{R}},$$

$$H_{\mathcal{I}}=\sum _{i\in \mathcal{S},k\in \mathcal{L}\mathcal{R}}\hslash v_{ik}\left(c_i^{†}{a}_k+a_k^{†}{c}_i\right)$$

$$H_{\mathcal{L}(\mathcal{R})}=\sum _{k\in \mathcal{L}(\mathcal{R})}\hslash {\omega }_k{a}_k^{†}{a}_k,$$

$$H_{\mathcal{L}}=\hslash {\omega }_0\sum _{j\in \mathcal{L}}\left(b_j^{†}{b}_{j+1}+b_{j+1}^{†}{b}_j\right)+\hslash {\mu }_{\mathcal{L}}\sum _{j\in \mathcal{L}}{b}_j^{†}{b}_j$$

$$H_{\mathcal{S}}=\hslash {\omega }_{\mathcal{S}}\sum _{\sigma }{n}_{\sigma }+\hslash U{n}_{\uparrow }{n}_{\downarrow }$$

$$H_{\mathcal{I}}=\hslash v\sum _{\sigma ,l=1_{\mathcal{L}},1_{\mathcal{R}}}\left(c_{\sigma }^{†}{b}_{l\sigma }+b_{l\sigma }^{†}{c}_{\sigma }\right),$$

$${\omega }_k=2{\omega }_0\cos \left(k\pi /\left(N+1\right)\right)+{\mu }_{\mathcal{L}},$$

$$I(t)=-⟨d{n}_{\mathcal{L}}/dt⟩=-2\sum _{\sigma ,k\in \mathcal{L}}{v}_k\Im ⟨a_k^{†}{c}_{\sigma }⟩,$$

$${\mathcal{C}}_{mn}=\mathrm{tr}[d_m^{†}{d}_n\rho ],$$

$$H=\sum _{m,n\in \mathcal{LS}\mathcal{R}}{\bar{H}}_{nm}{d}_m^{†}{d}_n$$

$$\dot{\mathcal{C}}=-ı[\bar{H},\mathcal{C}]/\hslash .$$

$$I_{\mathrm{ss}}^{\mathrm{ex}}$$

$$\approx$$

$$T(\omega )={\Gamma }_L{\Gamma }_R|G(\omega ){|}^2,$$

$$H=\prod _{{\omega }_k<0}{W}^k\cdot W^{\mathcal{S}}\cdot \prod _{{\omega }_k>0}{W}^k,$$

$$W^{{\omega }_k<0}=\left(\begin{array}{cccc}\text{mathds}1& v_{ik}{a}_k& v_{ik}{a}_k^{†}& {\omega }_k{a}_k^{†}{a}_k\\ & \text{mathds}1& & \\ & & \text{mathds}1& \\ & & & \text{mathds}1\end{array}\right),$$

$$W^{\mathcal{S}}=\left(\begin{array}{cccc}\text{mathds}1& -c_i& c_i^{†}& {\omega }_{\mathcal{S}}{n}_i\\ & & & -c_i^{†}\\ & & & c_i\\ & & & \text{mathds}1\end{array}\right),$$

$$W^{{\omega }_k>0}=\left(\begin{array}{cccc}\text{mathds}1& & & {\omega }_k{a}_k^{†}{a}_k\\ & \text{mathds}1& & v_{ik}{a}_k^{†}\\ & & \text{mathds}1& v_{ik}{a}_k\\ & & & \text{mathds}1\end{array}\right).$$

$$S_{\mathrm{eff}}=\ln \sqrt[3]{\frac{1}{L-1}\sum _n{e}^{3S_n}},$$

$$\begin{array}{r}\text{includegraphics}\end{array},$$

$$\begin{array}{r}\text{includegraphics}\end{array}.$$

$$U_n^{[0]}(dt)C_n$$

$$U_n^{[1]}(dt)A_n^c$$

$$U_{n,n+1}^{[2]}(dt)A_n^c{A}_{n+1}^c$$

$$U(dt)\approx T^{[1]}(dt/2)T^{[1]\ast }(dt/2)\approx T^{[2]}(dt/2)T^{[2]\ast }(dt/2).$$

$$T^{[1]}(y)=U_L^{[1]}(y)U_{L-1}^{[0]}(-y)U_{L-1}^{[1]}(y)\dots U_2^{[0]}(-y)U_2^{[1]}(y)U_1^{[0]}(-y)U_1^{[1]}(y),$$

$$T^{[2]}(y)=U_{L-1,L}^{[2]}(y)U_{L-1}^{[1]}(-y)U_{L-2,L-1}^{[2]}(y)\dots U_3^{[1]}(-y)U_{2,3}^{[2]}(y)U_2^{[1]}(-y)U_{1,2}^{[2]}(y),$$

$$\begin{array}{r}\text{includegraphics}\end{array}.$$

Full text

Tensor networks enable the systematic search for ground states of certain many-body Hamiltonians, as well as numerical time evolution, provided that there is a limited amount of entanglement present Orus (2019); Ran et al. (2020); Orús (2014); Eisert (2013); Schollwöck (2011); Verstraete et al. (2008). Quantum quenches – when a parameter of the Hamiltonian is suddenly changed – can, though, generate highly-entangled states, seen both experimentally Kaufman et al. (2016) and theoretically Alba and Calabrese (2017); Liu and Suh (2014); Kim and Huse (2013); Schachenmayer et al. (2013); Schuch et al. (2008); Calabrese and Cardy (2005). The large amount of entanglement creates a challenge for tensor network simulation and the efficient representation of the underlying quantum state Schuch et al. (2008); Eisert (2013). There are many approximate approaches under development to truncate further the description of the state and maintain control over the size of the tensor network White et al. (2018); Leviatan et al. (2017); Surace et al. (2019), but these rely on additional assumptions, such as the thermalizing nature of the dynamics. We will here develop a controllable approach to break the entanglement barrier for an important class of problems in transport.

Quantum transport through an impurity $\mathcal{S}$ is a paradigmatic example of a “pathological” quench. A global bias $\mu$ drives particles through an interface where they scatter, see Fig. 1(a). For particles around the Fermi level ($\omega=0$), each scattering event gives rise to an entangled electron-hole pair Beenakker (2006); Klich and Levitov (2009)

[TABLE]

across the left ($\mathcal{L}$) and right ($\mathcal{R}$) reservoir regions. The left component of the state represents a particle transmitted from $\mathcal{L}$ to $\mathcal{R}$ with transmission probability $T(0)$ and the right component the reflected particle (the phase is unimportant here). Given that the attempt frequency is $\mu/2\pi$ Levitov and Lesovik (1993), the entanglement entropy $S$ increases as

[TABLE]

where $H[T(0)]$ is the binary entropy of the transmission probability and $t$ the time Beenakker (2006); Klich and Levitov (2009).

This linear growth of spatial entanglement – and its “light cone” spread Chien et al. (2014) (see the heat map in Fig. 1(a) – is due to the linear increase in the number of entangled electron-hole pairs, as expressed by Eqs. (1,2).

This growth results in the failure of one-dimensional tensor networks – matrix product states (MPS) – beyond a “hard wall”: The required matrix product dimension $D$ increases exponentially with the timescale. Figure 2 shows this spectacular failure for the non-interacting Anderson impurity model (see the caption for its definition). This intrinsic, physically-based limitation restricts MPS to short timescales and small/moderately-sized lattices Cazalilla and Marston (2002); Zwolak and Vidal (2004); Gobert et al. (2005); Schneider and Schmitteckert (2006); Schmitteckert and Schneider (2006); Al-Hassanieh et al. (2006); Dias da Silva et al. (2008); Heidrich-Meisner et al. (2009); Branschädel et al. (2010); Chien et al. (2013); Gruss et al. (2018), or linear response via an equilibrium correlation function Bohr et al. (2006); Bohr and Schmitteckert (2007). Simulating time-dependent problems (artificial gauge fields, Floquet states, etc.), more complex many-body regions, or long relaxation times, requires a new approach.

The linear growth in entanglement entropy – or its consequence, the uncontrolled growth in matrix product dimension – is deceptive: The paradigmatic impurity model will, in the long time limit, have particles go from a state of frequency $\omega$ (on the left) to a state of the same frequency $\omega$ on the right, albeit with some characteristic spread. This entails that if one instead works with the single-particle eigenbasis of $\mathcal{L}$ and $\mathcal{R}$, ordered on a lattice as shown in Fig. 1(b), the entanglement should be limited (in higher dimensions, momentum conservation can play this role). In fact, recently it was shown, in the context of dynamical mean-field theory, that for real-time single-particle correlation functions (and equilibrium) the so-called star geometry, where the energy basis for the bath is used, suppresses entanglement from logarithmic into a localized structure with smaller overall magnitude Wolf et al. (2014). For quenches in the Anderson impurity model, it was shown that energy basis ordering naturally delineates the bad (linear entanglement growth) and good (limited entanglement) scenarios He and Millis (2017).

Unlike these cases, we address simulating the bad scenario and show that it can be transformed to a scenario with logarithmic growth, and thus intermediate between the bad and good. To do so, we use a mixed energy and spatial basis, reflecting the entanglement structure in Eq. (1) and incorporating the energy basis in two separate spatial regions. Figure 1 shows the steps leading to this mixed basis (diagonalizing the separate $\mathcal{L}$ and $\mathcal{R}$ spatial regions and then ordering them). Entanglement in this mixed basis is localized to the bias window and mostly between pairs of (iso- or nearly iso-energetic) sites, see the heat map in Fig. 1(b). The strength of the couplings to the impurity, as well as many-body interactions and inelastic scattering, determine the spread of entanglement. At the same time, the low dimensionality of $\mathcal{S}$ and the scattering nature of the states limits the amount of entanglement between the impurity and the reservoirs. We will comment on alternative structural representations later.

We note here that various approaches can perform computations with matrix product states, such as the density matrix renormalization group (DMRG) White (1992), the time-dependent variational principle (TDVP) Haegeman et al. (2016), or Krylov-based methods Zaletel et al. (2015). Unlike the schemes based on the Trotter decomposition of the Hamiltonian into local gates, they allow treating any Hamiltonian represented as a matrix product operator (MPO). Thus, we use DMRG to find the initial ground state in the preferred basis and TDVP for the subsequent time evolution. Since we work with the MPO of the Hamiltonian, its dimension is important since the formal scaling for time evolution in, e.g., 1D is $\mathcal{O}(D^{3}Md+D^{2}M^{2}d^{2})$, where $M$ is the MPO bond dimension and $d$ is the local Hilbert space dimension. When the reservoirs are non-interacting, the MPO of the mixed basis has a small, fixed $M$ for both the initial state and time evolution regardless of bias Sup .

Since there is overhead associated with the presence of long-range interactions, we also work under a guiding principle that both the Hamiltonian MPO and the state MPS should have limited $D$. The mixed basis, in contrast to the spatial basis (exponentially large MPS $D$) and the global single particle basis (extensive MPO $D$ when interactions are present), respects this principle in addition to capturing the natural structure of impurity transport Sup . Optimality questions aside, it permits an advantageous extension to open systems Wójtowicz et al. (2019) where a bias is maintained by external contacts to $\mathcal{L}$ and $\mathcal{R}$ eigenstates separately Gruss et al. (2016, 2017); Elenewski et al. (2017), as well as a suitable structure for fine–graining the reservoirs Zwolak (2008).

Figure 2(a) shows the result of employing this mixed-basis MPS. An inhomogeneous and modest $D\leq 128$ already gives excellent results out to a time extensive in the system size (a time equal to the reservoir length divided by the Fermi velocity, $2\omega_{0}$, at which the current “front” hits the open boundary and travels backwards toward the impurity Chien et al. (2014, 2012)). This mixed basis captures the linear growth in spatial entanglement entropy, Fig. 2(b), and allows for very large systems, Fig. 2(c).

As a consequence of naturally representing the entanglement structure, the majority of the lattice in the mixed basis has little entanglement across bipartite cuts with correlations predominantly between modes in the bias window (see Fig. 1(b)). Thus, the computational speedup is not just a consequence of a reduced $D$, but also an inhomogeneous $D$. As a point of comparison, the mixed-basis simulations in Fig. 2(a) took only 15 hours, whereas the spatial-basis simulation with the same $D=128$ took 44 hours, both on the same single core computer. While implementation choices affect this comparison, the empirical scaling follows from Fig. 2: To bring the spatial-basis simulation out to $t=128\,\omega_{0}$ requires five more doublings of $D$ just to move the breaking point (forgetting about overall error). The dominant $D^{3}Md$ contribution to the computational cost then indicates an approximate computational time of $(2^{5})^{3}\cdot 44$ hours, or 165 years.

Even for a large bias, the mixed-basis MPS still out performs the spatial basis, a direct consequence of the more local nature of entanglement in energy. Figure 3(a) shows the conductance and current versus time traces for $\mu=0$ to $4\,\omega_{0}$ (after which the $\mathcal{L}$ and $\mathcal{R}$ bands go out of alignment and there is no steady state current and entanglement saturates). In all cases, the mixed basis yields accurate results. Figure 3(b) shows that the effective entanglement entropy Sup grows logarithmically opposed to linear. The advantages of the mixed basis also extend to interacting reservoir models, where Fig. 3(c,d) show the current versus time and effective entanglement entropy for the same $\mathcal{L}\mathcal{S}\mathcal{R}$ system but with interactions in $\mathcal{L}\mathcal{R}$. The Hamiltonian’s MPO dimension grows with the system size in this case, and in some parameter regimes we expect a spatial basis may be more suitable (such as in a localized regime), but the mixed basis enables stable time evolution. Further work will be necessary to assess the efficiency gains across different parameter regimes, since many-body interactions in $\mathcal{L}$ and $\mathcal{R}$ modify both the MPO and entanglement structure of the problem.

The above can be straightforwardly applied to higher dimensional non-interacting reservoirs (since only their energy/momentum bases matter) and to interacting spinless fermions, including larger dimensional systems $\mathcal{S}$. For many-body systems of typical interest, though, one has to have spin, which requires simulating multiple channels. The Anderson impurity problem Anderson (1961) with electron-electron interaction $Un_{\uparrow}n_{\downarrow}$ at the impurity (where $n_{\uparrow(\downarrow)}$ is $\mathcal{S}$’s spin up (down) number operator), and its extension to larger $\mathcal{S}$ Bauer et al. (2013); Iqbal et al. (2013), is the paradigmatic example.

In the non-interacting limit the two spin channels are fully independent. In the presence of interactions, $U\neq 0$, the interchannel entanglement originates from the many-body contact at $\mathcal{S}$, which is in addition to the intrachannel entanglement around the bias window Sup . This suggests an X-shape MPS (i.e., a tree tensor network Shi et al. (2006); Murg et al. (2010)) depicted in Fig. 4(a) as a natural ansatz (a structure also supported by results of other recent works Bauernfeind et al. (2017)). Tree tensor networks, similarly to a one-dimensional MPS, possess a normal form. As such, the TDVP integration scheme of Ref. Haegeman et al. (2016) directly extends to such tensor network geometries Sup .

Figure 4(b) shows the X-lattice simulated in both the spatial and mixed bases. Just as with the non-interacting case, the spatial basis abruptly fails and increasing $D$ only gives logarithmic increase in the achievable simulation time. The mixed basis, though, enables the simulation to go out to a time extensive in the reservoir size. When $D$ is too small, it will lose accuracy, but it does not abruptly fail. This is reflected in the limited growth in entanglement, Fig. 4(c), which behaves similarly to the non-interacting case.

Figure 4(d) shows the conductance diagram of the paradigmatic Anderson impurity problem. For negative $U$, $\mathcal{S}$ is approximately half occupied in each channel, giving a $U/2$ contribution to $\mathcal{S}$’s on-site energy in the other channel. This results in a single high conductance state when the level energy is pushed into the bias window at $\omega_{\mathcal{S}}\approx-U/2$. As $\omega_{\mathcal{S}}$ becomes negative, the conductance peak bifurcates into two particle-hole dual, correlated high-conductance states. For one, there is a correlated state between one channel being occupied and current flowing in the other channel, giving $\omega_{\mathcal{S}}\approx-U$ to effectively push the current-carrying channel state into the bias window. The other is a correlated state between one channel being empty and current flowing in the other channel. This occurs at $\omega_{\mathcal{S}}\approx-\omega_{0}$ instead of $\omega_{\mathcal{S}}=0$ due to residual many-body correlations increasing the energy (a residual also present in the $\omega_{\mathcal{S}}\approx-U$ state). The accurate calculation of the whole conductance diagram enables the identification of these features.

Finally, we comment on the computational speedup. The spatial basis requires exponentially large $D$ in the total simulation time $T$, already alluded to above: Each doubling of $D$ increments the breaking point by $\Delta t$ (independent of the value of $D$), giving $\Delta t\ln D\approx T$ or $D\approx e^{T/\Delta t}$. For the mixed-basis simulation of the single channel model, we examine the error versus $D$ for several simultaneous multiples of the reservoir size and time Sup . The simultaneously, e.g., doubling of size and time is the method by which one achieves the long-time limit. The error decay suggests that doubling of the simulation time (and size) requires increasing $D$ to $\alpha D$, where $\alpha$ is bounded, to keep the overall error fixed. Thus, the computational cost is brought from $e^{3T/\Delta t}$ to $T^{p}$, where $p\approx 3$ Sup . We note, however, that for fermions with spin, the X-lattice configuration requires also evolution between the two channels. The entropy across this bond is the same in the spatial and mixed basis, and can itself increase linearly in time. For the range of parameters here, it is still quite small. In principle, this will dominate the scaling for long times. However, one will still get an exponential improvement in the prefactor of this contribution, since that prefactor depends in the intrachannel entanglement and thus is suppressed when going to the mixed basis. Other structures besides the X-lattice may improve this further.

The above general considerations demonstrate that difficult computational problems can be broached so long as the natural entanglement structure is recognized – here, by changing the canonical basis. This enables the accurate simulation of quantum transport that underlies many applications, from quantum dot platforms for computing to molecular and nanoscale electronic devices, and fundamental studies with cold-atom emulators. The long times achievable will be conducive to simulating transport through systems undergoing time-dependent driving to generate artificial gauge fields or Floquet states. Combining the approach here with other recent methods Dorda et al. (2015); Schwarz et al. (2018); Fugger et al. (2018) will push further the limits of simulation, as will developing algorithms to locally optimize the canonical basis Krumnow et al. (2016). As such, our results open new avenues to study the behavior and simulation of non-equilibrium many-body systems, from fermionic impurities to bosonic baths to the inherent structure of tensor networks.

I Supplementary Material

Quantum transport through an impurity or interface is typically approached via the Hamiltonian Jauho et al. (1994); Caroli et al. (1971)

[TABLE]

where $H_{\mathcal{S}}$ is the many-body Hamiltonian of the impurity region – the system $\mathcal{S}$ – which may include electron-electron interactions, electron-phonon/vibrational coupling, etc. The remaining terms are the coupling of the system and reservoirs, and the isolated reservoir Hamiltonians. Reflecting the non-interacting nature of the Fermi sea and recognizing that the partitioning into $\mathcal{S}$ and $\mathcal{L}\mathcal{R}$ can be done so that the relevant, spatially localized interaction region is in $\mathcal{S}$, it is standard Jauho et al. (1994); Caroli et al. (1971) to take these other Hamiltonians as quadratic forms

[TABLE]

and

[TABLE]

where $c_{i}^{\dagger}$ ($c_{i}$) and $a_{k}^{\dagger}$ ($a_{k}$) are the creation (annihilation) operators in $\mathcal{S}$ and $\mathcal{L}\mathcal{R}$, respectively, and $v_{ik}$ is the coupling for modes $i\in\mathcal{S}$ and $k\in\mathcal{L}\mathcal{R}$. Spin (when present) is implicit in the labels. This “impurity” Hamiltonian is the same general structure as that addressed with the Numerical Renormalization Group (of course, one can modify Eq. (S2) to have some other $\mathcal{S}$ operator but retaining the linearity in $\mathcal{L}\mathcal{R}$ operators). There, a logarithmic discretization of the energy basis in the reservoir(s) gives a finite number of modes with a more fine structure at low energy. After a transformation of this Hamiltonian with a finite number of degrees of freedom to a one-dimensional lattice (a Wilson chain), an iterative diagonalization process yields the low energy states Wilson (1975); Bulla et al. (2008).

Matrix product state simulations require a (quasi-) 1D lattice. Prior simulations thus considered either explicitly a lattice in one spatial dimension Cazalilla and Marston (2002); Zwolak and Vidal (2004); Schneider and Schmitteckert (2006); Schmitteckert and Schneider (2006); Al-Hassanieh et al. (2006); Dias da Silva et al. (2008); Heidrich-Meisner et al. (2009); Branschädel et al. (2010); Chien et al. (2013); Gruss et al. (2018) or some other real-space-like construction (e.g., one spatial dimension with energetically tapered boundaries Dias da Silva et al. (2008); Branschädel et al. (2010)). Essentially, this amounts to considering the reservoir Hamiltonian

[TABLE]

and similarly for $H_{\mathcal{R}}$. For simplicity, we take the hopping and chemical potential to be uniform within each reservoir, and take the same $\omega_{0}$ for both $\mathcal{L}$ and $\mathcal{R}$. The $b_{j}$ ($b_{j}^{\dagger}$) are the creation (annihilation) operators in $\mathcal{L}\mathcal{R}$ at the real-space site $j$.

When the system is a single Anderson impurity with equal coupling $v$ to both reservoirs, the remaining Hamiltonians are

[TABLE]

and

[TABLE]

where all indices now explicitly include spin $\sigma$ [Eqs. (S1)–(S4) have spin implicit], $n_{\sigma}=c_{\sigma}^{\dagger}c_{\sigma}$ is the number operator on the system site with spin $\sigma$, and $l=1_{\mathcal{L}(\mathcal{R})}$ is the site in $\mathcal{L}$ ($\mathcal{R}$) that contacts $\mathcal{S}$. For the specific simulations in this work, we will use this model and vary $\omega_{S}$, $v$, and $U$. However, we will work in the single-particle eigenbasis of each of these reservoirs separately. Thus, the spatial nature of this lattice will be inconsequential to the general considerations in our work (it only will change the band structure and the dispersion of the coupling). The computational approach can thus handle non-interacting reservoirs in any dimension, 1D, 2D, 3D, etc., and with long-range hopping.

The canonical transformation that defines the eigenbasis for the 1D reservoir model is $a_{k}=\sum_{j\in\mathcal{L}}\mathcal{U}^{\dagger}_{kj}b_{j}$ with $\mathcal{U}^{\dagger}_{kj}=\sqrt{2/\left(N+1\right)}\sin\left(jk\pi/\left(N+1\right)\right)$ and $k=1,\ldots,N$ for an $N$ site reservoir, yielding

[TABLE]

and $v_{k}=v\,\mathcal{U}^{\dagger}_{k1}$ in Eq. (S2) (and similarly for $\mathcal{R}$). The simulation technique will work in the more general setting where $\mathcal{S}$ is an arbitrary interacting impurity with many electronic sites and vibrational modes, although the computational cost will depend on the Hilbert space dimension and structure of $\mathcal{S}$.

To drive a current, we consider the $\mathcal{L}\mathcal{S}\mathcal{R}$ system initially in contact and in its ground state at zero temperature. At time $t=0$, a bias $\mu_{\mathcal{L}}=-\mu_{r}=\mu/2$ turns on, generating a current. An alternative case is to have $H_{\mathcal{I}}=0$ and $\mu_{\mathcal{L}}=-\mu_{r}=\mu/2$ initially, then turn on $H_{\mathcal{I}}$ and off $\mu$. This starts the system with a density imbalance that drives the current when the applied chemical potential no longer sustains the imbalance. A third case is to have $H_{\mathcal{I}}$ on initially and also the chemical potential drop, letting the latter go to drive the current. These lead to different time dynamics and initial entanglement, but they yield the same steady state and asymptotic entanglement growth Chien et al. (2014). The current from left to right is

[TABLE]

where $n_{\mathcal{L}}$ is the total particle number on the left reservoir, $\Im$ is the imaginary component, and again spin is implicit in the label $k$. In all simulations, the overall filling is determined by the initial state. It is set at (almost) half-filling with $N$ electrons for a system of $N+1+N$ modes (per spin channel).

For any finite system, one can directly simulate the dynamics of non-interacting electrons by evolving the correlation matrix Elenewski et al. (2017)

[TABLE]

where $\rho$ is the full density matrix and $d^{\dagger}_{m}$ ($d_{n}$) are the creation (annihilation) operators at mode $m,n\in\mathcal{L}\mathcal{S}\mathcal{R}$. Defining the single-particle Hamiltonian $\bar{H}$ through

[TABLE]

and using that $\mathrm{tr}\left(d^{\dagger}_{m}d_{n}[H,\rho]\right)=[\bar{H},\mathcal{C}]_{mn}$, the evolution of the correlation matrix is

[TABLE]

This equation can be evolved directly. The dynamics can also be simulated by diagonalizing the “small” dimensional $\bar{H}$ and then transforming the correlation matrix into its eigenbasis Chien et al. (2014).

We first consider the fully non-interacting model [$U=0$ in Eq. (S5)], dropping also the spin since there is no interaction between spin channels (we will multiply the current by an additional factor of two to account for spin, which is not the factor of two already appearing in Eq. (S8)). We then consider the case with interactions.

Figure 2 of the main text shows a matrix product state (MPS) and an exact (via the correlation matrices) simulation of transport in the non-interacting model, Eq. (S5) with $U=0$, using the spatial basis, as has been done in prior work. The steady-state particle current, $I_{\mathrm{ss}}$, is given by Landauer’s formula (regardless of the protocol for driving the current),

[TABLE]

where we explicitly include the factor of two out front (cancelling a factor of 1/2) to account for both spin channels and the second line is in linear response. This equation also defines the exact value of conductance $g^{\rm ex}$ for non–interacting case. The $f_{\mathcal{L}(\mathcal{R})}$ are the Fermi-Dirac distributions in the left (right) reservoir. The transmission function is

[TABLE]

with the retarded Green’s function of the impurity $G(\omega)=1/[\omega-\omega_{S}-\Sigma_{\mathcal{S}\mathcal{L}}-\Sigma_{\mathcal{S}\mathcal{R}}]$, the spectral function of the couplings $\Gamma_{\mathcal{L}(\mathcal{R})}=-2\Im\Sigma_{\mathcal{S}\mathcal{L}(\mathcal{S}\mathcal{R})}$, the self-energies $\Sigma_{\mathcal{S}\mathcal{L}(\mathcal{S}\mathcal{R})}=v^{2}G_{\mathcal{L}(\mathcal{R})}(\omega)$, and the reservoir Green’s functions $G_{\mathcal{L}(\mathcal{R})}(\omega)=1/[\omega-\mu_{\mathcal{L}(\mathcal{R})}-\Sigma_{\mathcal{L}(\mathcal{R})}(\omega)]$, and self-energies $\Sigma_{\mathcal{L}(\mathcal{R})}=(1/2)\left[\omega-\mu_{\mathcal{L}(\mathcal{R})}-\imath\sqrt{4\omega_{0}^{2}-(\omega-\mu_{\mathcal{L}(\mathcal{R})})^{2}}\right]$.

The response of the total system to the driving force results in a rapid rise of the current from zero as particles flow from one reservoir to the other, going into oscillations (due to the presence of Gibbs phenomenon) that decay as the current goes into a quasi-steady state Zwolak (2018). With a large matrix product dimension, the current from the MPS simulation will match the exact solution reasonably well until it abruptly fails for the spatial basis. The origin of the failure is the scattering nature of the problem: Particles come in from, e.g., the left, scattering off the interface at the impurity, generating entanglement between the two reservoirs in the process.

Figure S1 shows the maximum bipartite entanglement entropy across the lattice for position basis and different orderings in mixed basis. The “natural” ordering, with reservoir modes paired with their closest frequency mode on the opposing side and with the system placed around the Fermi level, has the smallest entanglement. We choose to consider canonical transformations that are only permutations to ensure that the matrix product operator that defines the evolution is low dimensional (see below), and, in particular, does not grow with time. We note that for the particular case with a fully non–interacting model, the global single-particle eigenbasis of the Hamiltonian for the time dynamics has zero entanglement growth during evolution. However, the MPO dimension for the initial Hamiltonian grows linearly with the total lattice size (and rotating the interacting component further adds to the MPO dimension). Having interacting systems in mind, including ones with larger interacting regions $\mathcal{S}$, we limit ourselves to the most natural mixed basis, which allows both for a simple MPO and – via a proper ordering – representation of the entanglement structure, in correspondence with the guiding principle mentioned in the main text. In Fig. S2, we show the results, including the entanglement entropy, of the mixed basis simulations for other values of $\omega_{S}$, $v$, $\mu$, and $U=0$ (complementary to Fig. 3 of the main text).

Many-body simulation. Various approaches can perform computations with matrix product states, such as the density matrix renormalization group (DMRG) White (1992) for obtaining the ground state, and the time-evolving block decimation algorithm (TEBD) Vidal (2003, 2004), Krylov- and expansion-based methods Schmitteckert (2004); Zaletel et al. (2015) and the time-dependent variational principle (TDVP) Haegeman et al. (2011, 2016) for simulating time evolution. In order to simulate the time evolution, we employ TDVP for matrix product states (MPS) Haegeman et al. (2011, 2016), where for simulating the effects of local gates we apply the Krylov-based method of Ref. Niesen and Wright, 2012 that is adaptive both in the timestep and in Krylov-subspace dimension. TDVP provides a means to tackle a broad class of Hamiltonians represented as a matrix product operator (MPO). For instance, the Hamiltonian in Eq. (S1) – limited to a single spin channel ($U=0$) in the mixed energy-spatial basis, Eqs. (S3) and (S5), and a single system site interacting with the reservoirs – can be expressed as an MPO with a small bond dimension, $M=4$ (for the spatial Hamiltonian in Eq. (S4) and without interactions in $\mathcal{S}$, the bond dimension is also 4). The latter is independent on the particular ordering of the energy modes. Each additional system site interacting with the reservoir would increase this bond dimension by 2. For a single channel and a single site in $\mathcal{S}$, the exact form for the mixed basis is

[TABLE]

where the sites are ordered according to $\omega_{k}$ form Eq. (S7) (jointly for $\mathcal{L}$ and $\mathcal{R}$). The initial state is the ground state of the Hamiltonian with $\mu_{L}=\mu_{R}=0$, which has an MPO of the same form (with sites ordered using the nonzero value of $\mu$ for $t>0$ so that exactly the same basis is used for the initial state and the subsequent time evolution). Finally, the MPO matrices read

[TABLE]

[TABLE]

and

[TABLE]

The terms equal to zero have been left blank to show the sparsity of the these matrices. Finally, the first and the last matrix in Eq. (S15) (i.e., the smallest and largest energy modes) are limited only to the first row and column, respectively. For such a setup, one can also employ the Jordan-Wigner transformation to the pseudo-spin operators. The fermionic nature of the problem is then reflected by, among other things, $\sigma^{z}$ operators replacing the identities $\mathds{1}$ that connect separate creation and annihilation operators. For the case of interacting reservoirs in Fig. 3(c,d) of the main text, the additional contribution to the Hamiltonian reads $\hbar U_{r}\sum_{j\in\mathcal{L}\mathcal{R}}(n_{j}-1/2)(n_{j+1}-1/2)$. Its MPO in the mixed basis is generated starting from the simple MPOs representing elementary operators. The full Hamiltonian is obtained by subsequent multiplication and addition (as well as bond dimension compression of the resulting MPOs) of the elemental components, using the standard calculus of matrix product states, see e.g., Ref. Schollwöck (2011). Its bond dimension grows approximately as $\sim 24N$.

For the simulations in Fig. 3 of the main text, we set a threshold on the Schmidt values kept, $s_{min}$. The bond dimension is limited by a number of Schmidt values larger than $s_{min}$ and a maximal bond dimension $D$ – whichever is smaller. We use $s_{min}=10^{-6}$ in most of the simulations, which we check is small enough not to influence the results. In the energy representation, the modes outside of the bias window ($-\mu/2$ to $\mu/2$) remain weakly entangled. For that reason, setting $s_{min}$ leads to a small bond dimension in that region and greatly speeds up the simulations, as indicated in the main text. Only the modes in the bias window are getting entangled, and the precision of the simulation for longer times is ultimately controlled by $D$. We show typical behavior of the errors in Fig. S3 and provide some data on convergence of the conductance in Fig. 4(d) of the main text in Fig. S4. We discuss further details of the simulations at the end of this Supplementary Material.

There are alternative setups/structures to handle the time dynamics in the energy basis. Multi-configuration time-dependent Hartree methods employ conditional states on the impurity Meyer et al. (1990); Beck et al. (2000). This was recently employed with MPS for bosonic baths Kloss et al. (2019). Thus, in addition to the physical structure of entanglement addressed here, there are also questions regarding the optimal implementation, which we will examine in a further contribution. We employ a standard MPS structure as we conjecture it will be the most scalable when going to larger-dimensional systems $\mathcal{S}$.

Effective entanglement entropy. To facilitate the comparison between the computational cost of the spatial and mixed bases, we define an effective measure of entanglement of the lattice of a single channel,

[TABLE]

where the sum is over all bipartite cuts (see Eq. (S17) and comment below it for $S_{n}$), and $L-1$ is the total number of relevant cuts. The definition stems from the fact that maximal amount of describable entanglement scales logarithmically with the bond dimension $D_{n}$, and that the leading computational cost of simulation local step of the time evolution is proportional to $D_{n}^{3}$ (simplifying here that the neighboring bonds have the same $D_{n}$ and limiting ourselves to the 1D ordering). As such, equation (S16) incorporates that the required $D$ is inhomogeneous across the cuts and includes, heuristically, how the entanglement entropy contributes to computational cost.

Comments on TDVP simulations. The TDVP procedure, which we rely on to simulate the time evolution of MPS, is masterfully explained in Ref. Haegeman et al., 2016. We briefly summarize it here in order to outline two variations which we employ in this article: combining 1-site and 2-site TDVP updates to, for efficiency, enlarge the bond dimension only when necessary and simulating the time evolution on the $X$-lattice in Fig. 4(a) (or, more generally, on a tree).

Let’s consider a quantum state represented as an MPS of length $L$,

[TABLE]

with dangling legs corresponding to local physical degrees of freedom and connected lines corresponding to virtual degrees of freedom of the tensor network. Above, we depict two representations of the same state in different mixed-canonical forms. All MPS tensors $A_{i}$ to the left (right) of the $n^{\mathrm{th}}$ bond are in left (right) canonical form Schollwöck (2011) – marked here using right-pointing (left-pointing) triangles. For instance, the entanglement of bipartite cuts appearing in Eq. (S16) is fully encoded in the singular values of $C_{n}$ – which we mark here as $\Lambda_{n,i}$ – and reads $S_{n}=-2\sum_{i=1}^{D_{n}}\Lambda_{n,i}^{2}\log_{2}{\Lambda_{n,i}^{2}}$ (for a normalized state and bond dimension $D_{n}$ of the $n^{\mathrm{th}}$ cut). We also define $s_{n}=\min_{i}\Lambda_{n,i}$ as the smallest singular value for bond dimension $D_{n}$.

We can now consider the Hamiltonian $H$ and its expectation value in the state $|\Psi[A]\rangle$,

[TABLE]

The 0-site effective Hamiltonian $H_{\rm eff}^{C_{n}}$, related with the central block $C_{n}$, follows from the expectation value of $H$ in $|\Psi[A]\rangle$, which is calculated/contracted all the way except the contribution from $C_{n}$. All the MPS tensors which are contracted to form $H_{\rm eff}^{C_{n}}$ are in proper left and right canonical forms with respect to the position of the central block $C_{n}$. Similarly, one introduces the 1-site effective Hamiltonian $H_{\rm eff}^{A_{n}^{c}}$ related with the (central) MPS tensor $A_{n}^{c}$, and 2-site effective Hamiltonian $H_{\rm eff}^{AA^{c}_{n,n+1}}$ for two adjacent MPS tensors $A_{n}^{c}$ and $A_{n+1}^{c}$ blocked together.

In order to simulate the Schrödinger equation, TDVP projects the action of $H$ on $|\Psi[A]\rangle$ on the tangent space of $|\Psi[A]\rangle$ Haegeman et al. (2011). To efficiently integrate it, the evolution operator is approximately decomposed into the set of local gates, which are used to update central blocks/sites Haegeman et al. (2016),

[TABLE]

We note, again, that above the MPS is in a correct mixed canonical form with respect to the updated elements. In practice, one does not calculate the matrix representation of the effective Hamiltonian, but, for efficiency, employs a Krylov-based procedure (we use the method in Ref. Niesen and Wright (2012)), which requires only the action of the effective Hamiltonian on a trial vector. The latter can be efficiently calculated by combining smaller building blocks (environments) when $H$ is given as an MPO (or as a sum of local terms or MPOs).

Ref. Haegeman et al., 2016 discusses two main decompositions to simulate the time evolution,

[TABLE]

For the 1-site TDVP scheme

[TABLE]

with central sites $A^{c}_{n}$ evolved forward in time, and central blocks $C_{n}$ evolved backward in time. $T^{[1]}$ constitute one sweep from left to right, where, before each local unitary update, MPS is put in a proper mixed canonical form. It is completed by its adjoint $T^{[1]*}$ with all the gates applied in the reverse order, i.e., a sweep from right to left, making it a $2^{\mathrm{nd}}$ order method in $dt$. It operates with fixed bond dimensions (at each cut) of the MPS.

Dynamical adjusting of the bond dimensions can be obtained by 2-site TDVP scheme

[TABLE]

(plus its adjoint $T^{[2]*}$ in the reverse order). Now, the 2-site gate are evolved forward in time, and 1-site gate are evolved backward in time. In case of the 2-site update $U^{[2]}$, two adjacent MPS matrices are blocked together and subsequently split using a singular value decomposition (SVD), truncating the virtual bond to given size/weights. It is, however, numerically significantly more costly both due to larger vectors appearing in the Krylov procedure and the additional SVD.

In this article, we employ a slight modification of the above procedures by combining the two schemes. The 2-site gates are employed only locally when both are necessary (all the Schmidt values of a given cut above a threshold $s_{min}$) and possible (bond dimension of a given cut is below the maximal $D$). Such an approach is consistent with the strongly inhomogeneous nature of the system we consider where entanglement is localized only in a part of the system, and the MPS bond dimensions between modes outside of the bias window can remain small. To that end, it is sufficient to note how to transition between parts of $T^{[1]}(y)$ and $T^{[2]}(y)$ during a sweep to build $T^{[{\rm mixed}]}(y)$. If the $n$-th bond is enlarged and the next one is not, one gets $T^{[{\rm mixed}]}(y)=\ldots U^{[0]}_{n+1}(-y)U^{[2]}_{n,n+1}(y)\ldots$. On the other hand, if $n^{\mathrm{th}}$ bond is not enlarged and the next one is, one has $T^{[{\rm mixed}]}(y)=\ldots U^{[2]}_{n+1,n+2}(y)U^{[0]}_{n}(-y)\ldots$. We perform the adjoint (sweep from right to left) using the exact reversal of the gates in $T^{[{\rm mixed}]}(y)$. We numerically observe, however, that finding new $T^{[{\rm mixed}]*}(y)$ based on the bond dimension/Schmidt weights does not reduce the order of the method in $dt$. For clarity, below we collect the full procedure as a pseudocode 1.

Finally, the tree-tensor network ansatz corresponding to the X-shape lattice in Fig. 4(a) of the main text is

[TABLE]

Again, we show two different mixed canonical representations of the same state with all the tensors to the left (right) of the central site [bond on the right-hand side] in the left (right) canonical form. On the right-hand side, we show the central block $C_{S}$ between the two spin channels, which in our case corresponds to the placement of the impurity $\mathcal{S}$ (hence the index).

The simulation of the time evolution is obtained using $U(dt)\approx T^{X}(dt/2)T^{X*}(dt/2)$. A one-way sweep is composed as $T^{X}(y)=T^{\uparrow[mixed]}(y)U^{[0]}_{\mathcal{S}}(-y)T^{\downarrow[mixed]}(y)$. Sweeps of the spin channels, $T^{\uparrow(\downarrow)[mixed]}(y)$, are done similarly as for the 1D chain above, and $U^{[0]}_{\mathcal{S}}$ describes the update of the central block $C_{S}$. If one wants to enlarge that bond, one can replace $U^{[0]}_{\mathcal{S}}(-y)$ with $U^{[1]}_{\downarrow\mathcal{S}}(-y)U^{[2]}_{\mathcal{S}}(y)U^{[1]}_{\uparrow\mathcal{S}}(-y)$, with $U^{[2]}_{\mathcal{S}}$ acting on two tensors connecting the spin channels, and $U^{[1]}_{\uparrow(\downarrow)\mathcal{S}}$ on each of them. We depict the order of the full sweep – combining to $U(dt)$ – with blue arrows on the right-hand side of Eq. (S25). Its symmetric form ensures that it is second order in $dt$. The total Hamiltonian (which generates the $U$ gates) is treated as a sum of (the contributions coming from) MPOs for each channel and the interacting term, $Un_{\downarrow}n_{\uparrow}$, which is coupling the system modes placed next to each other in the $X$-lattice geometry.