Lindblad equation, quantum computing, open quantum systems, completely positive and trace preserving, Gregory quadrature, nested Picard iteration, low-rank

1 Introduction↩︎

Open quantum systems refer to quantum systems that interact with their surrounding environment and therefore cannot be regarded as isolated. Such systems govern the hardware used in quantum computers. In this setting, the dynamics is governed not only by the Hamiltonian but also by dissipative and decoherence effects induced by the environment. A standard model for such dynamics is the Lindblad equation, which describes the evolution of the density matrix \(\rho(t)\).

The Lindblad equation plays a central role in quantum information processing, quantum control, and the modeling of realistic quantum devices, where dissipation and decoherence are unavoidable. Accurate and efficient numerical simulation of the density matrix is therefore of fundamental importance. However, this task is challenging for at least two reasons. First, the state space grows rapidly with the number of subsystems: for example, for a \(Q\)-qubit system, the density matrix has dimension \(2^Q\times 2^Q\). Second, a numerical method should preserve the key structural properties of the density matrix, such as Hermiticity, positive semidefiniteness, and unit trace. These requirements make the construction of efficient high-order numerical integrators challenging.

Existing methods that are designed to preserve the completely positive (CP) property often do so through constructions motivated by Choi’s theorem; see, for example, [1], [3]–[7]. In particular, a numerical update \[\rho^n \mapsto \rho^{n+1}=\mathcal{G}(\rho^n), \qquad \rho^n \approx \rho(t^n),\] is completely positive if the linear map \(\mathcal{G}\) admits a Kraus representation of the form \[\mathcal{G}(\rho)=\sum_l G_l \rho G_l^{\dagger}.\] The unit-trace condition can then be enforced, when necessary, by the renormalization \[\rho^{n+1}\leftarrow \frac{\rho^{n+1}}{\operatorname{Tr}(\rho^{n+1})}.\]

High-order methods for the Lindblad equation may be constructed either stochastically or deterministically. One stochastic approach is based on trajectory methods, often referred to as high-order unravelings. Such methods evolve random pure states, or random density matrices, associated with individual quantum trajectories determined by measurement records or jump events. The Lindblad solution is recovered only after averaging over these trajectories; see, for example, [8], [9]. Consequently, in addition to time discretization error, they typically incur sampling error. An alternative is provided by deterministic high-order integrators that evolve the density matrix directly. Among them, an important class is based on Duhamel’s principle; see, for example, [1], [3], [4], [10].

In our previous work [1], we proposed high-order low-rank CPTP schemes based on nested Picard iteration (NPI) with Gaussian quadrature rules. The use of Gaussian quadrature was motivated in part by the positivity of its weights, which is crucial for preserving complete positivity; see also [10] which uses Gaussian quadrature for high-dimensional integrals arising in the time-ordered methods of [4], [10]. While structure-preserving, the schemes in [1] are costly in particular in high order cases due to the need to repeatedly recalculate the density matrix at intermediate stages. To compute a solution of order \(p\), the Gaussian-quadrature-based NPI scheme requires \(O(\frac{p!}{2^p})\) flow-operator evaluations.

In this paper, we propose the Gregory NPI schemes to improve the efficiency of [1] by replacing Gaussian quadrature with Gregory quadrature on equidistant nodes. The resulting schemes remain high-order (up to ninth order), low-rank, and structure-preserving, while being significantly more efficient. In particular, to compute a solution of order \(p\), the scheme proposed here requires only \(N_q^{(p)}+1\) flow-operator evaluations, where \(N_q^{(p)}+1\) is the number of Gregory quadrature nodes, which is taken to be \(2p-2\). Therefore, the scheme only requires \(O(p)\) flow operator evaluations, which demonstrates the large reduction in computational cost particularly for high-order cases.

The remainder of this paper is organized as follows. In Section 2, we review the Lindblad equation and Gregory quadrature. In Section 3, we review the Gaussian-quadrature-based NPI schemes. In Section 4, we formulate the new and more efficient method based on Gregory quadrature. In Section 5, we study the stability of the Gregory NPI schemes without trace renormalization for a single-qubit model problem. In Section 6, we describe the low-rank truncation. In Section 7, we present numerical results for representative quantum systems. These include a two-qubit problem with a known exact solution, a controlled qudit-resonator system for examining convergence, low-rank truncation, and the effect of trace renormalization, and a CNOT-gate problem involving two qudits and a resonator. We conclude the paper in Section 8.

2 Preliminaries↩︎

In this section, we introduce the Lindblad equation and Gregory quadrature.

2.1 Governing Equations↩︎

In this work, the evolution of the density matrix \(\rho(t)\in\mathbb{C}^{d\times d}\) is governed by the Lindblad equation \[\begin{align} \label{eqn95master} \frac{d\rho(t)}{dt} &= -i\bigl(H(t)\rho(t)-\rho(t)H(t)\bigr)+\mathcal{L}(\rho(t)),\\ \rho(0) &= \rho_0. \nonumber \end{align}\tag{1}\] Here \(\rho_0\) is the initial density matrix, and \[H(t)=H_d+H_c(t),\] is the system Hamiltonian, where \(H_d\) is the drift Hamiltonian and \(H_c(t)\) is a time-dependent control Hamiltonian.

The dissipative part of the dynamics is described by the Lindbladian \[\label{eq:Lindbladian95general} \mathcal{L}(\rho) = \sum_{\alpha}\mathcal{L}_{\alpha}\rho \mathcal{L}_{\alpha}^{\dagger} -\frac{1}{2}\sum_{\alpha} \left( \mathcal{L}_{\alpha}^{\dagger}\mathcal{L}_{\alpha}\rho + \rho\,\mathcal{L}_{\alpha}^{\dagger}\mathcal{L}_{\alpha} \right),\tag{2}\] where \(\{\mathcal{L}_{\alpha}\}_{\alpha}\) is a family of jump operators and \(\dagger\) denotes the adjoint.

Our main interest is methods for simulating quantum computers and thus we consider a transmon-based model for a composite superconducting quantum device; see [11]–[14]. We adopt the notation used in [14]. For a system consisting of \(Q\) subsystems, where subsystem \(k\) has \(n_k\) energy levels, the drift Hamiltonian is given by \[H_d = \sum_{k=0}^{Q-1} \left( \omega_k a_k^\dagger a_k -\frac{\xi_k}{2}\,a_k^\dagger a_k^\dagger a_k a_k +\sum_{l>k} \Big( J_{kl}(a_k^\dagger a_l+a_k a_l^\dagger) -\xi_{kl} a_k^\dagger a_k a_l^\dagger a_l \Big) \right),\] where \(\omega_k\ge0\) is the \(0\to1\) transition frequency of subsystem \(k\), \(\xi_k\ge0\) is the self-Kerr coefficient, and \(J_{kl}\ge0\) and \(\xi_{kl}\ge0\) denote the coupling and cross-Kerr coefficients between subsystems \(k\) and \(l\), respectively. The lowering operator on subsystem \(k\) is defined by \[a_k = I_{n_0}\otimes\cdots\otimes I_{n_{k-1}} \otimes a^{(n_k)} \otimes I_{n_{k+1}}\otimes\cdots\otimes I_{n_{Q-1}},\] with \[a^{(n_k)} = \begin{bmatrix}0&1 & && \\ &0&\sqrt{2} && \\ &&\ddots&\ddots& \\ &&&0&\sqrt{n_k-1}\\ &&&&0\end{bmatrix}.\] The control Hamiltonian is given by \[H_c(t)=\sum_{k=0}^{Q-1} f^k(t)(a_k+a_k^\dagger),\] where \(f^k(t)\) is a real-valued control function. For the transmon model considered here, the jump operators in 2 are chosen to model relaxation and pure dephasing. More precisely, \[\mathcal{L}_{1k}=\frac{1}{\sqrt{T_1^k}}\,a_k, \qquad \mathcal{L}_{2k}=\frac{1}{\sqrt{T_2^k}}\,a_k^\dagger a_k,\] where \(T_1^k>0\) and \(T_2^k>0\) denote the characteristic relaxation and pure dephasing times for subsystem \(k\), respectively. Accordingly, the Lindbladian takes the form \[\mathcal{L}(\rho) = \sum_{k=0}^{Q-1}\sum_{l=1}^2 \mathcal{L}_{lk}\rho \mathcal{L}_{lk}^\dagger -\frac{1}{2} \left( \mathcal{L}_{lk}^\dagger \mathcal{L}_{lk}\rho + \rho \mathcal{L}_{lk}^\dagger \mathcal{L}_{lk} \right).\] The operators \(\mathcal{L}_{1k}\) and \(\mathcal{L}_{2k}\) represent energy dissipation and dephasing.

2.2 Classical Gregory Quadrature↩︎

