Accumulative reservoir construction: Bridging continuously relaxed and periodically refreshed extended reservoirs

The simulation of open many-body quantum systems is challenging, requiring methods to both handle exponentially large Hilbert spaces and represent the influence of (infinite) particle and energy reservoirs. These two requirements are at odds with each other: larger collections of modes can increase the fidelity of the reservoir representation but come at a substantial computational cost when included in numerical many-body techniques. An increasingly utilized and natural approach to control the growth of the reservoir is to cast a finite set of reservoir modes themselves as an open quantum system. There are, though, many routes to do so. Here, we introduce an accumulative reservoir construction—an ARC—that employs a series of partial refreshes of the extended reservoirs. Through this series, the representation accumulates the character of an infinite reservoir. This provides a unified framework for both continuous (Lindblad) relaxation and a recently introduced periodically refresh approach (i.e., discrete resets of the reservoir modes to equilibrium). In the context of quantum transport, we show that the phase space for physical behavior separates into discrete and continuous relaxation regimes with the boundary between them set by natural, physical timescales. Both of these regimes “turnover” into regions of over- and underdamped coherence in a way reminiscent of Kramers' crossover. We examine how the range of behavior impacts errors and the computational cost, including within tensor networks. These results provide the first comparison of distinct extended reservoir approaches, showing that they have different scaling of error versus cost (with a bridging ARC regime decaying fastest). Exploiting the enhanced scaling, though, will be challenging, as it comes with a substantial increase in (operator space) entanglement entropy.

Figure 1

Accumulative reservoir construction (ARC). ARC bridges two different approaches to simulate open quantum systems, from continuous relaxation (CR) to periodic refresh (PR). It consists of two steps of time $\tau$. The first is a unitary evolution with the composite Hamiltonian (highlighted in green on the left) of the system $\mathcal{S}$ and reservoirs $\mathcal{L}\mathcal{R}$, the reservoirs being discretized into a finite number of modes $N_{\mathcal{W}}$. The second step is a dissipative process that partially relaxes the reservoir modes towards their isolated equilibrium states (highlighted in blue on the right). In our implementation the dissipation is Lindblad damping with a homogeneous, intramode relaxation strength $\gamma$ and acts only on the reservoirs. ARC is shown in the context of quantum transport, where two reservoirs at different chemical potentials and/or temperatures drive a current $I$ through an impurity or extended system $\mathcal{S}$. Depending on $\gamma$, $\tau$, and $N_{\mathcal{W}}$, the NESS from ARC can be a good approximation to the NESS from macroscopic reservoirs. CR corresponds to the limit $\tau \to 0$ at fixed $\gamma$, whereas PR corresponds to $\gamma \to \infty$ at fixed $\tau$.

Figure 2

Resonant-level model. The schematic shows the example—the resonant level model—used to numerically characterize ARC. We will use the parameters ${\omega }_{\mathcal{S}}=v_{\mathcal{S}}={\omega }_0/2$, several different temperatures but with $T_{\mathcal{L}}=T_{\mathcal{R}}=T$, and the chemical potentials ${\mu }_{\mathcal{L}}=-{\mu }_{\mathcal{R}}={\omega }_0/4$ (i.e., a total potential drop of ${\omega }_0/2$ applied symmetrically). By going to the single particle eigenbasis of $\mathcal{L}$ and $\mathcal{R}$, the model conforms to the setup in Fig. 1. Even within tensor networks, we simulate the reservoirs using their eigenbasis and thus their spatial dimension is not directly relevant.

Figure 3

Time dynamics of ARC. The plots show the current between $\mathcal{L}$ and $\mathcal{S}$ versus time for $N_{\mathcal{W}}=128$ sites per reservoir. The ARC parameters are (a) $\tau \to 0$ at finite $\gamma \approx 0.084{\omega }_0$ (i.e., in the CR regime), (b) $\gamma \tau =0.1$ and $\tilde{\tau }\approx 37.90/{\omega }_0$, (c) $\gamma \tau =1$ and $\tilde{\tau }\approx 37.30/{\omega }_0$, and (d) $\gamma \to \infty$ at finite $\tau =64/{\omega }_0$ (i.e., in the PR regime). The red points indicate the current right before the application of dissipation, whereas CR is a continuous red line since relaxation is continuous. For all regimes of ARC except the CR limit, the stroboscopic points in red are the most relevant. The convergence of these at long times gives the estimate of the current within the exact NESS. We note that (a) $\gamma$ for CR is provided by the heuristic procedure of Ref. [55], [(b) and (c)] $\tilde{\tau }$ is given by the value that best estimates the current at fixed $\gamma \tau$, and (d) $\tau$ is given by the (physical) refresh time for PR, $\tau ={\tau }_{\mathcal{W}}=N_{\mathcal{W}}/2{\omega }_0$ with $2{\omega }_0$ the Fermi velocity.

