Digital quantum simulation framework for energy transport in an open quantum system

We can represent a system in interaction with the environment as an open quantum system, where the environment is modelled by a phonon bath. An open quantum system is a part of a larger closed system in Hilbert space $\mathcal{H}={\mathcal{H}}_{\mathrm{S}}\otimes {\mathcal{H}}_{\mathrm{B}}$, where ${\mathcal{H}}_{\mathrm{B}}$ is the phonon bath Hilbert space [42]. Assuming the initial state is represented by the separable density matrix ρ = ρ_S ⊗|0⟩⟨0|_B, the evolution of the total system is,

$$\begin{equation*}\rho \left(t\right)={U}_{\mathrm{S}\mathrm{B}}\left({\rho }_{\mathrm{S}}\otimes \vert 0\rangle \langle 0{\vert }_{\mathrm{B}}\right){U}_{\mathrm{S}\mathrm{B}}^{{\dagger}}.\end{equation*}$$

A partial trace over B gives the evolution of the open quantum system S,

$$\begin{equation*}\begin{aligned}\hfill {\rho }_{\mathrm{S}}\left(t\right)={\mathrm{T}\mathrm{r}}_{\mathrm{B}}\left(\rho \left(t\right)\right)& =\sum _{k}\langle k\vert {U}_{\mathrm{S}\mathrm{B}}\left({\rho }_{\mathrm{S}}\otimes \vert 0\rangle \langle 0{\vert }_{\mathrm{B}}\right){U}_{\mathrm{S}\mathrm{B}}^{{\dagger}}\vert k\rangle \hfill \\ \hfill & =\sum _{k}\langle k\vert {U}_{\mathrm{S}\mathrm{B}}\vert 0\rangle {\rho }_{\mathrm{S}}\left(0\right)\langle 0\vert {U}_{\mathrm{S}\mathrm{B}}^{{\dagger}}\vert k\rangle \hfill \end{aligned}\end{equation*}$$

which in the form of Kraus operators M_k will be,

$$\begin{equation}{\rho }_{\mathrm{S}}\left(t\right)=M\left({\rho }_{\mathrm{S}}\left(0\right)\right)=\sum _{k}{M}_{\mathrm{k}}{\rho }_{\mathrm{S}}\left(0\right){M}_{\mathrm{k}}^{{\dagger}}\end{equation} \tag{ 4 }$$

where

$$\begin{equation}{M}_{\mathrm{k}}=\langle k\vert {U}_{\mathrm{S}\mathrm{B}}\vert 0\rangle ={\mathrm{T}\mathrm{r}}_{\mathrm{B}}\left(\vert 0\rangle \langle k\vert {U}_{\mathrm{S}\mathrm{B}}\right).\end{equation} \tag{ 5 }$$

Here |k⟩ is an orthonormal basis for ${\mathcal{H}}_{\mathrm{B}}$ and ${\sum }_{k}\enspace {M}_{\mathrm{k}}^{{\dagger}}{M}_{\mathrm{k}}=\mathbb{1}$. This can be used to obtain operators for effects of environment on the system by introducing a bath to describe the dynamics and then tracing it out. However, ENAQT is a combination of quantum and dissipative effects, described by the Lindblad master equation, equation (3). We need to solve the master equation to obtain operators for the combined dynamics.

To solve master equation, one needs to obtain operators which capture the combined effect of quantum and dissipative processes. This can be obtained by combining the operators for different processes of ENAQT. We develop the following methodology to systematically derive the operators and the dynamical equation for ENAQT.

Methodology: To model both quantum and dissipative effects we introduce evolution operators which are a combination of Kraus operators and unitary quantum evolution. Using these operators, we derive the evolution equation for density matrix. We develop the general model by first taking a toy system and step by step adding different processes to it, to arrive at the final picture. We consider a system with one quantum jump and find a model for this toy system using the following procedure: we write an evolution for system + bath which appropriately captures the quantum jump. We trace out the bath from the evolution and find the Kraus operators for quantum jump in the system. Then, we introduce unitary quantum evolution to this system in addition to the quantum jump. Thus, the new evolution operators are obtained by appropriately combining Kraus operators and unitary quantum evolution operator. We use these operators to write the discrete time evolution equation for this setup, which appropriately captures the interplay of quantum and dissipative effects.

Next, we add another quantum jump to the system to generalise the toy model to simulate multiple quantum jumps. Again, we draw a parallel from the first case and follow the above procedure to develop a model for this setup. We write the evolution operators and arrive at the dynamical equation for this setup. This gives a general model for simulating ENAQT in open quantum systems with unitary quantum evolution and multiple environment induced quantum jumps.

3.2.1. Toy model: single quantum jump

Consider a two level system with states given by |0⟩_S and |1⟩_S. Say the bath B induces a quantum jump from |0⟩_S to |1⟩_S with probability p_0→1. The system + bath evolution can be formalised as follows,

$$\begin{equation*}\begin{aligned}\hfill \vert 0{\rangle }_{\mathrm{S}}\vert 0{\rangle }_{\mathrm{B}}& \to \sqrt{1-{p}_{0\to 1}}\vert 0{\rangle }_{\mathrm{S}}\vert 0{\rangle }_{\mathrm{B}}+\sqrt{{p}_{0\to 1}}\vert 1{\rangle }_{\mathrm{S}}\vert 1{\rangle }_{\mathrm{B}}.\hfill \\ \hfill \vert 1{\rangle }_{\mathrm{S}}\vert 0{\rangle }_{\mathrm{B}}& \to \vert 1{\rangle }_{\mathrm{S}}\vert 0{\rangle }_{\mathrm{B}}.\hfill \end{aligned}\end{equation*}$$

The Kraus operators for quantum jumps on the system, obtained from tracing out the bath are,

M₀ = $\left(\begin{matrix}\hfill \sqrt{1-{p}_{0\to 1}}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill \end{matrix}\right)\enspace \enspace$ and M₁ = $\left(\begin{matrix}\hfill 0\hfill & \hfill 0\hfill \\ \hfill \sqrt{{p}_{0\to 1}}\hfill & \hfill 0\hfill \end{matrix}\right)$.

The operators can also be written as,

$$\begin{equation}\begin{aligned}\hfill & {M}_{0}=\sqrt{1-{p}_{0\to 1}}\vert 0\rangle \langle 0\vert +\vert 1\rangle \langle 1\vert \quad \text{and}\quad \hfill & \hfill {M}_{1}=\sqrt{{p}_{0\to 1}}\vert 1\rangle \langle 0\vert .\end{aligned}\end{equation} \tag{ 6 }$$

Now, introduce quantum evolution to the system. The Lindblad equation, with the system being subject to free Hamiltonian H_S is,

$$\begin{equation}\frac{\partial \rho \left(t\right)}{\partial t}=-\frac{i}{\hslash }\left[{H}_{\mathrm{S}},\rho \left(t\right)\right]+L\left(\rho \left(t\right)\right),\end{equation} \tag{ 7 }$$

L(ρ(t)) is due to the effect of quantum jumps, where

$$\begin{equation}L\left(\rho \left(t\right)\right)=\sum _{k}\left[{L}^{k}\rho {L}^{k{\dagger}}-\frac{1}{2}{L}^{k}{L}^{k{\dagger}}\rho -\frac{1}{2}\rho {L}^{k}{L}^{k{\dagger}}\right].\end{equation} \tag{ 8 }$$