In this paper, we will consider the approximation of the solution to the initial value problem 1 from time 0 to time \(T\) on an equidistant grid \[\label{eq:time95grid} t^n = n \Delta t, \;\;n = 0,\ldots,N_T, \;\;\Delta t = T/N_T.\tag{3}\] In our scheme, we introduce \(\Delta T := N_q^{(p)} \Delta t\) as the “big" interval where the numerical quadrature will be performed. Here \(N_q^{(p)}+1\) is the number of Gregory nodes. In order to achieve \(p\)-th order accuracy \(N_q^{(p)}\) must be at least \(2p-3\). This is the value we always choose for \(N_q^{(p)}\) in this paper.

We now review the basics of classical Gregory quadrature rule. For a scheme of order \(p \geq 2,\) we use \(N_q^{(p)} + 1\) nodes to compute an integral \(\int_{t^n}^{t^{n+N_q^{(p)}}} g(t) dt\) which is an integral of a smooth function \(g(t)\) over the time interval \([t^n,t^{n+N_q^{(p)}}]\) of length \(\Delta T\). We let \(\tau^j = t^n + j\Delta t\) with \(j = 0,\ldots,N_q^{(p)}\) and let \[g^j:=g(\tau^j),\qquad j=0,1,\dots,N_q^{(p)}.\] The classical Gregory rule may be viewed as an endpoint-corrected composite trapezoidal rule, only the weights near the two endpoints are modified, while the interior weights are equal to one. The general definition of quadrature weights is provided in [15]. In this paper, we use \[\label{eq:Gregory-weighted} Q_{N_q^{(p)}}^{\mathrm{G}}[g] := \Delta t \sum_{j=0}^{N_q^{(p)}} \omega_j^{(p)} g^j,\tag{4}\] where the weights \(\omega_j^{(p)}, j=0, \ldots p-2\) are provided in Table 1. The weights for \(j=p-1, \ldots, 2p-3\) are the same but with the order reversed. Note that with the choice we make here, \(N_q^{(p)} := 2p-3\), there are no interior weights. For the classical Gregory weights, if \(p \ge 10,\) there will be negative weights associated with the Gregory quadrature, which will break the desired CP property. Therefore, in this paper, we only consider order below or equal to 9 (we note that even higher-order methods can straightforwardly obtained by using the weights proposed by Fornberg [15]).

3 CPTP Schemes via NPI↩︎

In this section, we review the numerical framework developed in [1] for constructing high-order completely positive trace-preserving schemes based on NPI.

Following [8], the Lindblad equation can be reformulated in the form \[\label{eqn:Lindblad95in95J95term} \cfrac{d\rho(t)}{dt} = J(t)\rho(t) + \rho(t)J^\dagger(t) + \sum_\alpha \mathcal{L}_\alpha \rho(t)\mathcal{L}_\alpha^\dagger ,\tag{5}\] where \[J(t) = -i H_{\mathrm{eff}}, \qquad H_{\mathrm{eff}} = H(t) + \cfrac{1}{2i}\sum_\alpha \mathcal{L}_\alpha^\dagger \mathcal{L}_\alpha .\] For convenience we introduce the notation \[\mathcal{ L}_L\rho(t) = \sum_\alpha \mathcal{L}_\alpha \rho(t)\mathcal{L}_\alpha^\dagger .\]

To illustrate the main idea from [1], we first consider the simplified equation without \(\mathcal{L}_{L}(\rho(t))\) (which is already on Kraus form) \[\frac{d\rho(t)}{dt} = J(t)\rho(t)+\rho(t)J^\dagger(t), \;\;\rho(t_0) = \rho_0.\] Since the right-hand side is not expressed in Kraus form, directly applying a standard ODE solver (for instance an explicit Runge–Kutta method) does not guarantee preservation of the CPTP structure. However, the density matrix \(\rho(t)\) is positive semidefinite and therefore admits a factorization \[\rho(t) = V(t)V^\dagger(t).\] Substituting this representation into the above equation shows that solving for \(\rho(t)\) is equivalent to solving the matrix differential equation \[\frac{dV(t)}{dt} = J(t)V(t), \quad V(t_0)=V_0,\label{eq:J40t41}\tag{6}\] where the initial factor \(V_0\) satisfies \[\rho_0 = V_0 V_0^\dagger .\] Now, the solution of 6 can be written in terms of the flow operator \(U(t,s)\) associated with the equation, namely \[V(t) = U(t,s)V(s).\] Consequently the density matrix evolves according to \[\label{eq:dflow} \rho(t) = U(t,s)\rho(s)U^\dagger(t,s).\tag{7}\]

Equation 7 in combination with treating the Kraus term in the full Lindblad equation as a forcing makes it natural to apply Duhamel’s principle. Then the solution of the complete equation 5 can be expressed as \[\label{eq:Duhamel95form} \rho(t) = U(t,\tau)\rho(\tau)U^\dagger(t,\tau) + \int_\tau^t U(t,s)\,\mathcal{L}_L\rho(s)\,U^\dagger(t,s)\,ds.\tag{8}\]

3.1 Nested Picard Iteration using Gaussian Quadrature↩︎

In [1] we discretize 8 using NPI. The basic idea is to gain one order of accuracy per iteration. Precisely, using the time grid 3 with timestep \(\Delta t\) a solution of \(p\)-th accuracy can be obtained by the following iteration \[\begin{gather} \label{eq:NPI95general} \rho^{n+1,k+1} = U^{(k+1)}(t^{n+1},t^n) \rho^n U^{(k+1)}(t^{n+1},t^n) ^\dagger \\ +\Delta t \sum_{j=0}^{N_q^{(k+1)}} \omega_j^{(k+1)} U^{(k)}(t^{n+1},t^{n+c^{(k+1)}_j}) \mathcal{L}_L \rho^{n+c^{(k+1)}_j,k} U^{(k)}(t^{n+1},t^{n+c^{(k+1)}_j}) ^\dagger, \end{gather}\tag{9}\] with \(k =0,1,\cdots, p-1\). The iteration is started by \(\rho^{n+1,0} = \rho^n\). Here, \(U^{(k)}(t,s)\) represents a \(k\)-th order approximation of the flow operator \(U(t,s)\), and \(\omega_j^{(k+1)}\) and \(c_j^{(k+1)}\) are the weights and nodes of the \((k+1)\)-th order Gaussian quadrature rule. Note here we have flexibility to choose various types of Gaussian quadrature (Gaussian, Gauss-Radau, Gauss-Lobatto) which will have a slight effect on the efficiency of the method.

As an example, below we show the details of a third order scheme, where the quadrature rule is taken as Gauss-Radau when \(k=3\) and midpoint rule for \(k=2.\) Figure 1 illustrates the details of the computations involved. The computation of \(\rho^{n+1,3}\) requires evaluating intermediate solutions \(\rho^{n+\frac{1}{3},1}\) and \(\rho^{n+\frac{2}{3},2}\) at the new nodes \(t^{n+\frac{1}{3}}\) and \(t^{n+\frac{2}{3}}\). This process involves six flow operators of different orders acting over different subintervals.

We can give an estimate on how the cost scales with order of accuracy. For Gaussian quadrature the number of quadrature nodes, \(N_q^{(p)}+1\), to obtain order \(p\) is approximately \(p/2\) (accounting for the option to include endpoints or not). Let \(M(p)\) be the number of Gaussian nodes and let \(F(p)\) be the number of flow-operator evaluations with the sum taken over iteration \(1\) to \(p\) (i.e. \(k = 0,\ldots,p-1\) in 9 ). Then \(M(p)\) and \(F(p)\) can be computed by the recurrence relations \[\begin{align} M(p+1) &=& \left(M(p)+1\right)\left(N_q^{(p+1)}+1\right) \approx \left(M(p)+1\right) \frac{p+1}{2} ,\\ F(p+1) &=& F(p)\bigl(N_q^{(p+1)}+1\bigr)+N_q^{(p+1)}+2 \approx F(p)\frac{p+1}{2}+\frac{p+3}{2}. \end{align}\] Here the approximate values were obtained by using \(N_q^{(p)}+1 \approx p/2\). We can see clearly the estimates for \(M, F\) scale as \(\frac{p!}{2^p}\) which is unfavorable for high order methods.

This motivates the new scheme in this paper by replacing Gaussian quadrature with Gregory quadrature, which retains high-order accuracy but enables the use of equispaced nodes for computational efficiency. This, in turn, allows for reuse of computations resulting in linear complexity for \(M\) and \(F\). The details are described in the next section.

4 Gregory Nested Picard Iteration Schemes↩︎