Figure 4

Phase diagram of ARC. (a) The schematic phase diagram (of the NESS current) shows several regimes as the total relaxation $\gamma \tau$ and the action $\tilde{\tau }=|\tau -ı2/\gamma |$ vary at fixed $N_{\mathcal{W}}$. In this parametrization, the PR limit ($\gamma \to \infty$ at fixed $\tau$) is on the upper boundary and the CR limit ($\tau \to 0$ at fixed $\gamma$) on the lower. Sweeping the action from small to large values shows three general regimes. On the left is an overdamped (Zeno) regime (i.e., large $\gamma$ and small $\tau$ corresponding to a small action). On the right is an underdamped regime. The weak relaxation regime on the bottom right (i.e., with CR-like behavior) leads to currents limited by the relaxation strength. In the overly coherent regime in the top right (i.e., with PR-like behavior), there are oscillations across the $\mathcal{L}\mathcal{S}\mathcal{R}$ system. The middle of the diagram contains the physical regimes, with a CR-like physical regime at the bottom and PR-like at the top. The physical regimes are separated by a region of an overly coherent $\mathcal{S}$. (b) The normalized current for $N_{\mathcal{W}}=256$ on the same diagram. The lines of constant $\gamma$—the ones coming vertically from the bottom—are ${\omega }_0/43.3$ and ${\omega }_0$ which are heuristic optimal for CR and transition to overdamped regime for CR, respectively. These are normal turnover points discussed extensively elsewhere [20, 21, 49, 54]. The solid lines of constant $\tau$—the ones coming vertically form the top—are from left to right, ${\tau }^{★}=\pi /\mathcal{W}$, $2{\tau }_{\mathcal{S}}$, and ${\tau }_{\mathcal{W}}$. The former two define the transition regime. The dashed line is ${\tau }_{\mathcal{S}}$, showing it also demarcates features. The ${\tau }_{\mathcal{W}}$ shows the transition to an overly coherent $\mathcal{L}\mathcal{S}\mathcal{R}$. To the right of ${\tau }_{\mathcal{W}}$, the current remains a good approximation until $2{\tau }_{\mathcal{W}}$ since that is when the density wavefront returns to the system after bouncing off the finite reservoir boundary. (c) The convergence time to NESS follows from Eq. (36). For another perspective on the phase diagram, Fig. 14 shows it within the $\gamma -\tau$ parametrization.

Figure 5

Convergence and computational cost of ARC. The three panels show the (a) current error, Eq. (44), (b) the correlation matrix error, Eq. (45), and (c) the estimated computational cost, Eq. (46), for our example system and for $N_{\mathcal{W}}=512$, $\beta {\omega }_0=40$, and $\mu ={\omega }_0/2$. The two error plots have many of the same features as the phase diagram for ARC, Fig. 4, as we expect since the different regimes directly impact the error. Using the current error as the quantity to be optimized, the best choice of action for a fixed $\gamma \tau$ is marked by the white points in (a). These points are replicated in the other plots [(b) and (c)] to see where they fall for other quantities. These points are close to where the correlation matrix error is minimum in (b), but they do not follow any features in the computational cost in (c). As is evident, CR- and PR-like regimes have very different estimated computational costs. Within the physical regimes, this cost is correlated with the error, with cost going up for more accurate simulations. Outside of these regimes, though, there is no correlation. The whites dots for $\gamma \tau =0.1,1$, and 10 are made thicker to guide the readers eye. Vertical solid line represents $\tilde{\tau }$ equal to heuristic choice for CR, and vertical dashed line represents $\tilde{\tau }={\tau }_{\mathcal{W}}$, which is the physical choice for PR. We employ these $\tilde{\tau }$ for the CR and PR limits.

Figure 6