Here L^k are quantum jump operators. The terms carry their usual meaning, first term represents the quantum jumps and the other two terms are normalisation terms for the case when the jump does not happen. The Kraus operators can be combined with unitary quantum evolution as follows to solve this master equation,

$$\begin{equation}{M}_{0}^{\prime }={M}_{0}U\qquad {M}_{1}^{\prime }={M}_{1},\end{equation} \tag{ 9 }$$

where $U={\mathrm{e}}^{-\frac{\mathrm{i}{H}_{\mathrm{S}}{\Delta}t}{\hslash }}$. It can be verified that ${M}_{0}^{\prime {\dagger}}{M}_{0}^{\prime }+{M}_{1}^{\prime {\dagger}}{M}_{1}^{\prime }=\mathbb{1}$. These evolution operators, equation (9), where M₀' has unitary quantum evolution and M₁' does not, will be useful to capture the appropriate dynamics, as will be proved in the following section. The operators can be interpreted as follows. Two processes are happening in the system at any time t: coherent evolution and quantum jumps. From a state, the population density can move out of the state via these two processes. From the total population, some goes out via quantum jumps. And from the population remaining after the quantum jumps, some goes out via quantum coherence. The second operator M₁' depicts the moving out through quantum jumps. M₀ in the first operator capture's the population that remains after the quantum jump. From this remaining population, some moves out via quantum coherence as captured by U in M₀' = M₀ U.

Thus, the discrete time density matrix evolution equation obtained from these operators is given by,

$$\begin{equation}\rho \left(t+{\Delta}t\right)={M}_{0}^{\prime }\left(\sqrt{{\Delta}t}\right)\rho \left(t\right){M}_{0}^{\prime {\dagger}}\left(\sqrt{{\Delta}t}\right)+{M}_{1}^{\prime }\left(\sqrt{{\Delta}t}\right)\rho \left(t\right){M}_{1}^{\prime {\dagger}}\left(\sqrt{{\Delta}t}\right).\end{equation} \tag{ 10 }$$

Here $\sqrt{{\Delta}t}$ is taken, so that considering contributions from M₀' and ${M}_{0}^{\prime {\dagger}}$ and assuming linear dependence on time, the net change is of the first order in Δt.

Substituting equation (6) in equation (10),

$$\begin{align}\hfill \rho \left(t+{\Delta}t\right)& =\left(1-{p}_{0\to 1}\right)\vert 0\rangle \langle 0\vert U\rho \left(t\right){U}^{{\dagger}}\vert 0\rangle \langle 0\vert +\sqrt{1-{p}_{0\to 1}}\vert 0\rangle \langle 0\vert U\rho \left(t\right){U}^{{\dagger}}\vert 1\rangle \langle 1\vert \hfill \\ \hfill & \quad +\sqrt{1-{p}_{0\to 1}}\vert 1\rangle \langle 1\vert U\rho \left(t\right){U}^{{\dagger}}\vert 0\rangle \langle 0\vert +\vert 1\rangle \langle 1\vert U\rho \left(t\right){U}^{{\dagger}}\vert 1\rangle \langle 1\vert +\left({p}_{0\to 1}\right)\vert 1\rangle \langle 0\vert \rho \left(t\right)\vert 0\rangle \langle 1\vert .\hfill \end{align} \tag{ 11 }$$

The second and the third terms of the preceding equation capture the decoherence due to the bath and contribute to the normalisation terms. The first term and the last term are due to the quantum jump. Correspondence of the different terms to the Lindblad equation given above can be seen here. U and U^† capture the quantum coherence.

3.2.2. General model: multiple quantum jumps

We can generalise the previous setup by including a quantum jump in the other direction, from |1⟩_S to |0⟩_S with probability p_1→0. This can be represented using a two qubit bath, where the first qubit (B₁) captures the first quantum jump and qubit, B₂ captures the second quantum jump. Then the system + bath evolution is given by,

$$\begin{equation}\begin{aligned}\hfill \vert 0{\rangle }_{\mathrm{S}}\vert 0{\rangle }_{{\mathrm{B}}_{1}}\vert 0{\rangle }_{{\mathrm{B}}_{2}}\to \sqrt{1-{p}_{0\to 1}}\vert 0{\rangle }_{\mathrm{S}}\vert 0{\rangle }_{{\mathrm{B}}_{1}}\vert 0{\rangle }_{{\mathrm{B}}_{2}}+\sqrt{{p}_{0\to 1}}\vert 1{\rangle }_{\mathrm{S}}\vert 1{\rangle }_{{\mathrm{B}}_{1}}\vert 0{\rangle }_{{\mathrm{B}}_{2}}\\ \hfill \vert 1{\rangle }_{\mathrm{S}}\vert 0{\rangle }_{{\mathrm{B}}_{1}}\vert 0{\rangle }_{{\mathrm{B}}_{2}}\to \sqrt{1-{p}_{1\to 0}}\vert 1{\rangle }_{\mathrm{S}}\vert 0{\rangle }_{{\mathrm{B}}_{1}}\vert 0{\rangle }_{{\mathrm{B}}_{2}}+\sqrt{{p}_{1\to 0}}\vert 0{\rangle }_{\mathrm{S}}\vert 0{\rangle }_{{\mathrm{B}}_{1}}\vert 1{\rangle }_{{\mathrm{B}}_{2}}\end{aligned}\enspace ;\end{equation} \tag{ 12 }$$

and the corresponding Kraus operators for the system, after tracing out the bath are,

$$\begin{equation}\begin{aligned}\hfill {M}_{00}=\sqrt{1-{p}_{0\to 1}}\vert 0\rangle \langle 0\vert & \qquad {M}_{01}=\sqrt{{p}_{0\to 1}}\vert 1\rangle \langle 0\vert \hfill \\ \hfill {M}_{11}=\sqrt{1-{p}_{1\to 0}}\vert 1\rangle \langle 1\vert & \qquad {M}_{10}=\sqrt{{p}_{1\to 0}}\vert 0\rangle \langle 1\vert .\hfill \end{aligned}\end{equation} \tag{ 13 }$$

Adding quantum evolution to Kraus operators when the system is also subject to free Hamiltonian H_S, similar to the previous case leads to the following evolution operators,

$$\begin{equation}\begin{aligned}\hfill {M}_{00}^{\prime }={M}_{00}U\qquad & {M}_{01}^{\prime }={M}_{01}\hfill \\ \hfill {M}_{11}^{\prime }={M}_{11}U\qquad & {M}_{10}^{\prime }={M}_{10}.\hfill \end{aligned}\end{equation} \tag{ 14 }$$

Drawing parallel from equation (11), the discrete time dynamical equation is given by,

$$\begin{align}\hfill \rho \left(t+{\Delta}t\right)& =\left(1-{p}_{0\to 1}\right)\vert 0\rangle \langle 0\vert U\rho \left(t\right){U}^{{\dagger}}\vert 0\rangle \langle 0\vert +\sqrt{1-{p}_{0\to 1}}\sqrt{1-{p}_{1\to 0}}\vert 0\rangle \langle 0\vert U\rho \left(t\right){U}^{{\dagger}}\vert 1\rangle \langle 1\vert \hfill \\ \hfill & \quad +\sqrt{1-{p}_{0\to 1}}\sqrt{1-{p}_{1\to 0}}\vert 1\rangle \langle 1\vert U\rho \left(t\right){U}^{{\dagger}}\vert 0\rangle \langle 0\vert +\left(1-{p}_{1\to 0}\right)\vert 1\rangle \langle 1\vert U\rho \left(t\right){U}^{{\dagger}}\vert 1\rangle \langle 1\vert \hfill \\ \hfill & \quad +\left({p}_{0\to 1}\right)\vert 1\rangle \langle 0\vert \rho \left(t\right)\vert 0\rangle \langle 1\vert +\left({p}_{1\to 0}\right)\vert 0\rangle \langle 1\vert \rho \left(t\right)\vert 1\rangle \langle 0\vert .\hfill \end{align} \tag{ 15 }$$