We now describe the Gregory quadrature based NPI scheme for the computation of a \(p\)th-order approximation. We consider discretization of the Duhamel formula \[\begin{gather} \label{eq:Duhamel95form95GREG} \rho(t^n+\Delta T) = U(t^n+\Delta T,t^n)\rho(t^n)U^\dagger(t^n+\Delta T,t^n) \\ + \int_{t^n}^{t^n+\Delta T} U(t^n+\Delta T,s)\,\mathcal{L}_L\rho(s)\,U^\dagger(t^n+\Delta T,s)\,ds. \end{gather}\tag{10}\] Assuming that the density matrix has been approximated at desired order of accuracy at times \(t^n,t^{n+1},\ldots,t^{n+N_q^{(p)}-1},\) we then compute the density matrix at time \(t^{n+N_q^{(p)}}\) by the following NPI \[\label{eq:NPI95general95new95alt} \begin{align} &\rho^{n+N_q^{(p)},k+1} = U^{(p)}(t^{n+N_q^{(p)}},t^n)\rho^n U^{(p)}(t^{n+N_q^{(p)}},t^n)^\dagger \\ & + \Delta t\sum_{j=0}^{N_q^{(p)}}\omega_j^{(p)} U^{(p)}(t^{n+N_q^{(p)}},t^{n+j})\mathcal{L}_L\rho^{n+j,k} U^{(p)} (t^{n+N_q^{(p)}},t^{n+j})^\dagger, \end{align}\tag{11}\] for \(k=0,1,\dots,p-1\). The iteration is started by using the approximation at the previous time step, that is, \[\rho^{n+N_q^{(p)},0}=\rho^{n+N_q^{(p)}-1}.\] Note that in the iteration \(\rho^{n+j,k} := \rho^{n+j}, j = 0,\ldots, N_q^{(p)}-1\) for all \(k\).

This scheme is different from the NPI scheme from [1] in multiple ways. First, the weights, \(\omega_j^{(p)},\) are now the Gregory quadrature weights of order \(p\) (which are positive up to order 9, see Table 1). Second, here the flow operators are always applied at the order of the scheme, while they were applied at successively increasing order in [1]. Third, and from a computational efficiency point of view this is perhaps the most important improvement, the formulation allows us to split the quadrature term into \[\begin{gather} \Delta t\omega_{N_q^{(p)}}^{(p)} \boxed{ \mathcal{L}_L\rho^{n+N_q^{(p)},k}}\\ + \Delta t\sum_{j=0}^{N_q^{(p)}-1}\omega_j^{(p)} U^{(p)}(t^{n+N_q^{(p)}},t^{n+j})\mathcal{L}_L\rho^{n+j, k} U^{(p)}(t^{n+N_q^{(p)}},t^{n+j})^\dagger. \end{gather}\] Here the key observation is that we let \(\rho^{n+j, k}=\rho^{n+j},\) therefore throughout the iteration on \(k,\) the boxed expression is the only term that needs to be updated, and it does not require the application of a flow operator. Introducing \(\tilde{\rho}^{n+N_q^{(p)}}\) as \[\begin{gather} \tilde{\rho}^{n+N_q^{(p)}} :=U^{(p)}(t^{n+N_q^{(p)}},t^n)\rho^n U^{(p)}(t^{n+N_q^{(p)}},t^n)^\dagger \\ +\Delta t\sum_{j=0}^{N_q^{(p)}-1}\omega_j^{(p)} U^{(p)}(t^{n+N_q^{(p)}},t^{n+j})\mathcal{L}_L\rho^{n+j} U^{(p)}(t^{n+N_q^{(p)}},t^{n+j})^\dagger \end{gather}\] the Gregory-NPI scheme can thus be written compactly as \[\rho^{n+N_q^{(p)},k+1} := \tilde{\rho }^{n+N_q^{(p)}}+ \Delta t\omega_{N_q^{(p)}}^{(p)} \mathcal{L}_L\rho^{n+N_q^{(p)},k}, \;\;k = 0,\ldots,p-1.\]

As just mentioned, \(\tilde{\rho}^{\,n+N_q^{(p)}}\) can be precomputed before the iteration starts, and that no flow operators are involved during the subsequent iteration stages. The algorithm for Gregory-NPI is described in Algorithm 2.

More generally, to obtain a solution of order \(p\), the proposed scheme requires only \((N_q^{(p)}+1)\) flow-operator evaluations. This is a substantial improvement over the Gaussian-quadrature-based NPI schemes.

We note that the Gregory NPI method belongs to the general linear methods (GLM) [16]. The first few steps \[\rho^1,\dots,\rho^{N_q^{(p)}-1}\] need to be initialized. In practice, one often uses either a higher-order one-step method or a finer time step size in this startup phase to ensure that the starting error does not contaminate the overall convergence order. Here we generate these starting values on an auxiliary fine grid with step size \(\Delta t/m\), where \(m=O(\Delta t^{-(p-3)/3})\) so that \((\Delta t/m)^3=O(\Delta t^p)\). The second-order scheme is first used to compute the approximations at the first two fine-grid nodes, not counting the initial value. The order is then increased successively: for \(q=3,\ldots,p-1\), the order-\(q\) scheme is used to compute the next two values, and the order-\(p\) scheme is applied from the fine-grid \((2p-3)\)-th node onward. Other approaches are of course also possible.

The approximation of the flow operator \(U(t^{n+k},t^n)\) associated with \[\frac{dV(t)}{dt}=J(t)V(t), \qquad V(t^n)=V_n.\] have been discussed in prior work [1]. In our numerical simulations, explicit approximations of the flow operator are constructed using Runge–Kutta methods. For implicit approximations, we use the trapezoidal rule for the second-order scheme and the method proposed in [17] for the fourth-order scheme. Other high-order implicit time schemes can be used but we do not consider those in this paper.

5 Stability Analysis↩︎

Theoretical results of GLM [16] have been discussed in the literature. The proposed Gregory NPI methods are not designed to be a GLM for general ODEs, but they are tailored for the Lindblad equation. Therefore, in this paper, we focus the analysis specifically on the Lindblad equation.

Although trace renormalization is a nonlinear operation that keeps the density matrix bounded, it may compromise the accuracy of the scheme. This effect is illustrated in the Qudit-cavity example of Section 7.2, where the numerical solution remains bounded but is not accurate. It is therefore necessary to study the stability of 11 without trace renormalization, which is what we will do in this section.

Similar to [1], we consider a single-qubit test equation with the Hamiltonian and the Lindblad operators of the model are given by (note that we use \(\Omega\) instead of \(\omega\) so as to not confuse it with the quadrature weights) \[H = \Omega \begin{bmatrix}0&0\\0&1\end{bmatrix}, \mathcal{L}_1= \frac{1}{\sqrt{T_1}}\begin{bmatrix}0&1\\0&0\end{bmatrix}, \;\;\mathcal{L}_2 = \frac{1}{\sqrt{T_2}}\begin{bmatrix}0&0\\0&1\end{bmatrix},\] This will give rise to \[J = -iH -\frac{1}{2} \mathcal{L}_1^\dagger \mathcal{L}_1 -\frac{1}{2} \mathcal{L}_2^\dagger \mathcal{L}_2 = \begin{bmatrix}0&0\\0&-i\Omega -\frac{1}{2T_1}-\frac{1}{2T_2}\end{bmatrix}.\]

Denote \(\mathcal{W} := \Delta t \Omega\), \(\mathcal{T}_1 := \frac{\Delta t}{T_1}\) and \(\mathcal{T}_2 := \frac{\Delta t}{T_2}\). All approximations of the flow operator \(U(t^{n+k},t^n)\) considered in the following subsections take the form

\[\begin{bmatrix}1&0\\0&\varsigma^{(p,j)}(\mathcal{W},\mathcal{T}_1, \mathcal{T}_2, k)\end{bmatrix},\] where \(k\) denotes the number of time steps from \(t^n\) to \(t^{n+k}\), \(p\) denotes the order of approximation, and \(j\in\{0,1\}\) indicates whether the approximation is explicit or implicit. More precisely, \[j= \begin{cases} 0, & \text{if the approximation is explicit},\\ 1, & \text{if the approximation is implicit}. \end{cases}\]

The quantity \(\varsigma^{(p,j)}(\mathcal{W},\mathcal{T}_1,\mathcal{T}_2,k)\) depends on the particular approximation of the flow. When there is no ambiguity, \(\varsigma^{(p,j)}(k)\) means \(\varsigma^{(p,j)}(\mathcal{W},\mathcal{T}_1, \mathcal{T}_2, k)\). We also use \(\varsigma^{(p)}(k)\) when whether the flow operator is implicit or explicit is not specified.