Mixed basis. For MPS simulations and the calculation of OSEE, we group the modes from reservoirs $\mathcal{L}$ and $\mathcal{R}$, in the frequency basis, and ordered them according to ${\omega }_k$ in Eq. (2). It reflects the emergent structure of correlations in the reservoirs [73]. We place the system modes around zero frequency. We indicate the bipartition of the lattice for calculating OSEE.

Figure 7

Error-cost scaling for ARC. The error is shown versus the estimated computational cost, Eq. (46), for (a) a low and (b) a high temperature both with bias $\mu ={\omega }_0/2$. For a given $\gamma \tau$, the reservoir size $N_{\mathcal{W}}$ increases from left to right (as indicated on the top of the figure), which is the only control parameter here since the estimated computational cost comes from noninteracting simulations. For ARC, the action is taken at the white points in Fig. 5. For CR, we choose the relaxation from heuristic approach [55]. For PR, we choose $\tau ={\tau }_{\mathcal{W}}$. While there are some transitions in the error-cost scaling, most of the data settles into a power law behavior, Eq. (48). The dashed lines show fitted scaling for CR (red), ARC $\gamma \tau =1$ (green), and PR (blue). Dotted lines mark related scaling laws. These are not fits but a reference lines of the same $\nu$ as their dashed counterparts. More information on the fit can be found in Sec. 5e and Appendix pp3-s5.

Figure 8

Thermal correlation length in the PR regime. We show the relative current error rescaled by $N_{\mathrm{th}}^2$ versus the number of reservoir modes divided by $N_{\mathrm{th}}$. We observe the collapse of the data for a range of temperatures onto a single curve, that exhibits two scaling regimes, one for $N_{\mathcal{W}}\lesssim N_{\mathrm{th}}$ and the other for $N_{\mathcal{W}}\gg N_{\mathrm{th}}$. In both regimes, the error vanishes as $N_{\mathcal{W}}^{-2}$, but with different prefactors.

Figure 9

Convergence with bond dimension. The plots show the error versus cost obtained using the MPS implementation as described in the main text. The data were obtained for (a) CR and (b) ARC $\gamma \tau =1$ at low and high temperature as indicated on the panels. The matrix product dimension, $D$, is increased at fixed $N_{\mathcal{W}}$ (shown in the figure), starting from $D=64$ (at the left) and increasing by 64 between consecutive points. While technically still an estimated computational cost, $C^{\mathrm{MPS}}={\tau }_c·N_{\mathcal{W}}·D^3$, this is the actual convergence of the MPS using $D$ from simulation (${\tau }_c$ is still from the exact evolution operator). The solid circle points are from the estimated costs using a scaled version of Eq. (46), $c·N_{\mathcal{W}}·{\tau }_c·2^{3·S_O}$. We take $c={10}^6$ as a uniform factor to match where the given $D$ is sufficient to achieve convergence at a set $N_{\mathcal{W}}$. The error of the current in MPS calculation is estimated by an average ${\sigma }_I$, as defined in Eq. (41), over a time window of duration $20/{\omega }_0$ for CR or 2 to 5 cycles for ARC. Time evolution uses an integration step of $0.125/{\omega }_0$. The refresh is applied in a single step rather than via time evolution.

Figure 10

CR: Kramers' turnover and estimating the optimal $\gamma$. The relaxation rate significantly influences the steady-state current. There are always three distinct regimes of behavior [20]. Examining $I/I^{\circ }$ versus $\tilde{\tau }$, these are seen here as the current initially increasing linearly with $\tilde{\tau }\sim {\gamma }^{-1}$ (overdamped), plateauing (physical), and then decreasing as $1/\tilde{\tau }\sim \gamma$ (relaxation-limited). There also can be anomalous regimes depending on the particular model [54], for which we see here a small Markovian anomaly on the left hand side of the plateau. In order to estimate the optimal $\gamma$ (or $\tilde{\tau }$), one can shift the alignment of modes in $\mathcal{L}$ and $\mathcal{R}$, which removes virtual tunneling between them. The $\mathcal{L}\mathcal{R}$ in alignment and shifted intersection provides a good estimate of the physical current for given $N_{\mathcal{W}}$. The heuristic approach avoids unphysical features by targeting a $\gamma$ that is sufficiently large that the discreteness of the modes is not influencing the transport but allows coherence to develop [55]. The plot uses $N_{\mathcal{W}}=512$ sites per reservoir, temperature $\beta {\omega }_0=40$, and bias $\mu ={\omega }_0/2$.