The last two terms, and the first and the fourth terms above represent the quantum jumps. The second and third terms capture the decoherence and contribute to the normalisation terms. U and U^† capture the unitary quantum evolution.

Substituting equation (14) in equation (15), we arrive at the general model for simulating quantum jumps and unitary quantum evolution,

$$\begin{align}\hfill \rho \left(t+{\Delta}t\right)& ={M}_{00}^{\prime }\left(\sqrt{{\Delta}t}\right)\rho \left(t\right){M}_{00}^{\prime {\dagger}}\left(\sqrt{{\Delta}t}\right)+{M}_{00}^{\prime }\left(\sqrt{{\Delta}t}\right)\rho \left(t\right){M}_{11}^{\prime {\dagger}}\left(\sqrt{{\Delta}t}\right)+{M}_{11}^{\prime }\left(\sqrt{{\Delta}t}\right)\rho \left(t\right){M}_{00}^{\prime {\dagger}}\left(\sqrt{{\Delta}t}\right)\hfill \\ \hfill & \quad +{M}_{11}^{\prime }\left(\sqrt{{\Delta}t}\right)\rho \left(t\right){M}_{11}^{\prime {\dagger}}\left(\sqrt{{\Delta}t}\right)+{M}_{01}^{\prime }\left(\sqrt{{\Delta}t}\right)\rho \left(t\right){M}_{01}^{\prime {\dagger}}\left(\sqrt{{\Delta}t}\right)+{M}_{10}^{\prime }\left(\sqrt{{\Delta}t}\right)\rho \left(t\right){M}_{10}^{\prime {\dagger}}\left(\sqrt{{\Delta}t}\right).\hfill \end{align} \tag{ 16 }$$

The second and third terms above are the crucial features of this dynamical equation. These asymmetric terms with M₀₀' and M₁₁' representing decoherence make this equation conceptually different from the evolution equation obtained from Kraus operators,

$$\begin{equation*}{\rho }_{\mathrm{S}}\left(t\right)=\sum _{k}{M}_{\mathrm{k}}{\rho }_{\mathrm{S}}\left(0\right){M}_{\mathrm{k}}^{{\dagger}}.\end{equation*}$$

This evolution equation does not have asymmetric terms of the form ${M}_{i}{\rho }_{\mathrm{S}}\left(0\right){M}_{j}^{{\dagger}}$, which are present in the equation we have arrived at, equation (16) by drawing a parallel, equation (15) from the case with single quantum jumps, equation (11). Now we can see that this parallelism gave us insight to add these asymmetric terms to equation (14), which would have been hard to see directly from the evolution equation for Kraus operators, equation (4). This shows that the methodology used in developing the model is effective to flesh out the finer details of the theoretical framework. Equation (14) can be applied to digitally simulate any open quantum system with environment assisted evolution. It is the discrete time solution of the Lindblad master equation.

The equivalence of the dynamical equations, equations (10) and (14) to the Lindblad formalism, equation (7) is proven in this section. Considering the case with single quantum jump, the following derivation shows that equation (10) is the discrete time solution to Lindblad master equation governing the system.

The master equation is in continuous time formalism. For continuous time limit of the discrete dynamics substituting Δt → ∂t in equation (10),

$$\begin{equation}\rho \left(t+\partial t\right)={M}_{0}^{\prime }\left(\sqrt{\partial t}\right)\rho \left(t\right){M}_{0}^{\prime {\dagger}}\left(\sqrt{\partial t}\right)+{M}_{1}^{\prime }\left(\sqrt{\partial t}\right)\rho \left(t\right){M}_{1}^{\prime {\dagger}}\left(\sqrt{\partial t}\right).\end{equation} \tag{ 17 }$$

Using equation (9),

$$\begin{align}\hfill {M}_{0}^{\prime }\left(\sqrt{\partial t}\right)& ={M}_{0}U={M}_{0}\left(\sqrt{\partial t}\right){\mathrm{e}}^{-\frac{\mathrm{i}H\partial t}{\hslash }}\hfill \\ \hfill & \approx {M}_{0}\left(\sqrt{\partial t}\right)\left[\mathbb{1}-\frac{iH\partial t}{\hslash }\right].\hfill \end{align} \tag{ 18 }$$

Now, ${M}_{0}{M}_{0}^{{\dagger}}+{M}_{1}{M}_{1}^{{\dagger}}=\mathbb{1}$ and ${M}_{0}={M}_{0}^{{\dagger}}$, so ${M}_{0}\left(\sqrt{\partial t}\right)$ can be written as,

$$\begin{align}\hfill {M}_{0}\left(\sqrt{\partial t}\right)& \approx \sqrt{\mathbb{1}-{M}_{1}\left(\sqrt{\partial t}\right){M}_{1}^{{\dagger}}\left(\sqrt{\partial t}\right)}\hfill \\ \hfill & \approx \mathbb{1}-\frac{1}{2}{M}_{1}\left(\sqrt{\partial t}\right){M}_{1}^{{\dagger}}\left(\sqrt{\partial t}\right).\hfill \end{align} \tag{ 19 }$$

Substituting equation (19) in equation (18),

$$\begin{align}\hfill {M}_{0}^{\prime }\left(\sqrt{\partial t}\right)& \approx {M}_{0}\left(\sqrt{\partial t}\right)\left[\mathbb{1}-\frac{iH\partial t}{\hslash }\right]\hfill \\ \hfill & \approx \left[\mathbb{1}-\frac{1}{2}{M}_{1}\left(\sqrt{\partial t}\right){M}_{1}^{{\dagger}}\left(\sqrt{\partial t}\right)\right]\left[\mathbb{1}-\frac{iH\partial t}{\hslash }\right].\hfill \end{align} \tag{ 20 }$$

Considering terms only up to first order in ∂t, we get,

$$\begin{equation}\begin{aligned}\hfill {M}_{0}^{\prime }\left(\sqrt{\partial t}\right)& \approx \left[\mathbb{1}-\frac{iH\partial t}{\hslash }-\frac{1}{2}{M}_{1}\left(\sqrt{\partial t}\right){M}_{1}^{{\dagger}}\left(\sqrt{\partial t}\right)\right]\hfill \\ \hfill {M}_{0}^{\prime {\dagger}}\left(\sqrt{\partial t}\right)& \approx \left[\mathbb{1}+\frac{iH\partial t}{\hslash }-\frac{1}{2}{M}_{1}\left(\sqrt{\partial t}\right){M}_{1}^{{\dagger}}\left(\sqrt{\partial t}\right)\right].\hfill \end{aligned}\end{equation} \tag{ 21 }$$

Substituting M₀', equation (21) and M₁', equation (9) in equation (17),