We also denote \(A(k):= -\frac{k}{2}\frac{\Delta t}{T_1} - \frac{k}{2}\frac{\Delta t}{T_2}-ik\Omega\Delta t = -\frac{k}{2}\mathcal{T}_1-\frac{k}{2}\mathcal{T}_2-ik\mathcal{W}\), which will be used in the stability analysis. Next, we derive matrix representations for the linear operators used in our numerical scheme, which will be needed in the stability analysis. Using columnwise vectorization and the identity \[\operatorname{vec}(AXB)=(B^T\otimes A)\operatorname{vec}(X),\] the map \(\rho\mapsto U\rho U^\dagger\) is represented by \[\operatorname{vec}(U\rho U^\dagger) = (\overline{U}\otimes U)\operatorname{vec}(\rho),\] where \(\overline{U}\) denotes the entrywise complex conjugate of \(U\). It follows that the matrices \[\begin{align} P_k^{(p)}&=\begin{bmatrix}1&0&0&0\\0&\overline{\varsigma^{(p)}(k)}&0&0\\0&0&\varsigma^{(p)}(k)&0\\0&0&0&|\varsigma^{(p)}(k)|^2\end{bmatrix},\\[1mm] Q_{U^{(p)}L,k} &= \begin{bmatrix}0&0&0&\frac{1}{T_1}\\0&0&0&0\\0&0&0&0\\0&0&0&\frac{|\varsigma^{(p)}(k)|^2}{T_2}\end{bmatrix},\\[1mm] Q_L &= \begin{bmatrix}0&0&0&\frac{1}{T_1}\\0&0&0&0\\0&0&0&0\\0&0&0&\frac{1}{T_2}\end{bmatrix}, \end{align}\] represent the operators \[\begin{align} & \operatorname{vec}(\rho) \mapsto \operatorname{vec}\left(U^{(p)}(t^{n+k},t^n)\rho \bigl(U^{(p)}(t^{n+k},t^n)\bigr)^\dagger\right),\\ &\operatorname{vec}(\rho) \mapsto \operatorname{vec}\left(U^{(p)}(t^{n+k},t^n)\mathcal{L}_L\rho \bigl(U^{(p)}(t^{n+k},t^n)\bigr)^\dagger\right),\\ &\operatorname{vec}(\rho) \mapsto \operatorname{vec}\left(\mathcal{L}_L\rho\right), \end{align}\] respectively. These three operators form the main components of our schemes. Their matrix representations indicate that \(\varsigma^{(p)}(k)\), \(T_1\), and \(T_2\) play essential roles in determining stability.

We begin with the stability analysis of the second-order scheme, which is essentially a one-step method. Similar analysis has been carried out in [1]. We then proceed to the higher-order schemes, whose stability must be analyzed in the framework of GLM.

5.1 Second-Order Scheme↩︎

The second-order approximation of the solution is defined by \[\begin{align} \rho^{n+1,0} & = \rho^n \nonumber\\ \rho^{n+1,1} &= U^{(2)}(t^{n+1},t^n) \rho^n U^{(2)}(t^{n+1},t^n) ^\dagger +\nonumber \\ &\qquad {\Delta t}\Big(\omega^{(2)}_0 U^{(2)}(t^{n+1},t^{n}) \mathcal{L}_L\rho^n U^{(2)}(t^{n+1},t^{n}) ^\dagger +\omega^{(2)}_1 \mathcal{L}_L {\rho^{n+1,0}} \Big),\nonumber\\ \rho^{n+1,2} &= U^{(2)}(t^{n+1},t^n) \rho^n U^{(2)}(t^{n+1},t^n) ^\dagger +\nonumber\\ &\qquad {\Delta t}\Big(\omega^{(2)}_0 U^{(2)}(t^{n+1},t^{n}) \mathcal{L}_L\rho^n U^{(2)}(t^{n+1},t^{n}) ^\dagger +\omega^{(2)}_1 \mathcal{L}_L {\rho^{n+1,1}} \Big),\label{eqn:32NPI95second-order} \end{align}\tag{12}\] where \(\omega^{(2)}_0 = \omega^{(2)}_1 = \frac{1}{2}\). For the second-order approximation of the flow operator, we apply the explicit midpoint method in the explicit case and the trapezoidal rule in the implicit case, i.e. \[U^{(2)}(t^{n+k},t^n)= \begin{cases} I+k\Delta t\,J+\dfrac{(k\Delta t)^2}{2}J^2, & \text{if } j=0, \qquad{\rm i.e., explicit}, \\[1.2ex] \left(I-\dfrac{k\Delta t}{2}J\right)^{-1}\left(I+\dfrac{k\Delta t}{2}J\right), & \text{if } j=1, \qquad{\rm i.e., implicit}. \end{cases}\] With this flow operator we get \[\varsigma^{(2,j)}(k)= \begin{cases} 1+A(k) +\dfrac{1}{2}A(k)^2, & \text{if } j=0, \\[2ex] \dfrac{ 1+\dfrac{1}{2}A(k) }{ 1-\dfrac{1}{2}A(k) }, & \text{if } j=1. \end{cases}\]

To analyze the stability of the scheme, we first rewrite the matrix equation in vectorized form. This yields an amplification matrix representation of the method, allowing the stability analysis to be reduced to a spectral study of the corresponding matrix, in particular through its eigenvalues and spectral radius. Throughout the following discussion, boldface symbols denote vectorized matrices. In particular, \(\boldsymbol{\rho}\) denotes the vectorized form of \(\rho\). Vectorizing 12 yields \[\begin{align} \boldsymbol{ \rho}^{n+1,1} & =P_1^{(2)} \boldsymbol{\rho}^n + \Delta t\omega_0^{(2)} Q_{U^{(2)}L,1} \boldsymbol{\rho}^n + \Delta t\omega_1^{(2)} Q_L {\boldsymbol{\rho}^{n+1,0}} , \\ \boldsymbol{\rho}^{n+1,2} & = P_1^{(2)} \boldsymbol{\rho}^n + \Delta t\omega_0^{(2)} Q_{U^{(2)}L,1} \boldsymbol{\rho}^n + \Delta t\omega_1^{(2)} Q_L {\boldsymbol{\rho}^{n+1,1}}. \end{align}\]

Substituting the expression for \(\boldsymbol{\rho}^{n+1,1}\) into the equation for \(\boldsymbol{\rho}^{n+1,2}\) yields \[\begin{gather} \boldsymbol{\rho}^{n+1,2} = \Bigl( P_1^{(2)} + \Delta t\omega_0^{(2)} Q_{U^{(2)}L,1} + {\Delta t}\omega_1^{(2)} Q_L P_1^{(2)} +\\(\Delta t)^2 \omega_0^{(2)} \omega_1^{(2)} Q_L Q_{U^{(2)}L,1}+ (\Delta t)^2(\omega_1^{(2)} )^2( Q_L )^2 \Bigr)\boldsymbol{\rho}^n. \end{gather}\] The amplification matrix of the second-order approximation is \[\begin{align} & P_1^{(2)} + \Delta t\omega_0^{(2)} Q_{U^{(2)}L,1} + {\Delta t}\omega_1^{(2)} Q_L P_1^{(2)} +(\Delta t)^2 \omega_0^{(2)} \omega_1^{(2)} Q_L Q_{U^{(2)}L,1}\\ &+ (\Delta t)^2(\omega_1^{(2)} )^2( Q_L )^2 =\\ & \begin{bmatrix}1&0&0& | \varsigma^{(2)}(1)|^2 \left(\frac{\Delta t}{T_1}\omega_1^{(2)} +\omega_0^{(2)}\omega_1^{(2)} \frac{(\Delta t)^2}{T_1T_2} \right) +\frac{\Delta t}{T_1}(\omega_0^{(2)} + \frac{\omega_1^{(2)} \Delta t}{T_2}) \\ 0&\overline{\varsigma^{(2)}(1)}&0&0 \\ 0&0&\varsigma^{(2)}(1)&0 \\ 0&0&0&|\varsigma^{(2)}(1)|^2 \left( 1+ \frac{\Delta t}{T_2} + \omega_0^{(2)}\omega_1^{(2)} (\frac{\Delta t}{T_2})^2 \right) + (\frac{\omega_1^{(2)}\Delta t}{T_2})^2\end{bmatrix}. \end{align}\] So the stability condition for the second-order scheme is \[\left| |\varsigma^{(2)}(1)|^2 \left( 1+ \frac{\Delta t}{T_2} + \omega_0^{(2)}\omega_1^{(2)} (\frac{\Delta t}{T_2})^2 \right)+ (\frac{\omega_1^{(2)}\Delta t}{T_2})^2 \right| \le 1.\] Note that, as \(\Delta t/T_2 \to 0\), the stability condition approaches \[\left| \varsigma^{(2)}(1)\right|\le 1,\] that is, the stability condition associated with the second-order flow solver. For illustration, Figure 3 shows the stability regions of the second-order schemes for three representative relations between the pure dephasing time and the decay time, namely, \[T_2=2T_1,\qquad T_2=T_1,\qquad T_2=\frac{T_1}{2}.\] In all three cases, the scheme with the implicit flow solver has a much larger stability region than the corresponding scheme with the explicit flow solver.

5.2 Fourth-Order Scheme↩︎