Figure 11

Typical regimes for PR. (a) The stroboscopic steady-state current vs $\tau /{\tau }_{\mathcal{S}}$ for three system sizes (giving three different ${\tau }_{\mathcal{S}}$). For very small $\tau$, the current rises to the steady-state level. If the refresh period $\tau$ is comparable to any timescale in $\mathcal{L}\mathcal{S}\mathcal{R}$, then coherent oscillations become apparent. For $\tau <{\tau }_{\mathcal{S}}$, for ${\tau }_{\mathcal{S}}$ the $\mathcal{S}$ timescale, we observe features in the steady-state current related to evolution inside $\mathcal{S}$, and distinct features for multiples of ${\tau }_{\mathcal{S}}$ that match refresh time with wavefront traversing the $\mathcal{S}$ from one side to the other. The features become less apparent as $\tau$ increases, eventually becoming negligible. The vertical dotted and dashed lines in (a) mark ${\tau }_{\mathcal{S}}$ and $2{\tau }_{\mathcal{S}}$, respectively. (b) The stroboscopic steady-state current vs $\tau /{\tau }_{\mathcal{W}}$ for the same set of system sizes. The current remains stable even after reaching the reservoir timescale ${\tau }_{\mathcal{W}}$. It is only when the wave front returns to the interface at $2{\tau }_{\mathcal{W}}$ that the current reverses, for which there will be further reversals at integer multiples of $2{\tau }_{\mathcal{W}}$ (i.e., $4{\tau }_{\mathcal{W}}$, etc.). The vertical dashed lines in (b) mark multiples of ${\tau }_{\mathcal{W}}$. After reversal, the timescale ${\tau }_{\mathcal{S}}$ determine when the current stabilizes. Unlike in the main text, $\mathcal{S}$ is a uniform lattice of $N_{\mathcal{S}}$ sites with ${\omega }_{\mathcal{S}}=0$ and with nearest-neighbor hopping $v_{\mathcal{S}}={\omega }_0/2$. The data are obtained for reservoirs with $N_{\mathcal{W}}=256$ modes each, temperature $\beta {\omega }_0=40$, and symmetrically applied bias $\mu ={\omega }_0/2$.

Figure 12

PR errors for increasing ${\mathit{N}}_{\mathcal{W}}$. (a) The curves show the relative error versus $\tau$ for many different $N_{\mathcal{W}}$. As $N_{\mathcal{W}}$ increases, the plateau for physical behavior increases. As well, more reservoir modes increase the accuracy of the current. The vertical dotted line shows the impurity timescale ${\tau }_{\mathcal{S}}$ and the vertical dashed line shows thermal time ${\tau }_{\mathrm{th}}=N_{\mathrm{th}}/2{\omega }_0$ for a reservoir with thermal relaxation time $\beta {\omega }_0=40$. (b) The curves show the same results as in (a) but versus the computational cost, Eq. (46).

Figure 13

Steady-state current for ARC. (a) The curves show the stroboscopic steady-state current for ARC vs $\tilde{\tau }$ for a few $\gamma \tau$. All the parameter regimes show a Kramers'-like increase and then plateau for small to moderate $\tilde{\tau }$. However, for larger $\tilde{\tau }$ and $\gamma \tau$ the curves escape the CR-like phase and coherent oscillations become apparent. For very large $\tilde{\tau }$, all curves enter the underdamped regime, which can be either directly similar to Kramers' rate-limiting process (i.e., a linear dependence on $\gamma$) or reflect the stroboscopic, PR-like current, which breaks down and eventually reverses. For $\gamma \tau =1$, however, the underdamped regime combines these two features and the current goes down with $\tilde{\tau }$ but still undergoes some rapid oscillation and eventually reverses when the wavefront returns back to the interface. (b) The curves show ARC with the same three $\gamma \tau$ versus $\tau$ instead of $\tilde{\tau }$ (the $\gamma$ is changing with $\tau$ since $\gamma \tau$ is fixed). The coherent oscillations for the two larger values of $\gamma \tau$ are seen to occur with the same frequency. Both plots take data from Fig. 4.