$$\begin{align}\hfill \rho \left(t+\partial t\right)& \approx \left[\left(\mathbb{1}-\frac{iH\partial t}{\hslash }-\frac{1}{2}{M}_{1}\left(\sqrt{\partial t}\right){M}_{1}^{{\dagger}}\left(\sqrt{\partial t}\right)\right){\ast}\rho \left(t\right)\left(\mathbb{1}+\frac{iH\partial t}{\hslash }-\frac{1}{2}{M}_{1}\left(\sqrt{\partial t}\right){M}_{1}^{{\dagger}}\left(\sqrt{\partial t}\right)\right)\right]\hfill \\ \hfill & \quad +{M}_{1}\left(\sqrt{\partial t}\right)\rho \left(t\right){M}_{1}^{{\dagger}}\left(\sqrt{\partial t}\right).\hfill \end{align} \tag{ 22 }$$

Considering terms only up to first order in ∂t and rearranging the terms,

$$\begin{align}\hfill \rho \left(t+\partial t\right)& =\rho \left(t\right)-\frac{i}{\hslash }\partial t\left[H,\rho \left(t\right)\right]+{M}_{1}\left(\sqrt{\partial t}\right)\rho \left(t\right){M}_{1}^{{\dagger}}\left(\sqrt{\partial t}\right)-\frac{1}{2}{M}_{1}\left(\sqrt{\partial t}\right){M}_{1}^{{\dagger}}\left(\sqrt{\partial t}\right)\rho \left(t\right)\hfill \\ \hfill & \quad -\frac{1}{2}\rho \left(t\right){M}_{1}\left(\sqrt{\partial t}\right){M}_{1}^{{\dagger}}\left(\sqrt{\partial t}\right).\hfill \end{align} \tag{ 23 }$$

We can see that the evolution operators, equation (9) helped us get the quantum jump and the normalisation terms correctly. Padding M₁' as well with unitary quantum evolution (U) would lead to extra unwanted terms.

Set ${M}_{1}\left(\sqrt{\partial t}\right)={L}_{1}\sqrt{\partial t}$,

$$\begin{equation}\rho \left(t+\partial t\right)=\rho \left(t\right)+\partial t\left[-\frac{i}{\hslash }\left[H,\rho \left(t\right)\right]+{L}_{1}\rho \left(t\right){L}_{1}^{{\dagger}}-\frac{1}{2}{L}_{1}{L}_{1}^{{\dagger}}\rho \left(t\right)-\frac{1}{2}\rho \left(t\right){L}_{1}{L}_{1}^{{\dagger}}\right].\end{equation} \tag{ 24 }$$

We obtain the Lindblad equation,

$$\begin{equation*}\frac{\partial \rho \left(t\right)}{\partial t}=-\frac{i}{\hslash }\left[{H}_{\mathrm{S}},\rho \left(t\right)\right]+L\left(\rho \left(t\right)\right)\left(L\left(\rho \left(t\right)\right)={L}_{1}\rho \left(t\right){L}_{1}^{{\dagger}}-\frac{1}{2}{L}_{1}{L}_{1}^{{\dagger}}\rho \left(t\right)-\frac{1}{2}\rho \left(t\right){L}_{1}{L}_{1}^{{\dagger}}\right).\end{equation*}$$

We started from the dynamical equation, equation (10) and arrived at the Lindblad master equation, equation (7). Thus, for the case with single quantum jump, we showed that the dynamical model is the solution for its Lindblad equation in the Markov approximation. A similar proof for the more general case with multiple quantum jumps, equation (16) can be worked out.

Thus, equations (14) and (16) can be used to describe the complete dynamics of environment assisted quantum walk in general open quantum systems.

The dynamics of ENAQT depend on the level of noise the system is subject to. Noise depends on the strength with which the system couples to the bath. At lower level of coupling, the strength is a fraction of the full coupling mode. We can denote this fraction as χ ∈ [0, 1]. The coupling strength dependent master equation can be written as,

$$\begin{equation}\frac{\partial \rho \left(t\right)}{\partial t}=-\frac{i}{\hslash }\left[{H}_{\mathrm{c}},\rho \left(t\right)\right]+\chi {L}_{p}\left(\rho \left(t\right)\right).\end{equation} \tag{ 25 }$$

This can be modified to be written as,

