Hybrid quantum-classical simulation of quantum speed limits in open quantum systems

TL;DR

This paper introduces a hybrid quantum-classical approach to efficiently compute quantum speed limits in open quantum systems, revealing how coupling strength and temperature influence quantum evolution speed and energy transfer efficiency.

Contribution

The paper presents a novel mixed Wigner-Heisenberg method for calculating quantum speed limits in multi-level open systems, combining quantum and classical dynamics for improved analysis.

Findings

QSL time exhibits a turnover in strong coupling regimes.

Subsystem-bath interactions significantly affect quantum evolution speed.

Potential link between QSL time and energy transfer efficiency in biological systems.

Abstract

The quantum speed limit (QSL) provides a fundamental upper bound on the speed of quantum evolution, but its evaluation in generic open quantum systems still presents a formidable computational challenge. Herein, we introduce a hybrid quantum-classical method for computing QSL times in multi-level open quantum systems. The method is based on a mixed Wigner-Heisenberg representation of the composite quantum dynamics, in which the open subsystem of interest is treated quantum mechanically and the bath is treated in a classical-like fashion. By solving a set of coupled first-order deterministic differential equations for the quantum and classical degrees of freedom, one can compute the QSL time. To demonstrate the utility of the method, we study the unbiased spin-boson model and provide a detailed analysis of the effect of the subsystem-bath coupling strength and bath temperature on the QSL time. In particular, we find a turnover of the QSL time in the strong coupling regime, which is indicative of a speed-up in the quantum evolution. We also apply the method to the Fenna- Matthews-Olson complex model and identify a potential connection between the QSL time and the efficiency of the excitation energy transfer at different temperatures.

Tables

Table 1. Table 1: Site energies E n subscript 𝐸 𝑛 E_{n} (diagonal entries) and excitonic coupling strengths V m ​ n subscript 𝑉 𝑚 𝑛 V_{mn} (off-diagonal entries) in units of cm − 1 superscript cm 1 \mathrm{cm}^{-1} for a subunit of the apo-FMO complex

BChl	1	2	3	4	5	6	7
1	12410	-87.7	5.5	-5.9	6.7	-13.7	-9.9
2	-87.7	12530	30.8	8.2	0.7	11.8	4.3
3	5.5	30.8	12210	-53.5	-2.2	-9.6	6.0
4	-5.9	8.2	-53.5	12320	-70.7	-17.0	-63.3
5	6.7	0.7	-2.2	-70.7	12480	81.1	-1.3
6	-13.7	11.8	-9.6	-17.0	81.1	12630	39.7
7	-9.9	4.3	6.0	-63.3	-1.3	39.7	12440

Equations

$${\tau }_{QSL}=\max \{{\tau }_1,{\tau }_2,{\tau }_{\infty }\},$$

$${\tau }_{QSL}=\frac{1-\mathcal{F}({\hat{\rho }}_0,{\hat{\rho }}_{\tau })}{\frac{1}{\tau }{\int }_0^{\tau }dt\sqrt{{\mathrm{Tr}}_S{\left(\frac{d}{dt}{\hat{\rho }}_t\right)}^2}}.$$

$$\mathcal{F}({\hat{\rho }}_0,{\hat{\rho }}_{\tau })=⟨{\hat{\mathcal{P}}}_{nn}(\tau )⟩,$$

$${\mathrm{Tr}}_S{\left(\frac{d}{dt}{\hat{\rho }}_t\right)}^2$$

$${\rho }_t^{mn}={\mathrm{Tr}}_S\left[{\hat{\mathcal{P}}}_{nm}{\hat{\rho }}_t\right]=⟨{\hat{\mathcal{P}}}_{nm}(t)⟩.$$

$$\hat{H}={\hat{H}}_S(\hat{\mathbf{x}})+{\hat{H}}_I(\hat{\mathbf{x}},\hat{\mathbf{X}})+{\hat{H}}_B(\hat{\mathbf{X}}).$$

$${\hat{A}}_W(\mathbf{X})=\int d\mathbf{z}{e}^{i\mathbf{P}\cdot \mathbf{z}/\hslash }⟨\mathbf{R}-\frac{\mathbf{z}}{2}\left|\hat{A}\right|\mathbf{R}+\frac{\mathbf{z}}{2}⟩,$$

$${\hat{H}}_W(\hat{\mathbf{x}},\mathbf{X})={\hat{H}}_S(\hat{\mathbf{x}})+{\hat{H}}_{I,W}(\hat{\mathbf{x}},\mathbf{X})+H_{B,W}(\mathbf{X}).$$

$$\frac{d}{dt}\hat{\mathbf{x}}(t)$$

$$\frac{d}{dt}\mathbf{X}(t)$$

$$\frac{d}{dt}{\mathbf{x}}^{\alpha {\alpha }^{\prime }}(t)$$

$$\frac{d}{dt}{\mathbf{X}}^{\alpha {\alpha }^{\prime }}(t)$$

$$\frac{d}{dt}⟨{\hat{\mathcal{P}}}_{nm}(t)⟩=\sum _{\alpha {\alpha }^{\prime }}\int d\mathbf{X}(0){\rho }_{B,W}(\mathbf{X}(0))\frac{d}{dt}{\mathcal{P}}_{nm}^{\alpha {\alpha }^{\prime }}(t){\rho }_0^{{\alpha }^{\prime }\alpha },$$

$${\rho }_{B,W}(0)$$

$$\hat{H}=-\Delta {\hat{\sigma }}_x+\frac{1}{2}\sum _{j=1}^N\left({\hat{P}}_j^2+{\omega }_j^2{\hat{R}}_j^2-2C_j{\hat{R}}_j{\hat{\sigma }}_z\right),$$

$${\hat{H}}_W=-\Delta {\hat{\sigma }}_x+\frac{1}{2}\sum _{j=1}^N\left(P_j^2+{\omega }_j^2{R}_j^2-C_j{R}_j{\hat{\sigma }}_z-C_j{\hat{\sigma }}_z{R}_j\right).$$

$$\mathcal{F}({\hat{\rho }}_0,{\hat{\rho }}_{\tau })=⟨{\hat{\mathcal{P}}}_{++}(\tau )⟩=\frac{1}{2}(1+B_z(\tau )),$$

$${\mathrm{Tr}}_S{\left(\frac{d}{dt}{\hat{\rho }}_t\right)}^2=\frac{1}{2}\sum _{m=x,y,z}{\left(\frac{d}{dt}{B}_m(t)\right)}^2.$$

$${\tau }_{QSL}^{iso}=\frac{1}{2}\frac{1-\cos (2\Delta \tau )}{\sqrt{2}\Delta }.$$

$$B_m(t)$$

$$\frac{d}{dt}{B}_m(t)$$

$$C_j=\sqrt{\xi \hslash {\omega }_0}{\omega }_j,{\omega }_j=-{\omega }_c\ln \left(1-j\frac{{\omega }_0}{{\omega }_c}\right),$$

$$\frac{d}{dt}{B}_x(t)$$

$$\frac{d}{dt}{B}_y(t)$$

$$\frac{d}{dt}{B}_z(t)$$

$$M_1(t)+i{M}_2(t)$$

$$\frac{d}{dt}{B}_x(t)$$

$$\frac{d}{dt}{B}_y(t)$$

$$\frac{d}{dt}{B}_z(t)$$

$$\hat{H}$$

$${\hat{H}}_W$$

