Efficient construction of time-invariant process tensors for simulating high-dimensional non-Markovian open quantum systems

TL;DR

A new algorithm for computing process tensors--matrix product operator (MPO) representations that capture the environment influence--which achieves greatly enhanced computational scalings with system size, while maintaining linear scaling with simulation length is proposed.

Abstract

Numerical methods for obtaining exact dynamics of non-Markovian open quantum systems are mostly limited to either small systems or to short-time evolution only. Here, we propose a new algorithm for computing process tensors--matrix product operator (MPO) representations that capture the environment influence--which achieves greatly enhanced computational scalings with system size, while maintaining linear scaling with simulation length. We build on recent developments in the field which allow for long-time evolutions through process tensors which have a time-translational invariance. These can be built for general Gaussian environments and generic coupling operators with the system using infinite time-evolving block decimation (iTEBD). We introduce a modified iTEBD algorithm using intermediate compression steps which bring down the computation time scaling with system size $d$ from $\mathcal{O}(d^8)$ to $\mathcal{O}(d^4)$, as well as significantly lowering the required memory. To illustrate the power of this method, we apply it to the problem of dispersive qubit readout in circuit QED, which was previously out-of-reach numerically. The full treatment of the measurement resonator, which requires a large system space, combined with the long simulation times precipitated by the separation of timescales between the measurement drive and the environment dissipation, is now possible. The algorithm we introduce not only allows for capturing non-Markovian dynamics in large open quantum systems, but also further extends all the existing capabilities of process tensors, for example in quantum optimal control, or in computation of multi-time correlations or of steady states, to more complex systems with tens of levels.

Figures (10)

• Figure 1: Graphical representation of the tensor networks for the time-invariant process tensors. (a) Two-dimensional infinite tensor network for the infinite influence tensor $\mathcal{F}^{\cdots\mu_0\mu_1\mu_2\mu_3\cdots}$ from link_Open_2024. (b) Equivalent rotated tensor network on which can be applied iTEBD from link_Open_2024. (c) Explicit representation of the $b^{\mu\nu}_{ab}(k)$ tensor with the Kronecker deltas as black dots. (d) Tensor network illustrating the application of the iTEBD contraction scheme to the $\mathcal{F}$ tensor. The unit cell of the iTEBD step is highlighted by a red dashed rectangle. The step consists in applying the $\tilde{b}(k)$ gate to the infinite MPS in canonical form with its $A$, $B$ tensors in blue and the corresponding diagonal weight matrices as red diamonds.
• Figure 2: (a) Tensor network for the gate contraction of the enhanced iTEBD step. (b) Tensor network with SVD performed on the $b(k)$ tensor. (c) Tensor network displaying the partial SVD on the $\theta^A$ and $\theta^B$ blocks highlighted by red dashed rectangles, with respect to the red $d^2$ size leg. (d) Resulting tensor network from the partial SVDs. The final SVD is done on the $\Theta$ block highlighted by a red dashed rectangle.
• Figure 3: Bond dimensions for the last step of the TTI-PT computation. The final TTI-PT bond dimension $\chi$ is shown for the old method in orange and for the new method in blue. The intermediate bond dimensions of the new method $\alpha$, $\beta_1$ and $\beta_2$ are shown respectively with diamonds, down-facing triangles and up-facing triangles. These can be compared with the orange dashed $d^2$ line of the old method. The benchmarks are done for a harmonic oscillator coupled through the operator $\hat{a}+\hat{a}^\dagger$ to an ohmic bath with parameters $\eta=0.01$ and $\omega_c=3$ with iTEBD parameters $\tilde{k}=500$, $\epsilon_\mathrm{rel}=10^{-3}$.
• Figure 4: (a) Memory requirement to store the largest tensor as a function of Hilbert space size $d$. Comparison between the regular method (in orange circles) and the enhanced method (in blue circles). The stars represent the potential memory requirement with matrix-free methods for the $\theta$ SVDs. (b) Dependence of CPU computation time with Hilbert space size $d$ for the full algorithm. Comparison between the regular method (in orange) and the enhanced method (in blue). Approximate scaling for both methods are drawn in black dashed lines. The benchmarks are done for a harmonic oscillator coupled through the operator $\hat{a}+\hat{a}^\dagger$ to an ohmic bath with parameters $\eta=0.01$ and $\omega_c=3$ with iTEBD parameters $\tilde{k}=500$, $\epsilon_\mathrm{rel}=10^{-3}$. The computations were done using an Intel Xeon Gold 5218 Processor @ 3.9GHz with 32 cores.
• Figure 5: Circuit diagram for dispersive qubit readout representing the Hamiltonian \ref{['eq:Hamiltonian']}. The Rabi Hamiltonian $\hat{H}_\mathrm{R}$ given by \ref{['eq:HR']} is drawn in green, the drive Hamiltonian $\hat{H}_d$ given by \ref{['eq:Hd']} is drawn in red, and the readout line Hamiltonian $\hat{H}_\mathrm{E}$ given by \ref{['eq:HE']} is drawn in blue with the corresponding spectral densities pictured below in shades of green to blue for the Purcell filter strengths $p$ studied here: $0$, $0.25$, $0.5$, $0.75$, and $0.9$.