Figure 14

ARC phase diagram in $\mathit{\gamma }\text{−}\mathit{\tau }$ space. (a) The abstract phase diagram for ARC. The lines of constant $\tau$ and constant $\gamma$ are now vertical and horizontal, respectively. The PR and CR limits are now on the top and left, respectively. This provides a different perspective on the phase diagram in Fig. 4. [(b) and (c)] The steady-state current and convergence time both on the $\gamma \text{−}\tau$ space. The data and lines of constant $\tau$ and $\gamma$ are at the same values as in Fig. 4.

Figure 15

Phase diagram for ${\mathit{N}}_{\mathcal{S}}=\mathit{16}$. The lines demarcate the boundary and some features of the overly coherent $\mathcal{S}$ regime. Unlike the main text, $\mathcal{S}$ is a uniform lattice with onsite frequencies of zero, a constant nearest-neighbor coupling $v_{\mathcal{S}}={\omega }_0/2$, and system-reservoir coupling of (a) ${\omega }_0$ and (b) ${\omega }_0/2$. The line of fixed $\tau$ equal to ${\tau }^{★}=\pi /\mathcal{W}$ is plotted in white, ${\tau }_{\mathcal{S}}$ plotted in green, and $2{\tau }_{\mathcal{S}}$ plotted in yellow. The $\mathcal{S}$ timescale is ${\tau }_{\mathcal{S}}=N_{\mathcal{S}}/2v_{\mathcal{S}}$. Both examples use $N_{\mathcal{W}}=256$ sites per reservoir, temperature $\beta {\omega }_0=40$, and bias $\mu ={\omega }_0/2$.

Figure 16

Convergence and computational cost of ARC. The three panels show the (a) current error, Eq. (41), (b) the correlation matrix error, Eq. (45), and (c) the estimated computational cost, Eq. (46), for our example system and for $N_{\mathcal{W}}=256$, $\beta {\omega }_0=40$, and $\mu ={\omega }_0/2$. The two error plots have many of the same features as the phase diagram for ARC, Fig. 4, as we expect since the different regimes directly impact the error. Using the current error as the quantity to be optimized, the best choice of action for a fixed $\gamma \tau$ is marked by the white points in (a). These points are replicated in the plot showing the correlation matrix error (b), and estimated computational costs (c) to see where they fall for other quantities. The whites dots for $\gamma \tau =0.1,1$, and 10 are made thicker to guide readers eye. Vertical solid line represents $\tilde{\tau }$ equal to heuristic choice for CR, and vertical dashed line represents $\tilde{\tau }={\tau }_{\mathcal{W}}$, which is a physical choice for PR. We employ these $\tilde{\tau }$ for the CR and PR limits.

Figure 17

Optimal estimates for ARC. The data take ARC optimal points from Fig. 16 for $N_{\mathcal{W}}=256$ and from Fig. 5 for $N_{\mathcal{W}}=512$. From the white points we extract (a) relative error ${\overline{\sigma }}_I$, (b) error for $\mathcal{S}$ correlation matrix, (c) operator space entanglement entropy for central cut, and (d) the convergence time to reach a steady state ${\tau }_c$. The data show the transition between optimal estimates falling into CR-like regime to these which are already in PR-like regime. Vertical dashed lines are included to guide readers eye and mark $\gamma \tau =0.1,1$, and 10, which are highlighted points in Figs. 16 and 5.

Figure 18

Relative error without the moving average. The panels show the (a) current error, ${\sigma }_I$, (b) the correlation matrix error, Eq. (45), for $\mathcal{S}$ as in the main text. The reservoirs are $N_{\mathcal{W}}=512$ sites each, temperature $\beta {\omega }_0=40$, and symmetric bias $\mu ={\omega }_0/2$. The plot takes the data from Fig. 5 but does not apply moving average in Eq. (44) to ${\sigma }_I$. White points are taken from Fig. 5 and mark the position of lowest average error ${\overline{\sigma }}_I$ for each $\gamma \tau$. These plots demonstrate that the averaging procedure does not deform the phase diagram. They also show that the averaging itself does not bias the optimum. However, on the CR side, the optimum from ${\overline{\sigma }}_I$ tracks the accidental crossing. Those crossings are not fully visible on either the PR or CR side of the phase diagram due to the finite resolution of the heat map.