$$\mathcal{F}({\hat{\rho }}_0,{\hat{\rho }}_{\tau })=⟨{\hat{\mathcal{P}}}_{11}(\tau )⟩.$$

$$\mathrm{BChls}1\to 2\to 3\rightleftharpoons 4,$$

$$⟨{\hat{\mathcal{P}}}_{nm}(t)⟩$$

Full text

I Introduction

The quantum speed limit (QSL) Mandelstam and Tamm (1945); Margolus and Levitin (1998); Deffner and Campbell (2017) sets the maximum speed (or equivalently the minimum time) at which a quantum system can evolve from an initial state to a target state. For closed systems and orthogonal states, the QSL time is given by the relation $\tau_{QSL}=\max\{\pi\hbar/(2\Delta E),\pi\hbar/(2E)\}$, with the two bounds usually referred to as Mandelstam-Tamm (MT) Mandelstam and Tamm (1945) and Margolus-Levitin (ML) Margolus and Levitin (1998) types, respectively. The MT bound depends on the variance of the energy of the initial state, $\Delta E^{2}$, while the ML bound depends on the mean energy with respect to the ground state, $E\equiv\langle G|\hat{H}|G\rangle$ (with $|G\rangle$ the ground state and $\hat{H}$ the Hamiltonian). Generalizations of the MT and ML bounds to nonorthogonal states and mixed initial states can be found in Refs. Pfeifer (1993); Zwierz (2012); Barnes (2013); Poggi et al. (2013); Deffner and Lutz (2013a); Russell and Stepney (2014); Campaioli et al. (2018a). The QSL has proven to be a useful concept in the fields of quantum computation Lloyd (2000); Santos and Sarandy (2015), quantum control Caneva et al. (2009); Hegerfeldt (2013); Mukherjee et al. (2013); Hegerfeldt (2014); Deffner (2014); Campbell and Deffner (2017), quantum metrology Giovannetti et al. (2011); Chin et al. (2012); Tsang (2013), and quantum thermodynamics Deffner and Lutz (2010); Abah and Lutz (2017); Campaioli et al. (2017). Interestingly, recent studies on closed systems have shown that the QSL can have a classical counterpart Shanahan et al. (2018); Okuyama and Ohzeki (2018).

In recent years, there has been a great deal of interest in the QSLs of open quantum systems (OQSs) Taddei et al. (2013); del Campo et al. (2013); Deffner and Lutz (2013b); Xu et al. (2014); Xu and Zhu (2014); Liu et al. (2015); Sun et al. (2015); Hou et al. (2015); Cimmarusti et al. (2015); Marvian and Lidar (2015); Marvian et al. (2016); Mirkin et al. (2016); Mondal et al. (2016); Pires et al. (2016); Ektesabi et al. (2017); Campaioli et al. (2018b); Funo et al. (2018); Wu and Yu (2018). For an OQS whose reduced dynamics is governed by a generic time-convolutionless master equation for the reduced density matrix, $\frac{d}{dt}\hat{\rho}_{t}=\hat{\mathcal{L}}_{t}\hat{\rho}_{t}$ (where $\hat{\mathcal{L}}_{t}$ is the time-dependent evolution superoperator) Breuer and Petruccione (2007), the QSL time can be expressed as Deffner and Lutz (2013b); Deffner (2017)

[TABLE]

where $\tau_{p}=\sin^{2}[\mathcal{B}(\hat{\rho}_{0},\hat{\rho}_{\tau})]/E_{\tau}^{p}$. The numerator of $\tau_{p}$ contains the Bures angle, $\mathcal{B}(\hat{\rho}_{0},\hat{\rho}_{\tau})=\arccos(\sqrt{\mathcal{F}(\hat{\rho}_{0},\hat{\rho}_{\tau})})$, between an initial state $\hat{\rho}_{0}$ and a final state $\hat{\rho}_{\tau}$ at time $\tau$ Bures (1969), where $\mathcal{F}=\left[\mathrm{Tr}_{S}\sqrt{\sqrt{\hat{\rho}_{0}}\hat{\rho}_{\tau}\sqrt{\hat{\rho}_{0}}}\right]^{2}$ is the quantum fidelity Jozsa (1994). The denominator $E_{\tau}^{p}=\frac{1}{\tau}\int_{0}^{\tau}dt||\hat{\mathcal{L}}_{t}\hat{\rho}_{t}||_{p}$ is the time average of $||\hat{\mathcal{L}}_{t}\hat{\rho}_{t}||_{p}$ over an actual evolution time interval $\tau$, where $||\hat{A}||_{p}=\left(\sum_{k}\alpha_{k}^{p}\right)^{1/p}$ is the Schatten-$p$-norm of $\hat{A}$ with $\alpha_{k}$ the $k^{th}$ singular value of $\hat{A}$, i.e., the $k^{th}$ eigenvalue of the Hermitian operator $\sqrt{\hat{A}^{\dagger}\hat{A}}$ Bhatia (1997). Thus, Eq. (1) involves Schatten-$p$-norms with $p=1,2$, and $\infty$ (referred to as the familiar trace, Hilbert-Schmidt, and operator norms, respectively). It should be noted that $\tau_{QSL}$ is not a physical time, but rather an intrinsic characteristic timescale associated with a system’s dynamics that satisfies the bound $\tau_{QSL}\leq\tau$ Deffner and Campbell (2017).

Because $\tau_{QSL}$ explicitly depends on the physical time $\tau$, it can be evaluated in two possible ways. The first way, which involves fixing the evolution time $\tau$, has been used to identify memory effects on the speed of quantum evolution Deffner and Lutz (2013b); Xu et al. (2014). In contrast, the second way involves varying the evolution time $\tau$ and studying $\tau_{QSL}$ as a function of $\tau$. This approach has been employed in the study of entanglement-assisted speed-up of quantum evolution Giovannetti et al. (2003); Batle et al. (2005); Borrás et al. (2006); Fröwis (2012). In both cases, the smaller the value of $\tau_{QSL}$ is, the faster the quantum evolution can be realized. Therefore, $\tau_{QSL}$ could potentially be used as a performance metric in the design of quantum-based technologies.

It was previously shown that the operator norm $||\hat{\mathcal{L}}_{t}\hat{\rho}_{t}||_{\infty}$ yields the maximum value of $\tau_{QSL}$ in Eq. (1) Deffner and Lutz (2013b). However, computing the operator norm is far from being a trivial task, as extracting the largest singular value from a time-dependent reduced density matrix is computationally expensive, if at all feasible. Thus far, evaluations of $\tau_{QSL}$ in OQSs have been limited to either exactly solvable models del Campo et al. (2013); Deffner and Lutz (2013b); Xu et al. (2014); Xu and Zhu (2014); Liu et al. (2015); Cheng et al. (2018); Deffner (2017), or have involved approximated forms of time-local master equations with limited applicabilities Mukherjee et al. (2013); Cimmarusti et al. (2015); Funo et al. (2018). Recently, an alternative expression for $\tau_{QSL}$ in Wigner phase space was proposed in Ref. Deffner (2017), which circumvents the task of determining the singular values of a high-dimensional operator. Nevertheless, determining the Wigner function for an arbitrary OQS is very challenging. Intriguing aspects of the QSL time have been revealed in previous studies on small and exactly solvable systems del Campo et al. (2013); Deffner and Lutz (2013b); Xu et al. (2014); Xu and Zhu (2014); Liu et al. (2015); Cheng et al. (2018); Deffner (2017). However, moving forward, it is desirable to develop computationally efficient methods for evaluating $\tau_{QSL}$ in arbitrary OQSs, especially considering the fact that knowledge of the minimal duration of a process is of fundamental importance to virtually all areas of quantum physics.