The fourth-order scheme is based on the fourth-order Gregory quadrature rule, which requires at least six quadrature nodes. The fourth-order approximation is defined by \[\begin{align} \rho^{n+5,0} & = \rho^{n+4} \nonumber\\ \rho^{n+5,k} &= U^{(4)}(t^{n+5},t^n) \rho^n U^{(4)}(t^{n+5},t^n) ^\dagger +\\ &\Delta t\Big( \omega^{(4)}_0 U^{(4)}(t^{n+5},t^n) \mathcal{L}_L \rho^n U^{(4)}(t^{n+5},t^n) ^\dagger+\nonumber\\ & \omega^{(4)}_1 U^{(4)}(t^{n+5},t^{n+1}) \mathcal{L}_L \rho^{n+1} U^{(4)}(t^{n+5},t^{n+1}) ^\dagger+\nonumber\\ & \omega^{(4)}_2 U^{(4)}(t^{n+5},t^{n+2}) \mathcal{L}_L \rho^{n+2} U^{(4)}(t^{n+5},t^{n+2}) ^\dagger+\nonumber\\ &\omega^{(4)}_3 U^{(4)}(t^{n+5},t^{n+3}) \mathcal{L}_L \rho^{n+3} U^{(4)}(t^{n+5},t^{n+3}) ^\dagger+\nonumber\\ & \omega^{(4)}_4U^{(4)}(t^{n+5},t^{n+4}) \mathcal{L}_L \rho^{n+4} U^{(4)}(t^{n+5},t^{n+4}) ^\dagger \Big)\nonumber\\ & + \Delta t \omega^{(4)}_5\mathcal{L}_L { \rho^{n+5,k-1}}, \qquad \label{eqn:32approximation95order954} \end{align}\tag{13}\] for \(k=1,2,3,4.\) Here, the Gregory weights are \(\omega^{(4)}_0 =\omega^{(4)}_5 = \frac{3}{8}, \omega^{(4)}_1 =\omega^{(4)}_4 = \frac{7}{6}, \omega^{(4)}_2 =\omega^{(4)}_3 = \frac{23}{24}\). For the fourth-order approximation of the flow operator, we can apply the classical fourth-order Runge-Kutta method in the explicit case, namely, \[U^{(4)}(t^{n+k},t^n)= I + k\Delta t J + \frac{(k\Delta t)^2}{2}J^2+\frac{(k\Delta t)^3}{6}J^3+\frac{(k\Delta t)^4}{24}J^4,\] or the fourth-order implicit scheme from [17] in the implicit case, namely, \[U^{(4)}(t^{n+k},t^n)= (I + i\frac{k\Delta t}{4}\overline{d} J)^{-1}(I + i\frac{k\Delta t}{4} d J)\, (I - i\frac{k\Delta t}{4}d J)^{-1}(I - i\frac{k\Delta t}{4}\overline{d} J),\] with \(d=\frac{1}{\sqrt{3}}-i\). It follows that \[\varsigma^{(4,j)}(k)=\begin{cases} 1+ A(k)+\frac{1}{2}A(k)^2 +\frac{1}{6}A(k)^3 +\frac{1}{24}A(k)^4,& \text{if } j=0,\\ \dfrac{( 1+i\frac{d}{4}A(k))(1-i\frac{\bar{d}}{4}A(k)) }{ (1+i\frac{\bar{d}}{4}A(k)) ( 1-i\frac{d}{4}A(k))}, & \text{if } j=1. \end{cases}\] We present the following stability result whose proof can be found in the Appendix.

Note that, for sufficiently large dephasing time, or equivalently, for sufficiently small \(\frac{\Delta t}{T_2}\), the stability region approaches the one defined by \[\left|\varsigma^{(4)}(5)\right|\le1.\]

In other words, when the dephasing time is large, the stability region is determined primarily by the stability of the flow operator. This phenomenon is also observed for the second-order method. In Figure 4, we plot the stability regions of the second- and fourth-order schemes for time step sizes \(\Delta t=0.2, 0.5, 0.8, 1.0\), under the constraint \(T_2=2T_1\). For current quantum devices, the relaxation and dephasing times are typically large enough that \(1/(2T_1)\) remains small, and practical choices of \(\Delta t\) therefore fall within the stability region.

5.3 General \(p\)-th Order Scheme↩︎

We can generalize the stability results for higher order.

The proof is similar to the proof of Theorem 1 and is skipped for brevity.

Here, similar to previous cases, when \(\frac{\Delta t}{T_2}\rightarrow 0\), i.e., when the pure dephasing time is sufficiently large, the stability condition is governed almost entirely by the approximation of the flow operator. Consequently, for higher-order schemes in which the flow operator is approximated by a Runge–Kutta method, the stability of the overall scheme is essentially inherited from that of the underlying Runge–Kutta method.

6 Low-Rank Truncation↩︎

We now describe the low-rank truncation procedure for scheme 11 . We represent the density matrix in factored form as \[\rho(t)=V(t)V(t)^\dagger,\] where \(\rho(t)\in\mathbb{C}^{d\times d}\) and \(V(t)\in\mathbb{C}^{d\times r(t)}\) with \(r(t)\ll d\). Thus, instead of storing and evolving \(\rho(t)\) directly, we work with its low-rank factor \(V(t)\).

At each update, the truncation is performed in two stages. The first stage is applied after the offline step. We write the intermediate approximation \(\widetilde{\rho}^{\,n+N_q^{(p)}}\) as \[\widetilde{\rho}^{\,n+N_q^{(p)}} = \widetilde{V}^{\,n+N_q^{(p)}} \bigl(\widetilde{V}^{\,n+N_q^{(p)}}\bigr)^\dagger,\] where \[\widetilde{V}^{\,n+N_q^{(p)}} = \Bigl[ U^{(p)}(t^{n+N_q^{(p)}},t^n)V^n ~~ \cdots ~~ \sqrt{\omega_j^{(p)}\Delta t}\, U^{(p)}(t^{n+N_q^{(p)}},t^{n+j})\mathcal{L}_{\alpha}V^{n+j} \quad \cdots \Bigr].\] Here, \(j=0,\dots,N_q^{(p)}-1.\) The second stage is applied during the nested Picard iteration. We set \[\rho^{n+N_q^{(p)},0} = V^{n+N_q^{(p)}-1}\bigl(V^{n+N_q^{(p)}-1}\bigr)^\dagger,\] and, for \(k=1,\dots,p\), \[\rho^{n+N_q^{(p)},k} = V^{n+N_q^{(p)},k}\bigl(V^{n+N_q^{(p)},k}\bigr)^\dagger,\] where \[V^{n+N_q^{(p)},k} = \Bigl[ \widetilde{V}^{\,n+N_q^{(p)}} \qquad \sqrt{\omega_{N_q}^{(p)}\Delta t}\, \mathcal{L}_{\alpha}V^{n+N_q^{(p)},k-1} \Bigr].\] To maintain a low-rank representation, we truncate \[\widetilde{V}^{\,n+N_q^{(p)}},\qquad V^{n+N_q^{(p)},k},\quad k=1,\dots,p,\] by means of a truncated singular value decomposition. Following [1], we choose the truncation tolerance proportional to the local truncation error. For a method of order \(p\), we take \[\varepsilon_{\mathrm{tr}}=(\kappa\Delta t)^{p+1},\] where \(\kappa>0\) is a prescribed constant.

When the number of columns in \(V^{n+N_q,k}\) becomes large, a direct truncated SVD may be computationally expensive. In that case, we use the same hierarchical truncation procedure as in [1], namely, we group the columns and apply truncated SVD successively to each group. Since each truncation step is performed at the level of the factors \(\widetilde{V}^{\,n+N_q}\) and \(V^{n+N_q,k}\), the resulting density matrix after truncation remains completely positive.

After advancing to time \(t^{n+N_q^{(p)}}\), we apply the trace renormalization \[\rho^{n+N_q^{(p)}} \leftarrow \frac{\rho^{n+N_q^{(p)}}}{\operatorname{Tr}(\rho^{n+N_q^{(p)}})}.\] We note that the dominant cost of the truncation procedure arises in the offline stage, while the online stage is relatively inexpensive. Compared with the scheme in [1], the truncation procedure in the present work is computationally cheaper.

7 Numerical Experiments↩︎

Direct simulation in the laboratory frame becomes expensive when \(H(t)\) contains highly oscillatory terms, since these induce stiffness and force very small time steps. To reduce this difficulty, we work in a rotating frame and apply the rotating-wave approximation, thereby eliminating much of the rapidly oscillatory behavior.

Given rotation frequencies \(\omega_k^r\), the transformed drift and control Hamiltonians are \[\begin{align} \widetilde{H}_d(t) =& \sum_{k=0}^{Q-1} \left( (\omega_k-\omega_k^r)a_k^\dagger a_k -\frac{\xi_k}{2}\, a_k^\dagger a_k^\dagger a_k a_k -\sum_{l>k}\xi_{kl} a_k^\dagger a_k a_l^\dagger a_l \right) \\ &\quad +\sum_{k=0}^{Q-1}\sum_{l>k} J_{kl} \left( \cos(\eta_{kl}t)(a_k^\dagger a_l+a_k a_l^\dagger) +i\sin(\eta_{kl}t)(a_k^\dagger a_l-a_k a_l^\dagger) \right), \\ \widetilde{H}_c(t) =& \sum_{k=0}^{Q-1} \left( p^k(\vec{\alpha}^k,t)(a_k+a_k^\dagger) +i q^k(\vec{\alpha}^k,t)(a_k-a_k^\dagger) \right), \end{align}\] where \[\eta_{kl}:=\omega_k^r-\omega_l^r.\] The laboratory-frame control can then be approximated by \[f^k(t)\approx 2p^k(t)\cos(\omega_k^r t)-2q^k(t)\sin(\omega_k^r t).\]