$$\begin{equation}=-\left(1-\chi \right)\frac{i}{\hslash }\left[{H}_{\mathrm{c}},\rho \left(t\right)\right]+\chi \left[-\frac{i}{\hslash }\left[{H}_{\mathrm{c}},\rho \left(t\right)\right]+{L}_{p}\left(\right. \rho \left(t\right)\right].\end{equation} \tag{ 26 }$$

For χ = 1 this reduces to the normal master equation. For other values of χ, it captures the dynamics at different strengths of system–bath coupling. This can be digitally simulated by the following dynamical equation,

$$\begin{equation}{\rho }_{\chi }\left(t+{\Delta}t\right)=\left(1-\chi \right)U\rho \left(t\right){U}^{{\dagger}}+\chi \rho \left(t+{\Delta}t\right),\end{equation} \tag{ 27 }$$

where ρ(t + Δt) is given by equation (16). χ can be used to study variation of the dynamics with respect to changes in the system–bath interaction. Thus we obtain the discrete time dynamical equation for digital quantum simulation of environment assisted quantum walk with variable bath coupling. In the next section, we apply this model to simulate the FMO complex.

FMO complex is the quantum transport channel for EET in green sulphur bacteria. The complex contains chromophores, which act as sites for excitons and are held by a protein scaffold at the right distances and orientations for efficient energy transfer. The sites show quantum coherence. The Hamiltonian [22] for the multi chromophoric system is given by,

$$\begin{equation}{H}_{\mathrm{c}}=\sum _{m=1}^{{N}_{\mathrm{c}}}{\varepsilon }_{m}{a}_{m}^{{\dagger}}{a}_{m}+\sum _{n{< }m}^{{N}_{\mathrm{c}}}{V}_{mn}\left({a}_{m}^{{\dagger}}{a}_{n}+{a}_{n}^{{\dagger}}{a}_{m}\right).\end{equation} \tag{ 28 }$$

N_c = 7 is the number of chromophores in FMO complex. The ${a}_{m}^{{\dagger}}$ and a_m are the creation and annihilation operators for an electron–hole pair (exciton) at chromophore m and ɛ_m are the site energies. V_mn are Coulomb couplings of the transition densities of the chromophores. At any time there is one exciton in the complex. Initial excitation occurs at site 1 or 6 and is transported to the sink at sites 3 and 4. The structure of the FMO complex [43] is shown in figure 1. Dominant couplings are represented by edges in figure 1 for which V_mn is large. The channel is subjected to noise by the environment. The phonon bath (protein scaffold) induces quantum jumps, decoherence and dephasing of excitons without changing the number of excitations. The phonon coupling Hamiltonian is,

$$\begin{equation*}{H}_{\mathrm{p}}=\sum _{m,n}^{{N}_{\mathrm{c}}}{q}_{mn}^{p}{a}_{m}^{{\dagger}}{a}_{n}.\end{equation*}$$

Damping of the excitation due to the radiation field is given by the Hamiltonian,

$$\begin{equation*}{H}_{\mathrm{r}}=\sum _{m}^{{N}_{\mathrm{c}}}{q}_{m}^{r}\left({a}_{m}^{{\dagger}}+{a}_{m}\right).\end{equation*}$$

Lamb shifts due to phonon and photon bath coupling contribute negligibly to the dynamics [44], so are excluded from the equations. The Lindblad master equation in the Born-Markov and secular approximations is given by,

$$\begin{equation}\frac{\partial \rho \left(t\right)}{\partial t}=-\frac{i}{\hslash }\left[{H}_{\mathrm{c}},\rho \left(t\right)\right]+{L}_{p}\left(\rho \left(t\right)\right)+{L}_{r}\left(\rho \left(t\right)\right).\end{equation} \tag{ 29 }$$

The respective Lindblad superoperators L_p and L_r are given by,

$$\begin{equation}{L}_{k}\left(\rho \left(t\right)\right)=\sum _{\omega }{{\Gamma}}^{k}\left(\omega \right)\sum _{m,n}\left[{A}_{m}^{k}\left(\omega \right)\rho {A}_{n}^{k{\dagger}}\left(\omega \right)-\frac{1}{2}{A}_{m}^{k}\left(\omega \right){A}_{n}^{k{\dagger}}\left(\omega \right)\rho -\frac{1}{2}\rho {A}_{m}^{k}\left(\omega \right){A}_{n}^{k{\dagger}}\left(\omega \right)\right].\end{equation} \tag{ 30 }$$

Here (k = p, r), ${A}_{m}^{p}\left(\omega \right)={\sum }_{{\Omega}-{{\Omega}}^{\prime }=\omega }{c}_{m}^{{\ast}}\left({M}_{{\Omega}}\right){\ast}{c}_{m}\left({M}_{{{\Omega}}^{\prime }}\right)\vert {M}_{{\Omega}}\rangle \langle {M}_{{{\Omega}}^{\prime }}\vert$, where |M_Ω⟩ is the exciton with frequency Ω and |M_Ω⟩ = ∑_m c_m (M_Ω)|m⟩. The exciton states and their energies are the eigenvectors and eigenvalues, obtained by diagonalizing Hamiltonian H_c [given in equation (28)]. Excitons are delocalised over sites and the ${A}_{m}^{p}\left(\omega \right)$ represent delocalised exciton transport. The system Hamiltonian delocalises the excitons over sites due to coherence and phonon bath induces relaxation of these delocalised excitons. The jump operators constructed from exciton states capture the effect of coherence on dissipative transport. Quantum coherence causes excitons to delocalise over different sites and the quantum jumps act between these delocalised exciton states rather than between different sites. The equation has both quantum and dissipative effects and captures their interplay which leads to greater transport efficiency. Γ^p (ω) are the rates for quantum jumps [21]given by,

$$\begin{equation}{{\Gamma}}^{p}\left(\omega \right)=2\pi J\left(\omega \right)\left(1+n\left(\omega \right)\right),\end{equation} \tag{ 31 }$$

where J(ω) is Ohmic spectral density and n(ω) = 1/[exp(ℏω/kT) − 1]. The rate of quantum jumps increases with temperature. This can be understood in terms of spontaneous and induced relaxations caused by the environment. As the temperature increases, the probability of induced relaxation increases, which cause more quantum jumps and thus the phonon coupling is higher.Damping is of the order of 1 ns⁻¹ [48] and the transfer time for excitation across the complex is ∼4 ps [44]. Since the damping contribution is negligible for the duration of the quantum walk, it is neglected in the analysis. Overall, the dynamics are given by,

$$\begin{equation}\frac{\partial \rho \left(t\right)}{\partial t}\approx -\frac{i}{\hslash }\left[{H}_{\mathrm{c}},\rho \left(t\right)\right]+{L}_{p}\left(\rho \left(t\right)\right).\end{equation} \tag{ 32 }$$

The above is similar to the ENAQT quantum walk described in equation (3). Thus, the dynamics in the FMO complex are described by the master equation for ENAQT. We can treat the multichromophoric system in the FMO complex as an open quantum system. The fluctuations of correlated protein environment form a quantum bath which enhances the energy transport in the FMO complex [45, 46]. Due to coherence, excitons delocalize over multiple chromophores. This facilitates the quantum jumps between excitons. Quantum walks can give an exponential speedup over classical walks due to interference which speeds up the energy transfer [47]. In the next section, we give the theoretical framework for digital quantum simulation of the dynamics of the FMO complex using the general solution, equation (16), of the master equation developed in section 3.2.

Zoom In Zoom Out Reset image size

The quantum jumps, analogous to Lindblad operators, equation (30) show delocalised exciton transport and can be represented as follows. Say, the system is in the exciton state |M⟩. If the probability of quantum jump to state |N⟩ in time $\sqrt{{\Delta}t}$ is γ_M→N (obtained from Γ^k (ω)), for any |N⟩ ≠ |M⟩, then the Kraus operators for quantum jumps from |M⟩ are given by:

$$\begin{equation}\begin{aligned}\hfill {M}_{MN}& =\sqrt{{\gamma }_{M\to N}}\vert N\rangle \langle M\vert \quad \text{for}\enspace \text{all}\;\enspace N\ne M\hfill \\ \hfill {M}_{MM}& =\sqrt{1-\sum _{N\ne M}{\gamma }_{M\to N}}\vert M\rangle \langle M\vert .\hfill \end{aligned}\end{equation} \tag{ 33 }$$

Similar operators can be given for jumps from all such |M⟩. Using the general model, equations (14) and (16), the dynamical equation for the FMO complex is given by ($U=\mathrm{exp}\left(-\frac{i{H}_{\mathrm{c}}{\Delta}t}{\hslash }\right)$),

$$\begin{equation}\rho \left(t+{\Delta}t\right)=\sum _{M}\left[{M}_{MM}U\rho \left(t\right){U}^{{\dagger}}{M}_{MM}^{{\dagger}}+\sum _{N\ne M}\left({M}_{MN}\rho \left(t\right){M}_{MN}^{{\dagger}}+{M}_{MM}U\rho \left(t\right){U}^{{\dagger}}{M}_{NN}^{{\dagger}}\right)\right].\end{equation} \tag{ 34 }$$

As proven in section 3.3, this equation is the solution for the Lindblad equation for ENAQT in FMO complex. It should be noted that the last term in this equation is the decoherence term and the first and second terms are due to the quantum jumps. The evolution is trace preserving (${\sum }_{i,j}\enspace {M}_{ij}^{{\dagger}}{M}_{ij}=\mathbb{1}$). This is the theoretical model for digital quantum simulation of the FMO complex.

Equations (34) and (27) represent evolution of the density matrix in the FMO complex. To verify the effectiveness of these models, we simulate them numerically and present the results. Evolution of population densities at different sites is calculated using equation (34).The step size (Δt) is chosen to be 10 fs, as observed coherence time is ∼300 fs and exciton relaxation time is ∼70 fs [16]. For the numerical simulation, the value of the system Hamiltonian is taken from [21]. The quantum jump rates between excitons are calculated from relaxation data presented in [49], using theoretical analysis given in [21, 44].

4.3.1. Result: high efficiency of energy transfer

In figures 2 and 3 we show the evolution of population densities when the initial excitation is at site 1 (figure 2) and site 6 (figure 3). The efficiency attained (sum of population on site 3 and site 4 at t = 4 ps) is ∼98 percent. It is in agreement with theoretical evidence presented in [22], by using master equation model, equation (29). It can also be seen from the graphs that transfer happens faster if the initial excitation is at site 6, as observed experimentally in [44].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.3.2. Directionality in the quantum walk

In the FMO complex, directionality of energy transport is given by quantum jumps which cause the transfer of exciton from initial sites towards the sink. To show this, we simulate the evolution in absence of quantum jumps (figure 4 and 5). As we can see, the population just oscillates between initial strongly coupled sites, that is it stays at the initial exciton. This is caused by Anderson localization due to disorder in the Hamiltonian of the FMO complex. Quantum jumps help overcome the trapping of exciton due to Anderson localization by transferring the population between different excitons. The probabilities of quantum jumps are greater in the direction from initial site towards the final site, thus giving directionality to the walk. Compared to these figures, the simulations with quantum jumps (figures 2 and 3) show exciton transport. The role of noise in ENAQT is also illustrated in [50].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.3.3. Result: dependence on environment temperature

The dynamics of the FMO complex are temperature dependent. The phonon couplings vary with temperature and the coupling strength can be seen to increase as the temperature increases, equation (31). The dynamics in the FMO complex have been observed up to ambient temperatures. At lower temperatures, the coupling strength is a fraction of the maximum coupling constants. This can be digitally simulated by the dynamical equation with tunable bath coupling, equation (27). In figures 6 and 7 we show the results of simulations with different phonon couplings, equation (27). For χ = 0.06 the evolution of population densities at different sites are shown when the initial excitation at site 1 (figure 6) and site 6 (figure 7). In figure 6 it can be seen that there is an initial oscillation of population between site 1 and 2. This is due to the high coupling between these two sites which causes the exciton over site 1 to delocalise to site 2. Slowly this oscillation dies, as quantum jumps cause the population to move towards the sink at sites 3 and 4, whose population starts rising. There is some population at site 7 as well, which serves as the connecting link between different sites (as can be seen in figure 1). As compared to figure 2, there is an enhanced effect of quantum dynamics (coherence between site 1 and 2) visible in figure 6. This is expected since the strength of the bath coupling is lower for χ = 0.06, leading to suppression of environment induced quantum jumps. Similar dynamics can be seen in figure 7 where there is an initial oscillation of population between site 6 and site 5. This population is slowly transported towards the sink at site 3 and 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 6 and 7 represent the dynamics in the FMO complex at 77 K. These results can be matched with the dynamics obtained in previous theoretical study [51]. We can see that the results of calculations done using the discrete time evolution model of equations (34) and (27) match with theoretical [22, 51] and experimental studies [44] done on the FMO complex.

The usual method for implementation of the FMO complex is to treat each site (position basis) on a separate qubit, whose excitation from the ground state |0⟩ to excited state |1⟩ using creation operator (a^†) depicts population at that site. This requires n qubits with and Hilbert space of the dimension 2ⁿ for implementing a quantum system with n sites. However, the Hamiltonian of the system, equation (28) has dimension n. Mapping the Hamiltonian to the expanded basis for implementation requires complex calculations and often multiple techniques have to be employed for simulation. Additionally, for implementing the non unitary part of evolution, each qubit of the system uses an ancilla qubit to capture dephasing of the qubit. This approach gives an O(n) space complexity of implementation.

However, simulation of the framework presented in this paper can be done on a smaller qubit space. The FMO complex represents a seven site system and the state of exciton on the sites can be represented by a Hilbert space with dim = 7. Minimum number of qubits required for simulating this is ⌈log 7⌉ = 3. The seven exciton states can be represented as,

$$\begin{equation}\begin{aligned}\hfill & {E}_{1}=\vert 001\rangle \qquad {E}_{5}=\vert 101\rangle \hfill \\ \hfill & {E}_{2}=\vert 010\rangle \qquad {E}_{6}=\vert 110\rangle \hfill \\ \hfill & {E}_{3}=\vert 011\rangle \qquad {E}_{7}=\vert 111\rangle .\hfill \\ \hfill & {E}_{4}=\vert 100\rangle \hfill \end{aligned}\end{equation} \tag{ 35 }$$

We can implement one quantum jump, say from exciton state E_i to E_j as

$$\begin{equation}\rho \to {\gamma }_{i\to j}\vert j\rangle \langle i\vert \rho \vert i\rangle \langle \enspace j\vert ,\end{equation} \tag{ 36 }$$

where |i⟩ represents exciton E_i and γ_i→j is probability of quantum jump from state E_i to E_j .

Overall, the dynamics can be represented as,

$$\begin{equation}\vert i\rangle \vert 0{\rangle }_{\mathrm{B}}\to \sqrt{1-{\gamma }_{i\to j}}\vert i\rangle \vert 0{\rangle }_{\mathrm{B}}+\sqrt{{\gamma }_{i\to j}}\vert j\rangle \vert 1{\rangle }_{\mathrm{B}}.\end{equation} \tag{ 37 }$$

We can modify it to include another bath qubit as,

$$\begin{equation}\vert i\rangle \vert 0{\rangle }_{{\mathrm{B}}_{1}}\vert 0{\rangle }_{{\mathrm{B}}_{2}}\to \sqrt{1-{\gamma }_{i\to j}}\vert i\rangle \vert 0{\rangle }_{{\mathrm{B}}_{1}}\vert 0{\rangle }_{{\mathrm{B}}_{2}}+\sqrt{{\gamma }_{i\to j}}\vert j\rangle \vert 1{\rangle }_{{\mathrm{B}}_{1}}\vert 1{\rangle }_{{\mathrm{B}}_{2}}.\end{equation} \tag{ 38 }$$

This can also be written as,

$$\begin{equation}\vert 0{\rangle }_{{\mathrm{B}}_{1}}\vert i\rangle \vert 0{\rangle }_{{\mathrm{B}}_{2}}\to \sqrt{1-{\gamma }_{i\to j}}\vert 0{\rangle }_{{\mathrm{B}}_{1}}\vert i\rangle \vert 0{\rangle }_{{\mathrm{B}}_{2}}+\sqrt{{\gamma }_{i\to j}}\vert 1{\rangle }_{{\mathrm{B}}_{1}}\vert j\rangle \vert 1{\rangle }_{{\mathrm{B}}_{2}}.\end{equation} \tag{ 39 }$$

The above is for one quantum jump. We can trace out the second bath qubit and obtain the Kraus operators,

$$\begin{equation}{M}_{0}^{ij}=\left(\sqrt{1-{\gamma }_{i\to j}}\right)\vert 0,i\rangle \langle 0,i\vert +\vert 1,i\rangle \langle 1,i\vert +{I}_{{\mathrm{B}}_{1}}\otimes \sum _{m\ne i}\vert m\rangle \langle m\vert ,\end{equation} \tag{ 40 }$$

$$\begin{equation}{M}_{1}^{ij}=\sqrt{{\gamma }_{i\to j}}\vert 1,j\rangle \langle 0,i\vert ,\end{equation} \tag{ 41 }$$

where $\vert p,q\rangle =\vert p{\rangle }_{{\mathrm{B}}_{1}}\vert q\rangle ,$ for p ∈ {0, 1} and q ∈ [1,7]. The above Kraus operators for one quantum jump leave all states |m⟩ ≠ |i⟩ and the state |1, i⟩ unchanged and implement the quantum jump from |0, i⟩ → |1, j⟩. Each step of the evolution consists of multiple quantum jumps, which can be represented as,

$$\begin{equation}\begin{aligned}\hfill & {1}^{\mathrm{s}\mathrm{t}}\;\text{jump}\;\left(\text{from}\;\vert 1\rangle \to \vert 2\rangle \right):\hfill \\ \hfill & \qquad {M}_{0}=\left(\sqrt{1-{\gamma }_{1\to 2}}\right)\vert 0,1\rangle \langle 0,1\vert +\vert 1,1\rangle \langle 1,1\vert +{I}_{{\mathrm{B}}_{1}}\otimes \sum _{m\ne 1}\vert m\rangle \langle m\vert \qquad {M}_{1}=\sqrt{{\gamma }_{1\to 2}}\vert 1,2\rangle \langle 0,1\vert \hfill \\ \hfill & {2}^{\mathrm{n}\mathrm{d}}\;\text{jump}\;\left(\text{from}\;\vert 1\rangle \to \vert 3\rangle \right):\hfill \\ \hfill & \qquad {M}_{0}=\left(\sqrt{1-{\gamma }_{1\to 3}}\right)\vert 0,1\rangle \langle 0,1\vert +\vert 1,1\rangle \langle 1,1\vert +{I}_{{\mathrm{B}}_{1}}\otimes \sum _{m\ne 1}\vert m\rangle \langle m\vert \qquad {M}_{1}=\sqrt{{\gamma }_{1\to 3}}\vert 1,3\rangle \langle 0,1\vert \hfill \\ \hfill & {3}^{\mathrm{r}\mathrm{d}}\;\text{jump}\;\left(\text{from}\;\vert 1\rangle \to \vert 4\rangle \right):\hfill \\ \hfill & \dots \hfill \\ \hfill & \dots \hfill \\ \hfill & {7}^{\mathrm{t}\mathrm{h}}\;\text{jump}\;\left(\text{from}\;\vert 2\rangle \to \vert 1\rangle \right):\hfill \\ \hfill & \qquad {M}_{0}=\left(\sqrt{1-{\gamma }_{2\to 1}}\right)\vert 0,2\rangle \langle 0,2\vert +\vert 1,2\rangle \langle 1,2\vert +{I}_{{\mathrm{B}}_{1}}\otimes \sum _{m\ne 2}\vert m\rangle \langle m\vert \qquad {M}_{1}=\sqrt{{\gamma }_{2\to 1}}\vert 1,1\rangle \langle 0,2\vert \hfill \\ \hfill & {8}^{\mathrm{t}\mathrm{h}}\;\text{jump}\;\left(\text{from}\;\vert 2\rangle \to \vert 3\rangle \right):\hfill \\ \hfill & \qquad {M}_{0}=\left(\sqrt{1-{\gamma }_{2\to 3}}\right)\vert 0,2\rangle \langle 0,2\vert +\vert 1,2\rangle \langle 1,2\vert +{I}_{{\mathrm{B}}_{1}}\otimes \sum _{m\ne 2}\vert m\rangle \langle m\vert \qquad {M}_{1}=\sqrt{{\gamma }_{2\to 3}}\vert 1,3\rangle \langle 0,2\vert \hfill \\ \hfill & \dots \hfill \\ \hfill & \dots \hfill \\ \hfill & \dots \hfill \\ \hfill & 4{2}^{\mathrm{n}\mathrm{d}}\;\text{jump}\;\left(\text{from}\;\vert 7\rangle \to \vert 6\rangle \right):\hfill \\ \hfill & \qquad {M}_{0}=\left(\sqrt{1-{\gamma }_{7\to 6}}\right)\vert 0,7\rangle \langle 0,7\vert +\vert 1,7\rangle \langle 1,7\vert +{I}_{{\mathrm{B}}_{1}}\otimes \sum _{m\ne 7}\vert m\rangle \langle m\vert \qquad {M}_{1}=\sqrt{{\gamma }_{7\to 6}}\vert 1,6\rangle \langle 0,7\vert .\hfill \end{aligned}\end{equation} \tag{ 42 }$$

The following proof shows that sequentially implementing these jumps captures the Markovian evolution of the system under multiple quantum jumps. Consider any two quantum jumps |u⟩ → |v⟩ and |w⟩ → |x⟩. Sequential application of Kraus operators gives,

$$\begin{equation}\begin{aligned}\hfill \tilde {\rho }& ={M}_{0}^{uv}\left(\rho \left(0\right)\right){{M}_{0}^{uv}}^{{\dagger}}+{M}_{1}^{uv}\left(\rho \left(0\right)\right){{M}_{1}^{uv}}^{{\dagger}}\hfill \\ \hfill \rho \left(t\right)& ={M}_{0}^{wx}\left(\tilde {\rho }\right){{M}_{0}^{wx}}^{{\dagger}}+{M}_{1}^{wx}\left(\tilde {\rho }\right){{M}_{1}^{wx}}^{{\dagger}},\hfill \end{aligned}\end{equation} \tag{ 43 }$$

where ${M}_{0/1}^{uv}$ are Kraus operators for jump from |u⟩ → |v⟩. Thus,

$$\begin{align}\hfill \rho & \to {M}_{0}^{wx}{M}_{0}^{uv}\left(\rho \right){M}_{0}^{uv{\dagger}}{M}_{0}^{wx{\dagger}}+{M}_{0}^{wx}{M}_{1}^{uv}\left(\rho \right){M}_{1}^{uv{\dagger}}{M}_{0}^{wx{\dagger}}\hfill \\ \hfill & \quad +\enspace {M}_{1}^{wx}{M}_{0}^{uv}\left(\rho \right){M}_{0}^{uv{\dagger}}{M}_{1}^{wx{\dagger}}+{M}_{1}^{wx}{M}_{1}^{uv}\left(\rho \right){M}_{1}^{uv{\dagger}}{M}_{1}^{wx{\dagger}}.\hfill \end{align} \tag{ 44 }$$

Now for general {u, v} ≠ {w, x}, using equation (40)

$$\begin{equation}\begin{aligned}\hfill {M}_{0}^{wx}{M}_{0}^{uv}& =\left(\sqrt{1-{\gamma }_{w\to x}}\right)\vert 0,w\rangle \langle 0,w\vert +\left(\sqrt{1-{\gamma }_{u\to v}}\right)\vert 0,u\rangle \langle 0,u\vert +\sum _{m\ne w,u}{I}_{{\mathrm{B}}_{1}}\otimes \vert m\rangle \langle m\vert \hfill \\ \hfill {M}_{0}^{wx}{M}_{1}^{uv}& =\left(\sqrt{{\gamma }_{u\to v}}\right)\vert 1,v\rangle \langle 0,u\vert \hfill \\ \hfill {M}_{1}^{wx}{M}_{0}^{uv}& =\left(\sqrt{{\gamma }_{w\to x}}\right)\vert 1,x\rangle \langle 0,w\vert \hfill \\ \hfill {M}_{1}^{wx}{M}_{1}^{uv}& =0.\hfill \end{aligned}\end{equation} \tag{ 45 }$$

We can see that ${M}_{0}^{wx}{M}_{1}^{uv}$ and ${M}_{1}^{wx}{M}_{0}^{uv}$ represent the quantum jumps and ${M}_{0}^{wx}{M}_{0}^{uv}$ gives the population that remains for coherent evolution. This is exactly like the general model of ENAQT developed in the paper. The above operators are calculated for {u, v} ≠ {w, x}. Similar calculations can be done when this condition does not hold. Thus, sequentially applying Kraus operators gives the Markovian evolution of the system, as described by ENAQT for two quantum jumps. The same argument can be extrapolated for multiple quantum jumps implemented one after the other.

In the procedure described above, $\vert 0{\rangle }_{{\mathrm{B}}_{1}}$ is entangled with the population remaining after quantum jumps and $\vert 1{\rangle }_{{\mathrm{B}}_{1}}$ is entangled with the part of system capturing quantum jumps.

The coherent evolution due to Hamiltonian of the system can be implemented at the end of all quantum jumps using state $\vert 0{\rangle }_{{\mathrm{B}}_{1}}$ as,

$$\begin{equation}C=\vert 0{\rangle }_{{\mathrm{B}}_{1}}\langle 0{\vert }_{{\mathrm{B}}_{1}}\otimes U+\vert 1{\rangle }_{{\mathrm{B}}_{1}}\langle 1{\vert }_{{\mathrm{B}}_{1}}\otimes I,\end{equation} \tag{ 46 }$$

where C is the unitary operator that acts on B₁ and the system to evolve the coherent part of density matrix (entangled with $\vert 0{\rangle }_{{\mathrm{B}}_{1}}$) under unitary evolution $U={\text{e}}^{-\text{i}\frac{{H}_{\mathrm{c}}}{\hslash }t}$.

B₁ is traced out at the end to complete one step of evolution with both coherent evolution and multiple quantum jumps.

The main advantage of using two bath qubits is that it helps us implement simultaneous quantum jumps sequentially. This greatly simplifies the simulation procedure. Overall, the second bath qubit is used to implement quantum jumps and the first bath qubit is used to implement coherent evolution. The first qubit is traced out after one complete step of evolution consisting of coherent evolution and multiple quantum jumps (which are implemented by tracing out the second bath qubit after each quantum jump and resetting it to $\vert 0{\rangle }_{{\mathrm{B}}_{2}}$).

The circuit for simulating one quantum jump is given in figure 8. Here, an example is given for quantum jump from |1⟩ → |2⟩. ${R}_{y}^{1,2}\left(\theta \right)$ is the rotation about y-axis on the Bloch sphere such that ${R}_{y}^{1,2}\left(\theta \right)\vert 0{\rangle }_{{\mathrm{B}}_{2}}=\mathrm{cos}\enspace \frac{\theta }{2}\vert 0{\rangle }_{{\mathrm{B}}_{2}}+\mathrm{sin}\enspace \frac{\theta }{2}\vert 1{\rangle }_{{\mathrm{B}}_{2}}$ and $\mathrm{sin}\enspace \frac{\theta }{2}=\sqrt{{\gamma }_{1\to 2}}$. The evolution of state under these operations is,

$$\begin{equation}\begin{aligned}\hfill & \text{First}\;\text{gate}:\hfill \\ \hfill & \vert 0{\rangle }_{{\mathrm{B}}_{1}}\vert 001\rangle \vert 0{\rangle }_{{\mathrm{B}}_{2}}\to \sqrt{1-{\gamma }_{1\to 2}}\vert 0{\rangle }_{{\mathrm{B}}_{1}}\vert 001\rangle \vert 0{\rangle }_{{\mathrm{B}}_{2}}+\sqrt{{\gamma }_{1\to 2}}\vert 0{\rangle }_{{\mathrm{B}}_{1}}\vert 001\rangle \vert 1{\rangle }_{{\mathrm{B}}_{2}}.\hfill \\ \hfill & \text{Second}\;\text{gate}:\hfill \\ \hfill & \sqrt{{\gamma }_{1\to 2}}\vert 0{\rangle }_{{\mathrm{B}}_{1}}\vert 001\rangle \vert 1{\rangle }_{{\mathrm{B}}_{2}}\to \sqrt{{\gamma }_{1\to 2}}\vert 1{\rangle }_{{\mathrm{B}}_{1}}\vert 010\rangle \vert 1{\rangle }_{{\mathrm{B}}_{2}}\hfill \\ \hfill & \text{which}\;\text{can}\;\text{also}\;\text{be}\;\text{written}\;\text{as}\;\left(\text{using}\;\text{eq.}\left(35\right)\right)\hfill \\ \hfill & \sqrt{{\gamma }_{1\to 2}}\vert 0{\rangle }_{{\mathrm{B}}_{1}}\vert 1\rangle \vert 1{\rangle }_{{\mathrm{B}}_{2}}\to \sqrt{{\gamma }_{1\to 2}}\vert 1{\rangle }_{{\mathrm{B}}_{1}}\vert 2\rangle \vert 1{\rangle }_{{\mathrm{B}}_{2}}.\hfill \end{aligned}\end{equation} \tag{ 47 }$$

Zoom In Zoom Out Reset image size

The circuit for simulating one step (figure 9) of evolution of the FMO complex is given in figure 10, constructed by composing quantum jumps sequentially.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 11 shows the complete circuit diagram for multiple steps. Apart from the quantum jumps and coherent evolution, there is an initial gate (D) which transforms the basis from position to excitons which is inverted at the end before measurement of population at different sites. This is done to aid our calculations for quantum jumps which are done in the exciton basis. The gaps between steps in B₁ qubit are to indicate that it is traced out at the end of each step and reset to $\vert 0{\rangle }_{{\mathrm{B}}_{1}}$ at the beginning of the next step

Zoom In Zoom Out Reset image size

4.4.1Complexity analysis

Number of gates: To analyse the complexity in terms of number of gates, we need a to decompose multi qubit gates to elementary gates. Thus, for single quantum jump, figure 8 can also be drawn as figure 12. We can see that we require log n gates for the first gate [equation (47)] and log n gates for the second gate [equation (47)]. Similar decomposition can be done for all quantum jumps shown in figure 10.

Zoom In Zoom Out Reset image size

There are n(n − 1) quantum jumps, each of which requires 2 log n gates. This gives O(n² log n) gates for implementing quantum jumps. There is one gate for implementing the coherent evolution. Decomposition of unitary evolution in terms of elementary gates will take O(n) gates (decomposition of arbitrary unitary gate over q qubits takes O(exp(q)) gates, and we have q = log  n qubits). Thus the gate complexity for implementing one step of evolution is O(n² log n). The number of gates is also smaller than needed in traditional Stinespring dilation of quantum channels which often need O(n⁶) gates [52]. The total complexity for complete evolution for time = T is $O\left(\frac{T}{{\Delta}t}{n}^{2}\enspace \mathrm{log}\enspace n\right)$, where one time step = Δt.

Number of qubits: We have used log n qubits for simulating the system and 2 bath qubits. Additionally, we use log n ancillas for simulating multi qubit control gates. Thus we need O(log n) qubits for simulating the FMO complex. The extra qubits can be avoided by using Suzuki–Lie Trotter decomposition of the gates.

We have provided a model for digital quantum simulation of energy transfer in open quantum systems (the approach can be applied to systems other than FMO complex). The alternate to such simulations is to use a Stinespring dilation for getting a unitary quantum evolution of system and bath which can simultaneously implement the quantum jumps. An approximate decomposition of the unitary evolution as elementary one and two qubit gates is then obtained using techniques like Suzuki–Lie Trotter decomposition. Or each site is treated on a separate qubit along with an ancillary qubit and the effective evolution is mapped on a larger Hilbert space of 7 qubits (dim = 2⁷) + 7 ancilla qubits. However, our framework only uses ⌈log 7⌉ = 3 qubits thus giving a log-reduction in space. There is also a huge reduction is space of ancillary qubits from O(n) or O(log n) to just 2 qubits. This reduction in complexity from poly(n) to O(1) gives the absolute minimum space complexity that can implement the bath.

Thus, the discrete time evolution equation, together with a rather straightforward way of varying the coupling strength of environment provides an optimal method to digitally simulate environment assisted transport in open quantum systems.