To address this challenge, we propose herein to compute $\tau_{QSL}$ for multi-level OQSs using a hybrid quantum-classical method Liu and Hanna (2018), which relies on the mixed Wigner-Heisenberg representation Aleksandrov (1981); Gerasimenko (1982); Zhang and Balescu (1988); Kapral and Ciccotti (1999) of the composite quantum dynamics. In contrast to Refs. Deffner (2017); Shanahan et al. (2018), we apply the Wigner transform only to the bath degrees of freedom (yielding a classical-like description of the bath) and retain the operator character of the subsystem degrees of freedom. In doing so, one can treat OQSs with large numbers of bath degrees of freedom and complex bath models, i.e., situations where exact methodologies become computationally intractable.

In this work, without loss of generality, we focus on $\tau_{2}$ in Eq. (1) because the Hilbert-Schmidt norm $||\hat{\mathcal{L}}_{t}\hat{\rho}_{t}||_{2}$ is mathematically less involved than the other measures, and can be expressed as $\sqrt{\mathrm{Tr}_{S}\left(\frac{d}{dt}\hat{\rho}_{t}\right)^{2}}$ Deffner and Lutz (2013b). We note that $\tau_{2}$ can also be obtained as a limit of an improved bound for the QSL time, which relies on an alternative definition of the quantum fidelity between states Ektesabi et al. (2017). The time-dependence of $\frac{d}{dt}\hat{\rho}_{t}$ can be readily simulated using our hybrid quantum-classical method without calculating singular values and without explicitly constructing the superoperator $\hat{\mathcal{L}}_{t}$. Although $\tau_{2}$ does not constitute the “tightest” bound (i.e., the duration of the process always exceeds $\tau_{2}$), it can still yield the same qualitative information as the tightest bound $\tau_{\infty}$. In some cases, $\tau_{2}$ may not even be much smaller than $\tau_{\infty}$, e.g., $\tau_{2}$ is only a factor of $\sqrt{2}$ smaller than $\tau_{\infty}$ for a two-level atom in a photonic crystal cavity Xu et al. (2014).

We apply the hybrid quantum-classical method to two prototype models whose quantum dynamics have been extensively studied over the years: the spin-boson model (SBM) Leggett et al. (1987); Weiss (2012) and the Fenna-Matthews-Olson (FMO) complex model Fenna and Matthews (1975); Adolphs and Renger (2006); Ishizaki and Fleming (2009a). For the SBM, we focus on the influence of the bath temperature and subsystem-bath coupling strength on the behaviour of $\tau_{QSL}$. We also compare our results with those obtained from the non-Markovian Bloch-Redfield equation (NM-BRE) and its Markovian version (M-BRE) Aslangul et al. (1986); Thoss et al. (2001); Boudjada and Segal (2014). For the FMO complex model, we explore the relationship between the QSL time and the efficiency of the excitation energy transfer at a physiological temperature (300 K) and a cryogenic temperature (77 K).

The paper is organized as follows: In section II we present the working expressions for the QSL time and the hybrid quantum-classical method, and then illustrate how to compute the QSL time using this method. In sections III and IV, we apply the hybrid quantum-classical approach to the SBM and FMO complex model. Here, the detailed simulation results and discussions are provided. We summarize our findings in section V. The hybrid quantum-classical equations of motion for the SBM and FMO complex model are provided in the Appendix.

II Working Expressions and Methodology

II.1 Quantum speed limit time expression

As mentioned in the introduction, we will use $\tau_{2}$ in Eq. (1) as a measure of the QSL time $\tau_{QSL}$. Noting the relation between the Bures angle and quantum fidelity, we therefore have

[TABLE]

It should be noted that although the above expression is derived by making use of the time-local master equation, it can be applied to arbitrary OQSs as it only depends on the reduced density matrix and its time derivative. Therefore, one can use any method, approximate or exact, that yields time-dependent information to evaluate $\tau_{QSL}$ in Eq. (2).

In this work, we take the initial condition of the multi-level subsystem to be a pure state, $\hat{\rho}_{0}=|n\rangle\langle n|$, where $\{|n\rangle\}$ spans the subsystem Hilbert space. In this case, the quantum fidelity reduces to Jozsa (1994)

[TABLE]

where $\mathcal{\hat{P}}_{nn}=|n\rangle\langle n|$ is a subsystem projection operator. We further obtain that

[TABLE]

where $L$ is the number of subsystem levels, $\rho_{t}^{mn}=\langle m|\hat{\rho}_{t}|n\rangle$, $|\cdot|$ takes the norm of its argument, and we have used the fact that

[TABLE]

From Eqs. (3) and (4), we see that the QSL time in Eq. (2) is fully determined by ensemble averages of time-dependent projection operators and their time derivatives. Below, we show how to compute these ensemble averages using a hybrid quantum-classical dynamics method.

II.2 The hybrid quantum-classical formalism and the DECIDE implementation

To compute the QSL time in an OQS, one needs a methodology that can accurately capture the reduced dynamics of a quantum subsystem. Such time evolution methods range from numerically exact techniques, such as the quasiadiabatic propagator path integral (QUAPI) Makri and Makarov (1995a, b) approach, to weak-system-bath perturbative methods that are derived systematically from the Nakajima Zwanzig equation Breuer and Petruccione (2007). While path integral-based tools offer an exact numerical solution, they are limited to treating small systems due to their computational costs. As well, converging the dynamics in “difficult” parameter regimes (i.e., strong subsystem-bath coupling, low temperature) becomes increasingly challenging. On the other end, the perturbative Redfield equation is easy to implement and employ, but given its perturbative nature it can only accurately capture weak subsystem-bath coupling effects.

Hybrid quantum-classical methods Tully (1990); Billing (1993); Prezhdo and Kisil (1997); Martens and Fang (1997); Donoso and Martens (1998); Tully (1998); Kapral and Ciccotti (1999); Donoso and Martens (2000); Wan and Schofield (2000); Horenko et al. (2002); Wan and Schofield (2002); MacKernan et al. (2002); Horenko et al. (2004); Roman and Martens (2007); Kim et al. (2008); MacKernan et al. (2008); Bai et al. (2014); Kim and Rhee (2014a, b); Wang et al. (2015); Martens (2016); Wang et al. (2016); Agostini et al. (2016); Subotnik et al. (2016); Kapral (2015), which treat the subsystem of interest quantum mechanically and the bath in a classical-like fashion, are viable alternatives to fully quantum mechanical ones. These methods are particularly useful for modelling quantum dynamical processes occurring in condensed phases. To arrive at a hybrid quantum-classical description of the dynamics, one can first perform a partial Wigner transform Wigner (1932) over the bath degrees of freedom of the quantum Liouville equation, which introduces a phase space description of the bath variables while retaining the operator character of the subsystem degrees of freedom, and then make physically motivated approximations. For example, by linearizing the resulting equation in $\hbar$, one can obtain the quantum-classical Liouville equation Aleksandrov (1981); Gerasimenko (1982); Zhang and Balescu (1988); Kapral and Ciccotti (1999). In this work, we adopt a recently developed hybrid quantum-classical approach Liu and Hanna (2018) that solves the quantum-classical Liouville equation in an efficient manner; namely, solving a system of coupled first-order differential equations (FODEs) for the coordinates of the subsystem and bath. We now describe the approach and its implementation.

Let us consider a generic OQS described by the following Hamiltonian