The relaxation and dephasing parameters used in the following experiments are chosen to be consistent with experimentally reported coherence scales in superconducting circuits. Recent transmon devices report relaxation and coherence times exceeding \(100\,\mu\mathrm{s}\), with high-coherence platforms reaching several hundred microseconds, while cavity or resonator memories can exhibit lifetimes from the millisecond regime to several tens of milliseconds. Thus, the values considered here are intended to be representative of experimentally relevant regimes rather than tied to a single device realization.

In the following subsections, we use standard bra–ket notation, where \[\ket{j_1j_2\cdots j_Q}:=\ket{j_1}\otimes\ket{j_2}\otimes\cdots\otimes\ket{j_Q}\] denotes a tensor-product basis state, and \(\bra{\psi}=\ket{\psi}^\dagger\).

7.1 Two-Qubit System↩︎

In the first experiment, we investigate the accuracy of several numerical methods for a composite quantum system consisting of two identical qubits with the same transition frequency, that is, \(\omega_0=\omega_1\). In the rotating frame, we choose the rotation frequencies to coincide with the physical frequencies, namely, \(\omega_k^r=\omega_k\) for \(k=0,1\). With this choice, the terms \(\sum_k \omega_k a_k^\dagger a_k\) are eliminated from the system Hamiltonian. The dipole-dipole interaction is characterized by a single cross-resonance coefficient \(J_{01}=0.2\) (in nondimensionalized units), so that the Hamiltonian reduces to \[H = J_{01}(a_0^\dagger a_1 + a_0 a_1^\dagger).\] We consider Lindblad operators associated with decay and dephasing processes, with uniform rates \[\mu = \frac{1}{T_1^0}=\frac{1}{T_1^1}=\frac{1}{50}, \qquad \gamma_1=\frac{1}{T_2^0}=\frac{1}{50}, \qquad \gamma_2=\frac{1}{T_2^1}=\frac{1}{50}.\] Define \[a=\frac{\gamma_1+\gamma_2}{4}, \qquad b=\sqrt{4J_{01}^2-a^2}.\]

For the initial condition \(\rho(0)=|10\rangle\langle 10|\), the exact solution at time \(t\) is given by \[\rho(t) = \begin{pmatrix} \rho_{11}(t)&0&0&0\\ 0&\rho_{22}(t)&\rho_{23}(t)&0\\ 0&\rho_{32}(t)&\rho_{33}(t)&0\\ 0&0&0&0 \end{pmatrix},\] where \[\begin{align} \rho_{11}(t)&=1-e^{-\mu t},\\ \rho_{22}(t) &= \frac{1}{2}\left[ e^{-\mu t} - e^{-(\mu+a)t} \left( \cos(b t)+\frac{a}{b}\sin(b t) \right) \right],\\ \rho_{33}(t) &= \frac{1}{2}\left[ e^{-\mu t} + e^{-(\mu+a)t} \left( \cos(b t)+\frac{a}{b}\sin(b t) \right) \right],\\ \rho_{23}(t) &= -\,i\,\frac{J_{01}}{b}\,e^{-(\mu+a)t}\sin(b t), \qquad \rho_{32}(t)=\,i\,\frac{J_{01}}{b}\,e^{-(\mu+a)t}\sin(b t). \end{align}\] We simulate the density matrix up to the final time \(T=6\) (in nondimensionalized units) using methods of orders two through nine, and measure the error of the density matrix at the final time in the Frobenius norm. Table ¿tbl:tab:2Q1? reports the errors for the second-, third-, and fourth-order methods. Here, “explicit” and “implicit” indicate whether the flow operator is approximated by an explicit or an implicit method. The observed convergence rates agree with the expected orders, except for the implicit third-order scheme, which exhibits fourth-order convergence. This is likely due to the use of a fourth-order implicit flow operator approximation in the third-order implicit scheme. The same phenomenon was observed in [1].

Figure 5 shows the error versus the time step size for methods of orders two through nine, where the flow operator is approximated by the Runge–Kutta method associated with each scheme. The results again confirm that the observed convergence rates are consistent with the expected orders.

7.2 Qudit-Resonator System with Control↩︎

In this section, we consider a composite quantum system consisting of a three-level qudit and a twenty-level resonator. The device model is standard and was also employed in [14] and our earlier work [1]. We adopt the same parameter values as in [1]; these are reported in Table 2. The corresponding control functions are given by \[p(t) = \frac{A}{2} \left(1 + \tanh (\delta ( t - \tau)) \right), \;\;q(t) = 0.\] The values of \(A, \delta\), and \(\tau\) are chosen as 10 MHz, −0.05 GHz and 2000 ns for the qudit and 15 MHz, −0.1 GHz and 200 ns for the resonator. We simulate the problem up to the final time \(T=\SI{2500}{\nano\second}\). We first study the convergence of the proposed schemes, then demonstrate the low-rank method, and finally present an example showing that trace renormalization, while preventing blow-up, may alter the accuracy of the solution. To obtain the reference solution, we use a large number of time steps for fourth-order scheme with implicit flow operator.

To study the convergence behavior, we consider the third- and fourth-order schemes with both explicit and implicit flow operators. The error is measured in the Frobenius norm of the density matrix at the final time, and Figure 6 plots the error versus the time step size. The observed convergence rates are consistent with the theoretical orders, except for the fourth-order scheme with an explicit flow operator, whose rate is slightly below the expected one. This is likely due to the fact that the error is already very small and hence influenced by machine precision. In addition, for a prescribed error level, the schemes with implicit flow operators require substantially fewer time steps than the corresponding explicit schemes. We also note that, for each method, the largest time step size used in the figure is close to the empirical maximum value for which trace renormalization does not noticeably deteriorate the accuracy. Overall, Figure 6 shows that the schemes with implicit flow operators offer greater flexibility in the choice of time step size.

Figure 7 displays the rank history for the fourth-order scheme. Here we take \(\kappa=0.33\) for the implicit-flow schemes and \(\kappa=0.53\) for the explicit-flow schemes. These parameters are chosen empirically so that the truncation error is comparable to the discretization error before truncation. For both types of flow operators, the final ranks remain similar across different time step sizes. While finer time steps may lead to temporarily larger intermediate ranks, truncation reduces them effectively, and the final rank increases only slightly. Together with Figure 6, this indicates that the scheme with the implicit flow operator is more efficient.

Finally, we investigate the effect of renormalization on the numerical solution. Figure 8 compares the population of the qudit obtained with a very large number of time steps (left), which serves as a reference solution, and that produced by the fourth-order explicit scheme with \(10000\) time steps (right). The right figure shows that renormalization keeps the solution stable, but does not guarantee sufficient accuracy. This is one of the main reasons for studying the stability of the scheme without renormalization in Section 5.

7.3 Two-Qudit-Resonator System with a CNOT Gate↩︎

In this experiment, we consider a composite quantum system consisting of two qudits, each with four energy levels, and a resonator with ten energy levels, in both closed- and open-system settings. We adopt the same physical parameters and the same CNOT control pulse as in [2]. The parameter values are listed in Table [table:parameters952qudit95resonator], and the control pulse is shown in Figure 9. In the closed-system case, the dynamics are governed by the von Neumann equation, that is, without the Lindbladian term. Our primary goal is to assess whether the control pulse from [2] remains effective for the von Neumann equation and whether its performance is preserved when higher energy levels are included. We also examine how this pulse, originally designed for a closed system, performs in the corresponding open-system setting.

We first examine the performance of the CNOT gate in the closed-system setting governed by the von Neumann equation. We consider four pure initial states, \[\ket{000},\quad \ket{010},\quad \ket{100},\quad \ket{110}.\] Figure 10 shows that the gate flips qudit 1 between the ground and excited states when qudit 2 is in the excited state, while leaving qudit 1 unchanged when qudit 2 is in the ground state. This behavior is consistent with the intended action of the CNOT gate.

We next increase the number of resonator energy levels from 10 to 20 while keeping each qudit at four levels. We also consider a larger configuration with 6 energy levels for each qudit and 20 energy levels for the resonator. In all cases, the control pulse exhibits the same qualitative behavior. For brevity, we omit the corresponding figures, as they are qualitatively similar to Figure 10.