[TABLE]

Here, $\hat{H}_{S}$ is the subsystem Hamiltonian and $\boldsymbol{\hat{x}}$ collectively refers to the complete set of subsystem projection operators $\mathcal{\hat{P}}_{nm}$ (i.e., the generalized coordinates of the subsystem); $\hat{H}_{B}=\sum_{j=1}^{N}[\hat{P}_{j}^{2}/2+\omega_{j}^{2}\hat{R}_{j}^{2}/2]$ is the Hamiltonian of the bosonic heat bath (containing $N$ harmonic oscillators) characterized by a temperature $T$ (or inverse temperature $\beta\equiv 1/T$), where $\hat{P}_{j}$, $\hat{R}_{j}$, and $\omega_{j}$ are the mass-weighted momentum, position, and frequency of $j$th oscillator, respectively, and $\boldsymbol{\hat{X}}=(\boldsymbol{\hat{R}},\boldsymbol{\hat{P}})$ with $\boldsymbol{\hat{R}}=(\hat{R}_{1},\hat{R}_{2},\ldots,\hat{R}_{N})$ and $\boldsymbol{\hat{P}}=(\hat{P}_{1},\hat{P}_{2},\ldots,\hat{P}_{N})$; and $\hat{H}_{I}$ is the subsystem-bath interaction Hamiltonian. The extension to multiple heat baths is straightforward, as will be seen. In what follows, we set $\hbar=1$ and $k_{B}=1$.

Since the generalized coordinates $\boldsymbol{\hat{x}}$ can be used to construct the reduced density matrix through Eq. (5), their time evolution can be used to determine the QSL time $\tau_{QSL}$. To simulate the time evolution $\boldsymbol{\hat{x}}$, we adopt the mixed Wigner-Heisenberg representation of the composite quantum dynamics Aleksandrov (1981); Gerasimenko (1982); Zhang and Balescu (1988); Kapral and Ciccotti (1999), which involves performing partial Wigner transforms over the bath coordinates of bath-dependent operators $\hat{A}(\boldsymbol{\hat{X}})$, i.e.,

[TABLE]

where $\boldsymbol{X}=(\boldsymbol{R},\boldsymbol{P})$ denotes the bath phase space variables and the subscript $W$ indicates that a partial Wigner transform has been performed. After this transformation, the composite dynamics is equivalently governed by the following Weyl-ordered, mixed Wigner-Heisenberg form of the Hamiltonian

[TABLE]

Weyl-ordering involves replacing product terms such as $\boldsymbol{\hat{x}}\boldsymbol{X}$ in $\hat{H}_{W}$ with the expression $\frac{1}{2}(\boldsymbol{\hat{x}}\boldsymbol{X}+\boldsymbol{X}\boldsymbol{\hat{x}})$. In this work, we simulate the dynamics of $\boldsymbol{\hat{x}}(t)$ and $\boldsymbol{X}(t)$ using the so-called DECIDE (Deterministic Evolution of Coordinates with Initial Decoupled Equations) scheme recently proposed in Ref. Liu and Hanna (2018), which provides an efficient, albeit approximate, way of solving the quantum-classical Liouville equation for $\hat{A}_{W}(\boldsymbol{X},t)$. By assuming a factorized initial state for the composite system (namely, $\hat{\rho}_{tot}(0)=\hat{\rho}_{0}\otimes\hat{\rho}_{B}$ where $\hat{\rho}_{tot}$ and $\hat{\rho}_{B}=\frac{e^{-\beta\hat{H}_{B}}}{\mathrm{Tr}_{B}[e^{-\beta\hat{H}_{B}}]}$ are the density matrices of the total system and bath, respectively), it was shown that the quantum-classical Liouville equation for $\hat{A}_{W}(\boldsymbol{X},t)$ could be unfolded into the following approximate set of coupled equations of motion (EOMs) for the time-dependent coordinates $\boldsymbol{\hat{x}}(t)$ and $\boldsymbol{X}(t)$ (see Ref. Liu and Hanna (2018) for the details of the derivation and the approximations involved),

[TABLE]

where $[\cdot,\cdot]$ and $\{\cdot,\cdot\}$ denote a commutator and Poisson bracket, respectively. To treat the quantum operators and classical variables on an equal footing, we cast Eq. (II.2) in a convenient basis $\{|\alpha\rangle\}=(|\alpha_{1}\rangle,\ldots,|\alpha_{L}\rangle)$ that spans the Hilbert space of the $L$-dimensional quantum subsystem,

[TABLE]

where $D^{\alpha\alpha^{\prime}}\equiv\langle\alpha|D|\alpha^{\prime}\rangle$. For bilinear subsystem-bath interactions, $i\langle\alpha|\left([\hat{H}_{W},\boldsymbol{\hat{x}}]\right)(t)|\alpha^{\prime}\rangle$ becomes a functional of the matrix elements $\{\boldsymbol{x}^{\alpha\alpha^{\prime}}(t)\}$ and $\{(\boldsymbol{\hat{x}}(t)\boldsymbol{X}(t)+\boldsymbol{X}(t)\boldsymbol{\hat{x}}(t))^{\alpha\alpha^{\prime}}\}$, the latter arising from the bilinear interaction in the Weyl-ordered Hamiltonian $\hat{H}_{W}$. On the other hand, $-\langle\alpha|\left(\{\hat{H}_{W},\boldsymbol{X}\}_{a}\right)(t)|\alpha^{\prime}\rangle$ becomes a functional of the matrix elements $\{\boldsymbol{x}^{\alpha\alpha^{\prime}}(t)\}$ and $\{\boldsymbol{X}^{\alpha\alpha^{\prime}}(t)\}$. The notation $\{\boldsymbol{z}^{\alpha\alpha^{\prime}}\}$ denotes a particular set of matrix elements of $\boldsymbol{z}$ in the basis $\{|\alpha\rangle\}$, the contents of which depend on the model under investigation. Also, $(\hat{x}_{l}X_{k})^{\alpha\alpha^{\prime}}=\sum_{\beta}x_{l}^{\alpha\beta}X_{k}^{\beta\alpha^{\prime}}$. It should be noted that, due to the subsystem-bath interactions, the classical coordinates depend on the subsystem operators and $\boldsymbol{X}^{\alpha\alpha^{\prime}}(t)\neq\boldsymbol{X}(t)\delta_{\alpha\alpha^{\prime}}$ at finite times (where $\delta_{\alpha\alpha^{\prime}}$ is the Kronecker delta function). The superscript on $\boldsymbol{X}^{\alpha\alpha^{\prime}}(t)$ acts as a label to distinguish the various $\it{c}$-numbers and their corresponding EOMs at finite times.

Equation (II.2) constitutes a set of coupled FODEs for the c-numbers ($\boldsymbol{x}^{\{\alpha\alpha^{\prime}\}}(t),\boldsymbol{X}^{\{\alpha\alpha^{\prime}\}}(t)$), where $\{\alpha\alpha^{\prime}\}$ denotes all the combinations of basis indices. The maximum number of coupled FODEs is $L^{2}(L^{2}-1+2N)$: the subsystem is described by $L^{2}-1$ projection operators (because the identity operator is excluded), the $N$ harmonic oscillators are described by $N$ displacements and $N$ momenta, and each of the resulting $(L^{2}-1+2N)$ coordinates has $L^{2}$ equations. However, one could reduce the number of FODEs if the subsystem has some symmetry. It should be noted that DECIDE is not limited to bilinear interactions, provided that the interaction can be decomposed into a finite number of terms involving matrix elements of coordinates.