Next, we test the same CNOT control pulse in the open-system setting. The relaxation and dephasing times for the two qudits and the resonator are chosen as \[T_{1}^{1}=T_{1}^{2}=40~\si{\micro\second},\qquad T_{1}^{R}=200~\si{\micro\second},\qquad T_{2}^{1}=T_{2}^{2}=20~\si{\micro\second},\qquad T_{2}^{R}=500~\si{\micro\second}.\] These values are experimentally plausible. The initial density matrix is taken to be \[\rho_0=\ket{\psi}\bra{\psi}, \qquad \ket{\psi} = \frac{1}{2}\bigl(\ket{000}+\ket{010}+\ket{100}+\ket{110}\bigr).\] We compute the solution using the fourth-order scheme with implicit flow operator and \(\Delta t=0.01\). The corresponding population dynamics are shown in Figure 11. Let \(P_{ijk}(t)\) denote the population of the state \(\ket{ijk}\) at time \(t\). At the final time, the numerical results give \[|P_{000}(550)-P_{000}(0)| = 2.79\times 10^{-4},\qquad |P_{010}(550)-P_{110}(0)| = 5.57\times 10^{-3},\] \[|P_{110}(550)-P_{010}(0)| = 5.54\times 10^{-3},\qquad |P_{100}(550)-P_{100}(0)| = 1.37\times 10^{-4}.\]

These results indicate that, for the above parameter regime, the control pulse remains effective in realizing the CNOT gate in the open-system setting.

8 Conclusion↩︎

In this work, we developed novel high-order low-rank structure-preserving numerical schemes for the Lindblad equation. The proposed methods build on NPI framework with two key ideas. First, Gaussian quadrature is replaced by Gregory quadrature, which allows the use of equidistant nodes while maintaining high-order accuracy. Second, the approximation of the flow operator is fixed at the highest available order within the construction, which further improves efficiency in practical computations. As a result, the new schemes retain the main advantages of the previous framework, including high-order accuracy, low-rank structure, and preservation of the key physical properties of the density matrix, while reducing the overall computational cost.

We also investigated the stability of the resulting schemes without trace renormalization. For a representative single-qubit model problem, the analysis clarifies the role of the flow operator approximation and shows how the stability behavior depends on the dissipative parameters and the time step size. In particular, for sufficiently large dephasing times, the stability region is determined primarily by the stability of the flow operator. This provides useful guidance for the practical implementation of the proposed methods. Low-rank truncation is incorporated into the framework. The numerical experiments demonstrate that the proposed methods achieve the expected convergence orders and perform well on several representative problems, including a two-qubit problem with known exact solution, a controlled qudit-resonator system, and CNOT gate simulations for two qudits coupled with a resonator in both closed- and open-system settings.

Appendix↩︎

In the Appendix, we establish the proof of the stability results for the fourth-order scheme. The analysis is carried out by first vectorizing the matrices and then reformulating the multistep scheme as an equivalent one-step method. The stability is then characterized by the amplification matrix of the resulting one-step system.

9 Proof of Theorem 1↩︎

By vectorizing the matrices \(\rho^{n+5,k}\), we can rewrite the fourth-order approximation in 13 as follows:

\[\begin{align} \boldsymbol{\rho}^{n+5,k} &= P_5^{(4)}\boldsymbol{\rho}^n + \Delta t\Big( \omega^{(4)}_0 Q_{U^{(4)}L,5}\boldsymbol{\rho}^n + \omega^{(4)}_1 Q_{U^{(4)}L,4}\boldsymbol{\rho}^{n+1} + \omega^{(4)}_2 Q_{U^{(4)}L,3}\boldsymbol{\rho}^{n+2} \\ &\quad + \omega^{(4)}_3 Q_{U^{(4)}L,2}\boldsymbol{\rho}^{n+3} + \omega^{(4)}_4 Q_{U^{(4)}L,1}\boldsymbol{\rho}^{n+4} \Big) + \Delta t\,\omega^{(4)}_5 Q_L \boldsymbol{\rho}^{n+5,k-1}, \end{align}\] for \(k=1,2,3,4\). Introducing the augmented vector \[\begin{bmatrix}\boldsymbol{\rho}^{n+5}\\ \boldsymbol{\rho}^{n+4}\\ \boldsymbol{\rho}^{n+3}\\ \boldsymbol{\rho}^{n+2}\\ \boldsymbol{\rho}^{n+1}\end{bmatrix} := \begin{bmatrix}\boldsymbol{\rho}^{n+5,4}\\ \boldsymbol{\rho}^{n+4}\\ \boldsymbol{\rho}^{n+3}\\ \boldsymbol{\rho}^{n+2}\\ \boldsymbol{\rho}^{n+1}\end{bmatrix},\] the scheme can be formulated as the one-step system \[\begin{bmatrix}\boldsymbol{\rho}^{n+5}\\\boldsymbol{\rho}^{n+4}\\\boldsymbol{\rho}^{n+3}\\\boldsymbol{\rho}^{n+2}\\\boldsymbol{\rho}^{n+1}\end{bmatrix} = \begin{bmatrix}A_1&A_2&A_3&A_4&A_5\\ I&O&O&O&O\\ O&I&O&O&O\\ O&O&I&O&O\\ O&O&O&I&O\end{bmatrix} \begin{bmatrix}\boldsymbol{\rho}^{n+4}\\\boldsymbol{\rho}^{n+3}\\\boldsymbol{\rho}^{n+2}\\\boldsymbol{\rho}^{n+1}\\\boldsymbol{\rho}^{n}\end{bmatrix},\] where \[\begin{align} A_1 &=\left(I+\omega^{(4)}_5\Delta tQ_L+(\omega^{(4)}_5\Delta tQ_L)^2+(\omega^{(4)}_5\Delta tQ_L)^3\right) \omega^{(4)}_4 \Delta tQ_{U^{(4)}L,1}+ (\omega^{(4)}_5\Delta tQ_L)^4,\\ A_k &= \left(I+\omega^{(4)}_5\Delta tQ_L+(\omega^{(4)}_5\Delta tQ_L)^2 +(\omega^{(4)}_5\Delta tQ_L)^3\right) \omega^{(4)}_{5-k} \Delta tQ_{U^{(4)}L,k}, \quad \text{for }k = 2,3,4,\\ A_5 &=\left(I+\omega^{(4)}_5\Delta tQ_L+(\omega^{(4)}_5\Delta tQ_L)^2 +(\omega^{(4)}_5\Delta tQ_L)^3\right) \left( P^{(4)}_5 +\omega^{(4)}_0 \Delta tQ_{U^{(4)}L,5} \right) . \end{align}\]

Define \[\begin{align} C& := \omega^{(4)}_5 \frac{\Delta t}{T_1} \sum_{j=0}^2 \left(\omega^{(4)}_5\right)^j\left(\frac{\Delta t}{T_2}\right)^j, \\ B^{(4)}& := \sum_{j=0}^3 \left(\omega^{(4)}_5\right)^j\left(\frac{\Delta t}{T_2}\right)^j. \end{align}\]

Then we have \[\begin{align} A_1 &= \omega^{(4)}_4\begin{bmatrix}0&0&0& \frac{\Delta t}{T_1} +\frac{\Delta t}{T_2} |\varsigma^{(4)} ( 1)|^2 C +\cfrac{(\omega^{(4)}_5\Delta t )^4} {T_1(T_2)^3} \\0&0&0&0 \\0&0&0&0 \\0&0&0& \frac{\Delta t}{T_2}|\varsigma^{(4)} ( 1)|^2 B^{(4)}+\cfrac{(\omega^{(4)}_5\Delta t )^4} {(T_2)^4}\end{bmatrix} \\ A_k &= \omega^{(4)}_{5-k} \begin{bmatrix}0&0&0&\frac{\Delta t}{T_1} +\frac{\Delta t}{T_2} |\varsigma^{(4)} ( k)|^2 C\\0&0&0&0\\0&0&0&0\\0&0&0&\frac{\Delta t}{T_2}|\varsigma^{(4)} ( k)|^2 B^{(4)}\end{bmatrix}, \text{ for } k=2,3,4\\ A_5 & =\begin{bmatrix}1&0&0&\overline{\varsigma^{(4)} ( 5) }C +\omega^{(4)}_0 \frac{\Delta t}{T_1} +\omega^{(4)}_0\frac{\Delta t}{T_2} |\varsigma^{(4)} ( 5)|^2 C \\0&\overline{\varsigma^{(4)} ( 5)}&0&0 \\0&0&\varsigma^{(4)} ( 5)&0\\0&0&0&|\varsigma^{(4)} ( 5)|^2B^{(4)} +\omega^{(4)}_0\frac{\Delta t}{T_2}|\varsigma^{(4)} ( 5)|^2 B^{(4)}\end{bmatrix}. \end{align}\] The matrix \[\begin{bmatrix}A_1&A_2&A_3&A_4&A_5\\ I&O&O&O&O\\ O&I&O&O&O\\ O&O&I&O&O\\ O&O&O&I&O\end{bmatrix},\] is similar to the matrix \[\begin{bmatrix} O & I & O & O & O\\ O & O & I & O & O\\ O & O & O & I & O\\ O & O & O & O & I\\ A_5 & A_4 & A_3 & A_2 & A_1\end{bmatrix}.\] Consequently, the two matrices share the same eigenvalues. The latter matrix is the block companion matrix of the matrix polynomial \[L(\lambda)= \lambda^5 - A_1\lambda^4 - A_2\lambda^3 - A_3\lambda^2 - A_4\lambda - A_5 .\]