The domain of applicability of DECIDE warrants some comments. To arrive at Eq. (II.2), one must truncate the corresponding quantum Heisenberg equations for the coordinates (see Ref. Liu and Hanna (2018) for details). In doing so, one neglects higher order terms that contain derivatives of the time-dependent coordinates with respect to the initial bath coordinates. As a result, one can underestimate the back-action from the bath onto the subsystem. If the subsystem dynamics is highly non-Markovian, these higher order terms are important and DECIDE will yield inaccurate results in the long-time limit. We note that strong memory effects can be induced by strong subsystem-bath couplings, slow baths with characteristic timescales longer than that of the subsystem, and very low temperatures. Therefore, DECIDE should be used with caution in these regimes.

II.3 Evaluation of quantum speed limit time

The ensemble averages of time-dependent projection operators and their time derivatives found in the QSL time expression in Eq. (2) have the following forms in the mixed Wigner-Heisenberg representation Sergi et al. (2003)

[TABLE]

where $\rho_{B,W}(0)$ is the partially Wigner-transformed thermal equilibrium distribution Imre et al. (1967)

[TABLE]

satisfying $\int d\boldsymbol{X}(0)\rho_{B,W}(\boldsymbol{X}(0))=1$.

Equation (11) suggests a trajectory-based molecular dynamics (MD) approach to compute the QSL time $\tau_{QSL}$: One generates a swarm of independent classical-like trajectories starting from different $\boldsymbol{X}(0)$ sampled from $\rho_{B,W}(0)$, and the same initial values of the matrix elements of the projection operators. Each trajectory of $\mathcal{P}_{nm}^{\alpha\alpha^{\prime}}(t)$ is obtained by integrating the $L^{2}(L^{2}-1+2N)$ coupled FODEs in Eq. (II.2) using the standard fourth-order Runge-Kutta scheme Dormand and Prince (1980). Inserting the time-dependent coordinates back into Eq. (II.2) yields the corresponding time derivatives. Averaging the $\mathcal{P}_{nm}^{\alpha\alpha^{\prime}}(t)$ and their time derivatives over the ensemble of trajectories yields the required ensemble averages found in the QSL time. Finally, Simpson’s rule Süli and Mayers (2003) is used to perform the numerical integration in the denominator of Eq. (2). To illustrate our methodology, we study the SBM and FMO complex model in the following sections.

III The spin-boson model

The behaviour of the QSL time has been explored in several exactly-solvable OQS models, including the damped Jaynes-Cummings model and the pure dephasing model Deffner and Lutz (2013b); Cheng et al. (2018). Herein, using the DECIDE hybrid quantum-classical method, we study the behaviour of the QSL time in the spin-boson model Leggett et al. (1987), which exhibits a rich dynamics but lacks a closed analytic solution Weiss (2012).

The SBM consists of a two-level spin in contact with a bosonic heat bath and has the following Hamiltonian Leggett et al. (1987); Weiss (2012)

[TABLE]

where $\hat{\sigma}_{x/z}$ are the Pauli spin matrices, $\Delta$ is the tunnelling frequency between the two spin states, $\omega_{j}$ is the frequency of the $j$th harmonic oscillator, and $C_{j}$ is the coupling coefficient between the spin and the $j$th harmonic oscillator. To characterize the influence of the heat bath on the two-level subsystem, we employ an Ohmic spectral density with an exponential cutoff, namely $J(\omega)=\frac{\xi}{2}\pi\omega e^{-\omega/\omega_{c}}$, where the Kondo parameter $\xi$ characterizes the subsystem-bath coupling strength and $\omega_{c}$ is the cutoff frequency. The corresponding Weyl-ordered Hamiltonian in the mixed Wigner-Heisenberg representation is given by

[TABLE]

III.1 Quantum speed limit time

To evaluate the QSL time in the SBM, we set the initial spin state as the spin-up state, such that $\hat{\rho}_{0}=|+\rangle\langle+|$, where $|\pm\rangle$ are the eigenstates of $\hat{\sigma}_{z}$. In this case, the quantum fidelity in Eq. (3) becomes

[TABLE]

where $B_{m}(\tau)\equiv\langle\hat{\sigma}_{m}(\tau)\rangle$ ($m=x,y,z$). Using the fact that the reduced density matrix of a two-level system can be expressed as $\hat{\rho}_{t}=\frac{1}{2}(\mathcal{I}+\sum_{m=x,y,z}B_{m}(t)\hat{\sigma}_{m})$, with $\mathcal{I}$ the $2\times 2$ identity matrix, Eq. (4) may be rewritten as

[TABLE]

Thus, the QSL time is fully determined by the ensemble averages of the Pauli matrices and their corresponding time derivatives.

It is interesting to note that one may derive an analytical expression for $\tau_{QSL}$ in the isolated limit (i.e., $\hat{H}_{I}=0$) where the spin dynamics is governed by the Bloch equations Weiss (2012): $\frac{d}{dt}B_{x}(t)~{}=~{}0,~{}~{}\frac{d}{dt}B_{y}(t)~{}=~{}2\Delta B_{z}(t),~{}~{}\frac{d}{dt}B_{z}(t)~{}=~{}-2\Delta B_{y}(t)$. Choosing the spin-up state as the initial state, the solution turns out to be $B_{x}(t)=0$, $B_{y}(t)=\sin(2\Delta t)$ and $B_{z}(t)=\cos(2\Delta t)$. Therefore, the QSL time for an isolated two-level subsystem is

[TABLE]

This will be used as a reference point to assess the effect of the subsystem-bath interaction on the QSL time.

III.2 Computational details

We now describe the details of our DECIDE simulations. In addition, to pinpoint the effects of stronger subsystem-bath coupling and non-Markovianity on the QSL times, we compare our DECIDE results (which capture these effects to a certain extent) to those of the perturbative Bloch-Redfield master equation with and without the Markov approximation.

III.2.1 DECIDE simulations

The generalized coordinates for the subsystem are taken to be the Pauli matrices, viz., $\boldsymbol{\hat{x}}=(\hat{\sigma}_{x},\hat{\sigma}_{y},\hat{\sigma}_{z})$. The DECIDE EOMs for the subsystem and bath coordinates are derived from Eq. (II.2) using $\hat{H}_{W}$ from Eq. (14) (the resulting EOMs are listed in Eq. (.1) in the Appendix). To compute the QSL time, we must calculate the ensemble averages $B_{m}(t)$ and their time derivatives (taking $\{|\alpha\rangle\}=\{|+\rangle,|-\rangle\}$ for convenience and the initial subsystem state to be $\hat{\rho}_{0}=|+\rangle\langle+|$) according to,

[TABLE]

We evaluate the right-hand-sides of the above equations by averaging over a swarm of independent classical-like trajectories, with each trajectory starting from different values of the bath coordinates and the same values of the Pauli matrix elements. More specifically, the initial values of the bath coordinates are $\boldsymbol{X}^{\alpha\alpha^{\prime}}(0)=\boldsymbol{X}(0)\delta_{\alpha\alpha^{\prime}}$ (due to the factorized initial state), with $\boldsymbol{X}(0)$ sampled from $\rho_{B,W}(0)$ in Eq. (12). The Ohmic bath spectral density with an exponential cutoff is discretized as Thompson and Makri (1999); Wang et al. (2001)

[TABLE]

where $j$ runs from 1 to $N$, and $\omega_{0}=\frac{\omega_{c}}{N}(1-e^{-\omega_{max}/\omega_{c}})$ with $\omega_{max}$ the maximum frequency of the bath oscillators. It should be noted that although we employ an Ohmic spectral density here, DECIDE, just like any other hybrid quantum-classical dynamics method, can handle arbitrary bath spectral densities. Noting the forms of the Pauli matrices in the $\{|+\rangle,|-\rangle\}$ basis, it is easy to show that the non-vanishing initial values of the Pauli matrix elements are $\sigma_{x}^{-+}(0)=\sigma_{x}^{+-}(0)=1,\sigma_{y}^{-+}(0)=i,\sigma_{y}^{+-}(0)=-i,\sigma_{z}^{--}(0)=-1$, and $\sigma_{z}^{++}(0)=1$. Starting from the aforementioned initial conditions, we then numerically integrate Eq. (.1) using the fourth-order Runge-Kutta method Dormand and Prince (1980), which results in trajectories of $\sigma_{m}^{++}(t)$ and their time derivatives. Finally, averaging $\sigma_{m}^{++}(t)$ and their time derivatives over the ensemble of trajectories yields the required ensemble averages for constructing the QSL time, $\tau_{QSL}$.

III.2.2 Bloch-Redfield master equation simulations

To understand the impacts of strong subsystem-bath interactions and Markovian dynamics, we compare the results of our DECIDE simulations to those of the perturbative non-Markovian Bloch-Redfield equation (NM-BRE) and Markovian Bloch-Redfield equation (M-BRE). Perturbative methods offer a simple means of simulating the reduced dynamics of two-level quantum systems Breuer and Petruccione (2007). By treating the subsystem-bath interaction, $\hat{H}_{I}$, as a perturbation to the lowest nontrivial order, one can obtain the NM-BRE within the Born approximation Aslangul et al. (1986); Thoss et al. (2001); Boudjada and Segal (2014),

[TABLE]

The kernels of these integro-differential equations are $\Gamma_{x}(t)=-\sin(2\Delta t)M_{2}(t)$, $\Gamma_{xx}(t)=\cos(2\Delta t)M_{1}(t)$, and $\Gamma_{yy}(t)=M_{1}(t)$, with $M_{1}(t)$ and $M_{2}(t)$ satisfying the following relation,

[TABLE]

Next, introducing the Markov approximation, one arrives at the simpler M-BRE, which involves the following set of time-local equations Thoss et al. (2001),

[TABLE]

where the time-dependent coefficients are given by $G_{x}(t)=-\int_{0}^{t}dt^{\prime}\sin(2\Delta t^{\prime})M_{2}(t^{\prime})$, $G_{xx}(t)=G_{yy}(t)=\int_{0}^{t}dt^{\prime}\cos(2\Delta t^{\prime})M_{1}(t^{\prime})$, and $G_{yz}(t)=-\int_{0}^{t}dt^{\prime}\sin(2\Delta t^{\prime})M_{1}(t^{\prime})$. After integrating Eqs. (III.2.2) and (III.2.2) using the fourth-order Runge-Kutta scheme Dormand and Prince (1980), we compute the QSL time from the $B_{m}(t)$’s and their time derivatives.

III.3 Numerical results

Our goal is to investigate the interplay of subsystem-bath coupling strength and bath temperature on the behaviour of the QSL time. Noting the regime of validity of DECIDE, here we focus on cases with $\Delta<\omega_{c}$ and $T>\Delta$, while we vary the coupling strength $\xi$.

To benchmark the performances of the various methods for simulating the reduced dynamics of the two-level system, we first compare in Fig. 1 the spin polarization $\langle\hat{\sigma}_{z}(\tau)\rangle$ as obtained from NM-BRE, M-BRE, and DECIDE, with those calculated using the numerically exact QUAPI method Makri and Makarov (1995a, b).

From the comparison, it is evident that M-BRE can only be applied in the very weak coupling regime, as even for $\xi=0.001$ we observe deviations from the numerically exact QUAPI results. By taking the non-Markovianity into account, NM-BRE improves upon M-BRE in the weak coupling regime as seen in panels (a)-(b) of Fig. 1. While both M-BRE and NM-BRE cannot be applied in the strong coupling regime given their perturbative nature, as seen in panels (c)-(d) of Fig. 1, it is important to note that they predict distinct dynamical behaviours in the strong coupling regime, viz., M-BRE fails to capture the bath-induced spin polarization, while NM-BRE overestimates the spin polarization (see panel (d) of Fig. 1). We will return to this observation when interpreting the behaviour of the QSL time in Fig. 2. In contrast, DECIDE works well across all of the coupling regimes, demonstrating its reliability in simulating the reduced dynamics and hence the QSL time.

We now turn to the study of the QSL time. In Fig. 2, we depict $\tau_{QSL}$ with a fixed evolution time of $\tau=1$ as a function of the subsystem-bath coupling strength $\xi$. In addition, we consider two different spin tunnelling frequencies $\Delta=0.2$ (Fig. 2 (a)) and $\Delta=0.6$ (Fig. 2 (b)).

Since the spin-boson system starts from a factorized initial state, NM-BRE and M-BRE can qualitatively capture the reduced dynamics at a very short time scale across the whole coupling regime as seen in Fig. (1). Therefore, it is not surprising that we observe similar turnover behaviours of $\tau_{QSL}$ as a function of the coupling strength $\xi$ using all three methods. However, we should emphasize that only the predictions of DECIDE are reliable in the strong coupling regime, given the results in Fig. 1. When we compare the results of NM-BRE and M-BRE in the weak coupling regime, we find that the former predicts a smaller value of $\tau_{QSL}$ than the latter, for both values of the temperatures and spin tunnelling frequencies, which is a direct demonstration of non-Markovianity-induced speed-up Deffner and Lutz (2013b). We also note the perfect agreement between the DECIDE and NM-BRE results in the weak coupling regime, pointing to the importance of non-Markovianity even at high temperatures and weak subsystem-bath couplings, where it is often believed that non-Markovianity is minimal.

When we increase the spin tunnelling frequency from $\Delta=0.2$ to $\Delta=0.6$, we see that DECIDE predicts a sharper decrease in the QSL time in the strong coupling regime (cf. Figs. 2 (a) and (b)), and this leads to $\tau_{QSL}<\tau_{QSL}^{iso}$. This suggests that one could speed up the quantum evolution of an open quantum system by entering into the strong coupling regime with large tunnelling frequencies. This result is also a direct consequence of the non-Markovianity-induced speed-up Deffner and Lutz (2013b), since a larger spin tunnelling frequency corresponds to a larger ratio of $\Delta/\omega_{c}$ and consequently a stronger non-Markovian effect. As for the role of temperature in the QSL time, in Fig. 2 we observe that both M-BRE and DECIDE predict an increase in $\tau_{QSL}$ with increasing temperature for all couplings, while NM-BRE predicts an inverse temperature dependence in the very strong coupling regime. In the strong coupling regime, the spin is more polarized at the lower temperature, which implies a shorter QSL time $\tau_{QSL}$ at lower temperatures. Therefore, NM-BRE fails to capture the correct temperature dependence of $\tau_{QSL}$ in the strong coupling regime, while M-BRE does. This result demonstrates that since the perturbative and Markov approximations are linked, it is sometimes more consistent to enforce them together (as in the M-BRE) than separately.