Since each of the five matrices \(A_1, \ldots, A_5\) is upper triangular, the matrix polynomial \(L(\lambda)\) is also upper triangular. Consequently, the diagonal entries of \(L(\lambda)\) are \[\lambda^5-1,\qquad \lambda^5-\overline{\varsigma^{(4)}(5)},\qquad \lambda^5-\varsigma^{(4)}(5),\] and \[\lambda^5- \sum_{k=1}^{5}\lambda^{5-k}\omega^{(4)}_{5-k} \frac{\Delta t}{T_2}|\varsigma^{(4)}(k)|^2 B^{(4)} -\lambda^4 \left(\cfrac{\omega^{(4)}_5\Delta t}{T_2}\right)^4 -|\varsigma^{(4)}(5)|^2 B^{(4)}.\] Therefore, the stability condition is characterized by the roots of these diagonal entries, which proves Theorem 1.

10 Proof of Theorem 2↩︎

After vectorizing the matrices \(\rho^{n+n_p,k}\), we express the \(p\)-th-order approximation in [eq:NPI_general_new] in the following form: \[\begin{align} \boldsymbol{\rho}^{n+n_p,0} & = \boldsymbol{\rho}^{n+n_p-1}, \\ \boldsymbol{\rho}^{n+n_p,k} &= P_{n_p}^{(p)}\boldsymbol{\rho}^n + \Delta t \sum_{i=0}^{2p-4} \omega^{(p)}_iQ_{U^{(p)}L,n_p-i}\boldsymbol{\rho}^{n+i} + \Delta t\omega^{(p)}_{n_p} Q_L \boldsymbol{\rho}^{n+n_p,k-1}, \end{align}\] for \(k=1,2,\cdots,p\). Introducing the augmented vector \[\begin{bmatrix}\boldsymbol{\rho}^{n+n_p}\\\boldsymbol{\rho}^{n+n_p-1} \\\vdots \\\boldsymbol{\rho}^{n+2}\\\boldsymbol{\rho}^{n+1} \end{bmatrix} := \begin{bmatrix}\boldsymbol{\rho}^{n+n_p,p}\\\boldsymbol{\rho}^{n+n_p-1} \\\vdots \\\boldsymbol{\rho}^{n+2}\\\boldsymbol{\rho}^{n+1} \end{bmatrix},\] the scheme can be rewritten as the one-step system

\[\begin{align} \begin{bmatrix}\boldsymbol{\rho}^{n+n_p}\\\boldsymbol{\rho}^{n+n_p-1} \\\vdots \\\boldsymbol{\rho}^{n+3}\\\boldsymbol{\rho}^{n+2}\\\boldsymbol{\rho}^{n+1} \end{bmatrix} &= \begin{bmatrix}A_1&A_2&\cdots&A_{n_p-2}&A_{n_p-1}&A_{n_p}\\ I&O&\cdots&O&O&O\\ O&I&\cdots&O&O&O\\ \vdots&\vdots&\ddots &\vdots&\vdots&\vdots\\ O&O&\cdots&I&O&O\\ O&O&\cdots&O&I&O\end{bmatrix}\begin{bmatrix}\boldsymbol{\rho}^{n+n_p-1} \\\boldsymbol{\rho}^{n+n_p-2}\\\vdots \\\boldsymbol{\rho}^{n+2} \\\boldsymbol{\rho}^{n+1} \\\boldsymbol{\rho}^{n}\end{bmatrix}, \end{align}\]

where \[\begin{align} A_1 &=\left(I+\sum_{i=1}^{p-1} (\omega^{(p)}_{n_p}\Delta tQ_L)^i \right) \omega^{(p)}_{n_p-1} \Delta tQ_{U^{(p)}L,1} + \left(\omega^{(p)}_{n_p}\Delta tQ_L\right)^p, \\ A_k &= \left(I+\sum_{i=1}^{p-1} (\omega^{(p)}_{n_p}\Delta tQ_L)^k \right) \omega^{(p)}_{n_p-k} \Delta tQ_{U^{(p)}L,k}, \text{ for } k = 2,3,\cdots,n_p-1,\\ A_{n_p} &= \left(I+\sum_{k=1}^{p-1} (\omega^{(p)}_{n_p}\Delta tQ_L)^k \right) (P^{(p)}_{n_p}+ \omega^{(p)}_0 \Delta tQ_{U^{(p)}L,n_p}). \end{align}\] Define \[\begin{align} C^{(p)} : & = \sum_{k=0}^{p-2} \dfrac{\omega^{(p)}_{n_p}\Delta t}{T_1}\left(\dfrac{\omega^{(p)}_{n_p}\Delta t}{T_2} \right)^{k-1}, \\ B^{(p)}: & = \sum_{k=0}^{p-1} \left(\dfrac{\omega^{(p)}_{n_p}\Delta t}{T_2} \right)^k. \end{align}\] Then \[\begin{align} A_1 &=\omega^{(p)}_{n_p-1} \begin{bmatrix}0&0&0& \frac{\Delta t}{T_1} +\frac{\Delta t}{T_2} |\varsigma^{(p)} ( 1)|^2 C^{(p)} +\cfrac{(\omega^{(p)}_{n_p}\Delta t )^p} {T_1(T_2)^{p-1}} \\ 0&0&0&0\\0&0&0&0\\0&0&0& \frac{\Delta t}{T_2}|\varsigma^{(p)} ( 1)|^2 B^{(p)}+\cfrac{(\omega^{(p)}_{n_p}\Delta t )^p} {(T_2)^p}\end{bmatrix}, \\ A_k &=\omega^{(p)}_{n_p-k} \begin{bmatrix}0&0&0&\frac{\Delta t}{T_1} +\frac{\Delta t}{T_2} |\varsigma^{(p)} ( k)|^2 C^{(p)}\\ 0&0&0&0\\0&0&0&0\\0&0&0&\frac{\Delta t}{T_2}|\varsigma^{(p)} ( k)|^2 B^{(p)}\end{bmatrix}, \text{ for } k = 2,\cdots,n_p-1,\\ A_{n_p} &= \begin{bmatrix}1&0&0&\overline{\varsigma^{(p)} (n_p)}C^{(p)} +\omega^{(p)}_0 \frac{\Delta t}{T_1} +\omega^{(p)}_0\frac{\Delta t}{T_2} |\varsigma^{(p)} ( n_p)|^2 C^{(p)}\\0&\overline{\varsigma^{(p)} (n_p)}&0&0\\0&0&\varsigma^{(p)} (n_p)&0\\0&0&0&|\varsigma^{(p)} (n_p)|^2B^{(p)}+\omega^{(p)}_0\frac{\Delta t}{T_2}|\varsigma^{(p)} ( n_p)|^2 B^{(p)}\end{bmatrix} . \end{align}\]

Let \(\mathcal{M}^{(p)}= \begin{bmatrix}A_1&A_2&\cdots&A_{n_p-2}&A_{n_p-1}&A_{n_p}\\ I&O&\cdots&O&O&O\\ O&I&\cdots&O&O&O\\ \vdots&\vdots&\ddots &\vdots&\vdots&\vdots\\ O&O&\cdots&I&O&O\\ O&O&\cdots&O&I&O\end{bmatrix}\) be the amplification matrix. By a block permutation, \(\mathcal{M}^{(p)}\) is equivalent to \[\begin{bmatrix} O &I & O & \cdots&O&O\\ O & O &I & \cdots &O&O\\ \vdots & \vdots & \vdots & \ddots & O&O \\ O & O & O &\cdots &I & O \\ O & O & O &\cdots & O & I\\ A_{n_p} & A_{n_p-1} & A_{n_p-2} &\cdots & A_2 & A_1\end{bmatrix} ,\] which is the block companion matrix associated with the matrix polynomial \[L(\lambda) = \lambda^{n_p} - \sum_{k=1}^{n_p} A_k\lambda^{n_p-k}.\] Since each \(A_k\), \(k=1,2,\dots,n_p\), is upper triangular, the matrix polynomial \(L(\lambda)\) is also upper triangular. Its diagonal entries are \[\lambda^{n_p}-1,\qquad \lambda^{n_p}-\overline{\varsigma^{(p)}(n_p)},\qquad \lambda^{n_p}-\varsigma^{(p)}(n_p),\] and \[\lambda^{n_p} -\sum_{k=1}^{n_p} \lambda^{n_p-k}\Delta t\,\omega^{(p)}_{n_p-k}\frac{\Delta t}{T_2}|\varsigma^{(p)} ( k)|^2 B^{(p)} -\lambda^{n_p-1}\bigl(\omega^{(p)}_{n_p}\Delta t\bigr)^p\frac{1}{(T_2)^p} -\left| \varsigma^{(p)}(n_p)\right|^2 B^{(p)}.\] Therefore, the stability condition is determined by the roots of these diagonal polynomials. This completes the proof of Theorem 2.

References↩︎

---

1. Department of Mathematics, Virginia Tech, Blacksburg, VA 24060 ().↩︎

2. Submitted to the SISC Special Section on Quantum Computing on 2026-06-30.↩︎