Qualitatively speaking, the turnover behaviour of the QSL time at higher coupling values, which is reminiscent of the turnover behaviour of the thermal conductance observed in the SBM Boudjada and Segal (2014); Liu et al. (2017), results from the rapid energy exchange processes between the subsystem and bath in the strong coupling regime. We can gain quantitative insight into the reduction of the QSL time by looking at the time dependences of the spin polarization, $\langle\hat{\sigma}_{z}(\tau)\rangle$, and $\sqrt{\mathrm{Tr}_{S}(\dot{\rho}_{t})^{2}}$, as calculated by DECIDE, in the strong coupling regime (see Fig. 3). We find that the stronger the coupling strength, the larger is the spin polarization. This monotonically increasing behaviour in the spin polarization implies a monotonically decreasing trend in the numerator of Eq. (2) (see also Eq. (15) in the strong coupling regime. On the other hand, the integrand in the denominator of Eq. (2), $\sqrt{\mathrm{Tr}_{S}(\dot{\rho}_{t})^{2}}$, is almost independent of coupling strength in the strong coupling regime.

In Fig. 4, we present results of the QSL time at different temperatures as a function of the coupling strength, but now using longer evolution times $\tau$. We find that an increase in the evolution time from $\tau=1$ to $\tau=10$ and $\tau=20$, does not induce qualitative changes in the predictions of DECIDE compared to what we found in Fig. 2. However, $\tau_{QSL}$ obtained from M-BRE now exhibits a monotonically increasing trend, as M-BRE fails to accurately describe the reduced dynamics at long times in the strong coupling regime. Although NM-BRE can still predict a turnover behaviour, the temperature dependence of $\tau_{QSL}$ in the strong coupling regime, as determined by NM-BRE, is inaccurate.

IV The Fenna-Matthews-Olson complex

There has been much interest in studying the dynamics of excitation energy transfer (EET) in light-harvesting complexes, with a focus on the question of whether or not electronic quantum coherences assist the highly efficient transfer of energy in the photosynthetic process. In particular, the FMO complex has gained much attention, with early 2D electronic spectroscopy experiments suggesting the existence of long-lived coherences at 77 K Engel et al. (2007), and more recent studies disagreeing with the interpretation of the early experimental results Duan et al. (2017). Apart from the debate over the existence of these long-lived coherences and their potential impact on the energy transfer efficiency, there has been a multitude of theoretical and computational studies on the FMO complex model aimed at benchmarking new methodologies and questioning the role of coherences in quantum dynamics Ishizaki and Fleming (2009a, b); Nalbach and Thorwart (2010); Ke and Zhao (2016). Herein, we demonstrate the computation of QSL times for the FMO model and explore the relationship between these times and the EET efficiency at different temperatures.

IV.1 Model

The photosynthetic FMO complex can be described by the following Frenkel exciton Hamiltonian in the single-excitation subspace Ishizaki and Fleming (2009a)

[TABLE]

In the above Hamiltonian, $\mathcal{\hat{P}}_{nm}=|n\rangle\langle m|$, where the basis state $|n\rangle$ corresponds to the $n$th chromophoric site in its electronically excited state and the remaining sites in their electronic ground state; $E_{n}$ denotes the excited state energy corresponding to state $|n\rangle$; and $V_{mn}$ is the excitonic coupling strength between the $n$th and $m$th sites. Each site is coupled to an independent heat bath containing $M$ harmonic oscillators at temperature $T$. We assume that the spectral densities of all the baths are equivalent, and can be characterized by Debye-Drude functions $J(\omega)=2\lambda_{D}\frac{\omega\tau_{c}}{1+\omega^{2}\tau_{c}^{2}}$ Ishizaki and Fleming (2009a), where $\lambda_{D}$ is the bath reorganization energy and $\tau_{c}$ is the characteristic time. We adopt the values $\lambda_{D}=35~{}\mathrm{cm}^{-1}$ and $\tau_{c}=50~{}\mathrm{fs}$, as they yield excellent agreement between the experimental data and numerical simulations Read et al. (2008).

As an illustration, we consider the apo-FMO complex, which includes seven bacteriochlorophyll (BChl) pigment-proteins per subunit; the conventional numbering of the BChls has been used here. The values of the site energies $E_{n}$ and excitonic coupling strengths $V_{mn}$ for a subunit of the apo-FMO complex Adolphs and Renger (2006); Ishizaki and Fleming (2009a) are listed in Table 1.

The corresponding Weyl-ordered Hamiltonian in the mixed Wigner-Heisenberg representation takes the form

[TABLE]

where $V_{nn}=E_{n}+\sum_{j=1}^{M}\frac{C_{n,j}^{2}}{2\omega_{n,j}^{2}}$. For convenience and without loss of generality, we assume that the single excitation is initially located at site 1, such that $\hat{\rho}_{0}=|1\rangle\langle 1|$. As a result, the quantum fidelity in Eq. (3) reduces to the simple form,

[TABLE]

Equation (4) remains the same with $L=7$.

It was previously found that an initial excitation at BChl 1 rapidly transfers to BChls 3 and 4 (collectively known as the target region because they make contact with the reaction center complex) according to the following EET pathway Ishizaki and Fleming (2009a),

[TABLE]

where the double arrow implies that the excitation energy equilibrates between BChls 3 and 4 after BChl 3 is populated.

IV.2 Simulation details

It is well known that standard Redfield-type methods fail in describing EET even for dimers Ishizaki and Fleming (2009b); Nalbach and Thorwart (2010). Therefore, we only use DECIDE to simulate the FMO model. The DECIDE EOMs for the subsystem and bath coordinates are listed in Eq. (.2) in the Appendix. The only non-zero element of the initial subsystem density matrix is $\rho_{0}^{11}=1$. Based on Eq. (11) and taking $\{|\alpha\rangle\}=\{|1\rangle,|2\rangle,\ldots,|7\rangle\}$, the ensemble averages of the projection operators and their time derivatives, which are involved in the QSL time, are given by

[TABLE]

To compute these ensemble averages, we generate a swarm of independent classical-like trajectories, with each trajectory starting from different values of the bath coordinates and the same values of the Pauli matrix elements. More specifically, the initial values of the bath coordinates are $\boldsymbol{X}^{\alpha\alpha^{\prime}}(0)=\boldsymbol{X}(0)\delta_{\alpha\alpha^{\prime}}$ (due to the factorized initial state), with $\boldsymbol{X}(0)$ sampled from $\rho_{B,W}(0)$, which is now a product of seven partially Wigner-transformed Gaussian distributions (each given by Eq. (12)). The Debye-Drude spectral density is discretized as Wang et al. (2001)

[TABLE]

where $\omega_{n,j}=\tan\left[j\arctan(\omega_{max}\tau_{c})/M\right]/\tau_{c}$. In the basis $\{|\alpha\rangle\}=\{|1\rangle,|2\rangle,\ldots,|7\rangle\}$, the initial nonzero values of the subsystem projection operator matrix elements, $\mathcal{P}_{nm}^{nm}(0)$, are

[TABLE]

Starting from the aforementioned initial conditions, we then numerically integrate Eq. (.2) using the fourth-order Runge-Kutta scheme Dormand and Prince (1980), which results in trajectories of $\mathcal{P}_{nm}^{11}(t)$ and their time derivatives. Finally, averaging $\mathcal{P}_{nm}^{11}(t)$ and their time derivatives over the ensemble of trajectories yields the required ensemble averages for constructing the QSL time, $\tau_{QSL}$.

IV.3 Numerical results

The validity of the DECIDE method for treating the FMO complex model was established in Ref. Liu and Hanna (2018). In Fig. 3 of Ref. Liu and Hanna (2018), the time-dependent populations for the first four BChl pigment-proteins (the populations of the remaining pigments remain very small) obtained using DECIDE are compared with those obtained using the numerically exact forward-backward stochastic Schrödinger equation (FB-SSE) Ke and Zhao (2016). As can be seen, DECIDE works well in describing the EET at both the physiological temperature (300 K) and the cryogenic temperature (77 K).

As the parameters of the FMO complex model were determined by numerically fitting the experimental data, we only vary the evolution time $\tau$ to study the QSL time, unlike in our calculations for the SBM. We recall that the QSL time can be analyzed either by keeping $\tau$ fixed, or by studying the behaviour of the QSL time as a function of $\tau$. In Fig. 5, we present results for $\tau_{QSL}$ as a function of $\tau$ for both a physiological temperature (300 K) and a cryogenic temperature (77 K).

At $\tau=0$, $\tau_{QSL}=0$ due to a unit quantum fidelity. As time progresses, we first observe oscillations in $\tau_{QSL}$ at both temperatures, with the oscillations at 77 K being more pronounced. At later times, after the oscillations have decayed, $\tau_{QSL}$ is seen to grow almost linearly with $\tau$. Apart from the very short time regime, the FMO complex at 77 K has a shorter QSL time, or equivalently, a faster quantum evolution speed. Since the QSL time sets the bound on the minimal evolution time for an initial quantum state to reach a target state, a shorter QSL time should correspond to a faster EET process between BChl 1 and the target region. The fact that $\tau_{QSL}\ll\tau$ here is noteworthy and warrants further investigation.

To understand the significance of a faster quantum evolution speed at 77 K, we recall that BChls 3 and 4 define the target region in contact with the reaction center complex. The total population accumulating in time at sites 3 and 4 at both 77 K and 300 K is shown in the inset of Fig. 5. As can be seen from the inset, the results for the two temperatures are very similar up to $\sim 0.6$ ps, after which the total population at 77 K begins to exceed that at 300 K. By $\tau=1$ ps, the total population at 77 K is slightly larger than that at 300 K (0.4 vs. 0.35, respectively). On the other hand, the QSL times at these two temperatures begin to show deviations at a very short time and differ by about a factor of two at $\tau=1$ ps. These differences are a reflection of the fact that the populations only convey a portion of the dynamical information contained in the QSL time (as the QSL time also contains information about the coherences). Now, since only the excitation energy that has been transferred into the target region can be utilized in the reaction center, the total population accumulation at BChls 3 and 4 essentially determines the efficiency of the EET process. Therefore, by comparing our QSL times to the total populations of BChls 3 and 4 at the two different temperatures, our results suggest that a faster quantum evolution speed could lead to a more efficient EET in the long-time limit. For this reason, studying the behaviours of QSL times may allow one to discover useful design principles for creating efficient artificial light-harvesting devices.

V Summary

Hybrid quantum-classical dynamics methods have been primarily developed to simulate quantum processes in condensed phase systems, such as chemical reaction dynamics and vibrational energy relaxation. The main goal of our paper was to offer a hybrid quantum-classical method as a powerful tool for studying performance bounds in open quantum systems. Specifically, we simulated the quantum speed limit in multi-level open quantum systems, using a hybrid quantum-classical method that is based on the mixed Wigner-Heisenberg representation of the composite quantum dynamics. Within this approach, the information over the QSL time is encoded into a set of coupled first-order differential equations for the quantum and classical degrees of freedom. The flexibility inherent to this method allowed us to explore systems that go beyond very simple models and parameters regimes that could not be examined with low-order perturbative quantum master equations.

We studied QSL times in two models. First, we considered the spin-boson model, a prototype open quantum system. The effects of bath temperature and subsystem-bath coupling strength on the behaviour of the QSL time were analyzed. In particular, we found that the QSL time exhibits a turnover behaviour as a function of the subsystem-bath coupling strength, which we attributed to the strong dissipation effect in the strong coupling regime. Under certain circumstances, this turnover could even lead to a speed-up of the quantum evolution in the strong coupling regime as compared to in the isolated limit. By comparing the results from the hybrid quantum-classical method with those obtained from NM-BRE and M-BRE, we concluded that perturbative methods should not be used to calculate QSL times beyond their strict regimes of validity.

As a second example, we studied the QSL time in the Fenna-Matthews-Olson complex. By varying the evolution time, we found that the quantum evolution speed increases as we reduce the temperature. We suggested that this faster quantum evolution speed could lead to a faster excitation energy transfer (EET) to the reaction center complex.

Future studies will aim at gaining a better understanding of the role of the QSL time in open quantum systems using the hybrid quantum-classical method. We also anticipate the development of potential quantum control protocols for EET applications by identifying key characteristics of the QSL time in a given EET process.

Acknowledgements.

J. Liu and G. Hanna acknowledge support from the Natural Sciences and Engineering Research Council of Canada (NSERC). D. Segal acknowledges support from an NSERC Discovery Grant and the Canada Research Chair program.

Appendix: DECIDE equations of motion

The general EOMs in the DECIDE hybrid mixed quantum-classical method are given in Eq. (II.2). Here, we provide the specific DECIDE equations for the spin-boson and FMO complex models.

.1 Equations for the spin-boson model

From Eq. (II.2), one can write down the EOMs for the quantum and classical coordinates using the form of the spin-boson Hamiltonian in Eq. (14). This results in the following set of coupled FODEs for the matrix elements of the spin Pauli matrices and bath oscillator coordinates (taking $\{|\alpha\rangle\}=\{|+\rangle,|-\rangle\}$),

[TABLE]

where the terms of the form $[R_{j}(t)\hat{\sigma}_{m}(t)]^{\alpha\alpha^{\prime}}$ are evaluated as $\sum_{\beta}R_{j}^{\alpha\beta}(t)\sigma_{m}^{\beta\alpha^{\prime}}(t)$. As can be seen, the above set of EOMs consists of $4\times(3+2N)$ coupled FODEs for the matrix elements $(\sigma_{x}^{\{\alpha\alpha^{\prime}\}},\sigma_{y}^{\{\alpha\alpha^{\prime}\}},\sigma_{z}^{\{\alpha\alpha^{\prime}\}},\boldsymbol{X}^{\{\alpha\alpha^{\prime}\}})$, where $N$ is the number of bath oscillators.

.2 Equations for the FMO complex model

From Eq. (II.2), one can write down the EOMs for the quantum and classical coordinates using the form of the FMO complex Hamiltonian in Eq. (24). This results in the following set of coupled FODEs for the matrix elements of the chromophoric site projectors $\{\hat{\mathcal{P}}_{nm}\}$ (with $m,n\in\{1,2,\ldots,7\}$) and bath oscillator coordinates (taking $\{|\alpha\rangle\}=\{|1\rangle,|2\rangle,\cdots,|7\rangle\}$),

[TABLE]

where the terms of the form $[R_{n,j}(t)\hat{\mathcal{P}}_{nm}(t)]^{\alpha\alpha^{\prime}}$ are evaluated as $\sum_{\beta}R_{n,j}^{\alpha\beta}(t)\mathcal{P}_{nm}^{\beta\alpha^{\prime}}(t)$. Given the completeness condition $\sum_{n=1}^{7}\hat{\mathcal{P}}_{nn}=1$, there are $49\times(48+14M)$ coupled FODEs for the subsystem and bath matrix elements, where $M$ is the number of bath oscillators coupled to each site.