Time-dependent renormalized quantum master equation

The study of quantum dissipative systems is necessary for understanding the elementary processes of life and for the development and control of quantum devices [1, 2, 3]. Many theoretical developments have been reported for such systems, and there are two main investigative approaches: (1) nonperturbative analysis, where the system dynamics are calculated using computational resources [4, 5, 6, 7]; and (2) techniques based on analytical approximation, such as, a perturbative method with polaron representation or variational principles [8, 9, 10, 11, 12, 13, 14, 15, 16]. These analytical techniques limit the applicability of quantum dissipative systems, although they have the advantage that they require considerably lesser amount of computational resources.

To consider the dynamics of a reference system within a dissipative system, we must only construct the formalism for that system. The projection operator technique is an easy means of deriving the closed equation of the relevant system, which is called the “quantum master equation (QME)”. After the closed equation is derived, we can employ the appropriate approximation for the system. For example, when the total system consists of two electronic exciton states and nuclear vibrations, there are two limiting cases. First, when the excitonic coupling is much smaller than the other parameters, the exciton state is localized at one site but can be transferred to another site. The transfer rate can be expressed through the Förster formula [17] by treating the excitonic coupling as a perturbation. Second, in the opposite limiting case, the exciton–phonon interaction is treated as the perturbative term, and Redfield theory [18] describes the reaction dynamics between each delocalized exciton state.

The projection operator method developed by Nakajima and by Zwanzig [19, 20] reduces the information of the total system to facilitate the derivation of the QME. However, since the conventional projection operator simply reduces the information of the irrelevant system, its dynamics are not considered by this method. Willis and Picard proposed the time-dependent projection operator[21, 22], which also incorporates the dynamics of the irrelevant system. Their derived QME consists of parallel equations for the relevant and irrelevant systems and has nonlinear terms. The merit of this approach is that it is applicable when the irrelevant system is in a non-equilibrium state. Furthermore, Linden and May reported the derivation of a QME with nonlinear terms achieved using a time-dependent projection operator technique [23]. In their formalism, the nonlinear QME becomes a linear one in the Markovian limit. Recently, Degenfeld-Schonburg and Hartmann developed the self-consistent Mori projector technique[24], whereby equations of reduced density operators are analytically derived for a small number of closed parallel nonlinear equations. This indicates that the nonlinear terms in the nonlinear equation are important for analyzing the dynamics of the reduced density matrix in a many-particle system. However, the properties of the nonlinear term of the nonlinear QME in these theories remain unclear.

In this paper, we first utilize time-dependent projection operator techniques to derive a conventional nonlinear QME, in the manner of Shibata and Hashitsume[25]. Next, we apply time-dependent renormalization to this derived equation, drastically approximate the new equation to an appropriate perturbation order, and apply it to the weak electronic coupling model. Furthermore, we present the numerical results obtained with this new equation and compare them with those produced by the conventional theory and by a numerically exact calculation. Finally, we discuss the various properties of the new theory.

We start by defining the Hamiltonian and its Liouville operator. The total system is constructed as the sum of “relevant” (matter) and “irrelevant” (field) parts that can interact with each other. Hence, the total Hamiltonian is expressed as

where, for example, $H_{m}$ is the electronic Hamiltonian for the relevant (matter) system, $H_{f}$ is the irrelevant (field) Hamiltonian describing the environment nuclear vibration (phonon field), and $H_{m-f}$ is the interaction Hamiltonian between the relevant and irrelevant systems. Based on the Hamiltonians, we define each eigenvector as $|\phi_{i}^{f}\rangle$ ($|\phi_{j}^{m}\rangle$) of the $i$th ($j$th) quantum number for the relevant (irrelevant) Hamiltonian.

For perturbative approximation, we assume that the total Hamiltonian of Eq. (1) can be divided into two terms as follows:

where $H^{0}$ is the nonperturbative Hamiltonian, and $H^{\prime}$ is the perturbative Hamiltonian. The definition of the Liouville operator is completely analogous. Based on the Liouville operator, the Liouville–von Neumann equation in the interaction representation can be written as

where $W_{I}(t)$ is the density operator in the interaction representation and $L^{\prime}_{I}(t)$ is the interaction representation of $L^{\prime}$, defined as $U_{0}^{{\dagger}}(t)L^{\prime}\equiv e^{iL^{0}t/\hbar}L^{\prime}$.

To derive the QME, the time-dependent projection operator[21] is defined as

where $\mbox{Tr}_{f}$ and $\mbox{Tr}_{m}$ are the exclusive trace operators for the irrelevant and relevant systems, respectively (see Fig. 1),

and Tr is defined as $\mbox{Tr}_{f}\mbox{Tr}_{m}$. Finally, $r_{I}(t)$ and $R_{I}(t)$ are the reduced density operators for the relevant and irrelevant systems, respectively, defined as

The time-dependent projection operator satisfies the following relations: $[P(t),\partial/\partial t]\equiv-\dot{P}(t)$, $\dot{P}(t)W_{I}(t)=0$, $P(t_{1})P(t_{2})\equiv P(t_{1})$, $Q(t_{1})Q(t_{2})\equiv Q(t_{2})$, and $P(t_{1})Q(t_{2})\equiv 0$, where we define $Q(t)=1-P(t)$.

As an analytical application, we used the following spin-boson Hamiltonian for a two-electronic state system:

where $E_{i}$ is the energy of the electronic excited state $|i\rangle$ in the $i$th site, and $V_{ij}$ denotes the exciton coupling strength between the $i$th and $j$th sites. In addition to the electronic states, we considered the phonon bath for nuclear vibration only: $b_{k}^{{\dagger}}$ ($b_{k}$) is the creation (annihilation) operator for the $k$th phonon mode, and its frequency is $\omega_{k}$. Finally, $g_{ik}$ is the exciton–phonon coupling strength for the $i$th excited state and the $k$th phonon mode.

In the following section, we numerically analyze the property of Eq. (33) under secular approximation. The expression of the elements of the reduced density matrix $r_{\mathcal{I}}(t)$ is given in the Appendix. Hence, according to Eq. (57), we ultimately obtain the nonlinear QME for the weak electronic coupling limit as follows:

where $k_{da}$ and $K_{da}$ are respectively expressed as

where the bracket $\langle\cdots\rangle_{B}$ is defined as $\mbox{Tr}_{f}[\cdots\rho_{B}]$, $\rho_{B}$ is $e^{-\beta H_{B}}/\mbox{Tr}_{f}[e^{-\beta H_{B}}]$, and $H_{B}$ is bath Hamiltonian defined as $\sum_{k}\hbar\omega_{k}b_{k}^{{\dagger}}b_{k}$.

For Eqs. (40) and (41), we assume that the vibrational relaxation is very fast, and that the distribution expressed by $\rho_{B}$ may represent a thermal equilibrium state in which the excitation almost localizes at each electronic state as $e^{-iH_{i}t/\hbar}\rho_{B}e^{iH_{i}t/\hbar}\simeq\rho_{i}$, where $\rho_{i}=e^{-\beta H_{i}}/\mbox{Tr}_{f}[e^{-\beta H_{i}}]$ for $i=d$ or $a$. With this assumption, we approximate Eqs. (40) and (41) as follows:

where we define the brackets as $\langle\dots\rangle_{i}=\mbox{Tr}_{f}[\dots\rho_{i}]$. According to the definition of the Hamiltonian in Eq. (35), the correlation function in the integrand of Eqs. (42) and (43) is expressed as

where

for $i=d$ or $a$. Here, $J_{i}(\omega)$ is the density of states with the exciton–phonon interaction

and it has the property $J_{i}(-\omega)=-J_{i}(\omega)$[26, 27]. Here, $n(\omega)$ is the Bose–Einstein distribution function, and $\lambda_{i}$ is the reorganization energy of the $i$th state defined as $\hbar\int_{0}^{\infty}d\omega J_{i}(\omega)/\omega$.

The factors $k_{ad}$ and $k_{da}$ are equivalent to the rate constant expression in Förster theory[17]. The factor $K_{da}$ is a constant that is involved in nonlinear effects and the integrand has the form of the difference between the time correlation functions of forward and backward reactions. Owing to the expressions in Eqs. (38) and (39), if the off-diagonal elements $r_{ad}(t)$ and $r_{da}(t)$ of the density operator corresponding to coherence are zero as in the initial condition, the effect of the nonlinear term does not appear. In addition, even when the values of the diagonal elements of the density operators $r_{dd}(t)$ and $r_{aa}(t)$ are nearly equivalent, the effect of the nonlinear term in Eq. (39) is small. Taking this aspect into consideration, the next section describes numerical calculations of the nonlinear QME, insofar as the initial condition has some electronic coherence in the relevant system.

For the numerical analysis, we assumed that $J_{i}(\omega)$ is independent of the state $|i\rangle$ by denoting it as $J(\omega)$, and we used the overdamped Brownian oscillator model given by

Furthermore, we used the following parameters: exciton coupling $V_{ad}=20$ cm^-1, reorganization energy $\lambda=200$ cm^-1 (which is 10 times larger than the value of $V_{ad}$) to satisfy the perturbative approximation), cutoff energy $\gamma=500$ cm^-1 (which is proportional to the inverse of the relaxation time $\tau_{B}$ in the nuclear fluctuations at about 10 fs) to satisfy the Markov approximation, and temperature $T=300$ K. The energy gap $E_{d}-E_{a}$ between the donor and acceptor was 500 cm^-1. It was apparent that the value of $V_{ad}$ was smaller than the other parameters. In the site representation, the initial conditions for the numerical analysis were $r_{dd}(0)=0.9$, $r_{aa}(0)=0.1$, and $r_{da}(0)=r_{ad}(0)=\sqrt{r_{dd}(0)r_{aa}(0)}$.

We numerically calculated the time-dependency of the diagonal element of the density matrix by nonlinear QME (NLQME) in the case where $C=1$ up to $t_{max}=5$ ps, which was the order of the upper limit of the time region estimated as $\hbar\gamma/V_{ad}^{2}$ by the second-order perturbation. Fig. 2 shows the results by NLQME, linear QME (LQME) for the case of $C=0$, and hierarchical equation (HEOM)[6]. Indeed there were some differences between the results by NLQME (LQME) and HEOM in shorter intervals of time ($t<1$ ps) due to the approximation by the Markovian limit. However, the results by NLQME were similar to those by HEOM when $t>1$ ps.

To numerically analyze the energy gap $E_{d}-E_{a}$, the coupling strength $V_{ad}$, and the reorganization energy $\lambda$ dependencies of the diagonal density matrix element $r_{dd}(t)$ by NLQME, LQME, and HEOM, we define the average of difference $\Delta_{{\rm L-HEOM}}$ ($\Delta_{{\rm nonL-HEOM}}$) between the results by NLQME (LQME) and those by HEOM as:

where $r^{{\rm nonL}}_{dd}(s)$ is the diagonal density matrix element calculated by NLQME and $r^{\rm L}_{dd}(s)$ is that by LQME, and $r^{\rm HEOM}_{dd}(s)$ is the diagonal element of the reduced density matrix by HEOM.

The results of the coupling strength $V_{ad}$ and the energy gap $E_{d}-E_{a}$ dependencies of the average differences $\Delta_{{\rm nonL-HEOM}}$ and $\Delta_{{\rm L-HEOM}}$ are shown in Fig. 3. In the region where the energy gap is considerable, the NLQME results match the HEOM results relatively well.

In the region of small $V_{ad}$, which corresponds to the applicable region from Förster theory, the results by NLQME are more accurate than the results by LQME. Hence, the nonlinear term in NLQME plays an important role in such a case. On the other hand, as the value of $V_{ad}$ increases, both results by LQME and NLQME fail to match the HEOM results. However, the value of $\Delta_{{\rm nonL-HEOM}}$ is smaller than that of $\Delta_{{\rm L-HEOM}}$. This implies that the applicable parameter region for NLQME is more widened than that for LQME. In other words, since the applicability by NLQME improves in the direction of increasing $E_{d}-E_{a}$ and $V_{ad}$, NLQME can extend the applicability of the theory from the limit of Förster theory to the intermediate coupling region.

Fig. 4 shows the reorganization energy $\lambda$ and the energy gap $E_{d}-E_{a}$ dependencies of the average differences. The reproducibility of NLQME with respect to the reorganization energy is better than that of LQME when the energy gap is more than approximately 400 cm^-1. However, there is a region where the results by LQME are better than those by NLQME—namely, when the reorganization energy $\lambda$ is over 400 cm^-1 and the energy gap $E_{d}-E_{a}$ is smaller than $\lambda$.

At present, many theoretical studies employing QME are based on conventional projection operators, for which the dynamics of the irrelevant system must be in strong thermal equilibrium regardless of the state of the relevant system. Here, we focused on the dynamics of the relevant system interacting with the irrelevant system regardless of the equilibrium state, because the time-dependent projection operator techniques provide nonequilibrium dynamics in the irrelevant system.

We applied time-dependent renormalization to the conventional nonlinear QME. Although Linden and May also derived a time convolution QME using the time-dependent projection operator by Willis and Picard, they found that the nonlinear terms in the derived equation cancel identically in the Markov approximation[23]. Hence, their equation coincides with the conventional LQME. This cancellation can be attributed to simple perturbative expansion, corresponding to the case of $C=0$ in Eqs. (31) or (32) in our theory. By contrast, our proposed formula has a new nonlinear term and the derived nonlinear terms do not cancel in the Markov approximation. These differences could be due to the renormalization in the case of $C\neq 0$ in Eqs. (31) or (32) in our theory. Recently, the formalism for a closed equation of the reduced density operator based on the self-consistent Mori projector technique (c-MoP) was offered by Degenfeld-Schonburg and Hartmann[24]. Their numerical analysis in the stationary state demonstrated highly accurate results with the c-MoP comparable with mean field theory. The c-MoP derives closed nonlinear quantum master equations for the N-body system using traditional projection operators. Derivation of the equation involves focusing on each element in the reduced density matrix. Specifically, the derived equations relating to the relevant system are represented exclusively by each element of the reduced density matrix. By assuming translation invariance by the symmetry of lattice, the derived QME was reduced to equations for cluster of subsystems. In the following, to compare our method by the time-dependent renormalization to that by the c-MoP or by Willis and Picard, we briefly summarize the approximations in this paper.

The starting point is the time-convolutionless equations (11) and (12). We applied time-dependent renormalization $H^{0r}(t)$ and $H^{\prime r}(t)$ as Eqs. (13) and (14) to these equations in section 3.1. In deriving Eq. (28) from the renormalized time convolutionless Eq. (25), we performed second-order perturbative expansion with renormalized perturbative term $H^{\prime r}(t)$, which corresponds to Born approximation. The first term on the r.h.s of Eq. (30) which approximates the reduced density matrix of the irrelevant system $R_{\mathcal{I}}(t)$ with first order perturbation, corresponds to that of the initial condition $R_{\mathcal{I}}(0)$ by Born approximation, and the second term indicates that the dynamics of the relevant system is considered. Since the effect of the second term on the r.h.s of Eq. (30) is assumed to be perturbatively weak, it causes small fluctuations around the thermal equilibrium state. One of the reduced density matrices of the relevant system $r_{\mathcal{I}}(s)$ in the second and third terms of the right-hand side of Eq. (31) depends on the time variable $s$ indicating non-Markovian character. These originate because $r_{\mathcal{I}}(s)$ of the second term is due to the time-dependent projection operator $Q(s)$ in Eq. (28), and $r_{\mathcal{I}}(s)$ of the third term occurs under the influence of the dynamics in the irrelevant system of Eq. (30). By using Markov approximation, $r_{\mathcal{I}}(s)$ can be rewritten as $r_{\mathcal{I}}(t)$ and then the equation is approximately expressed as Eq. (32). In the Markov limit, the final formula is expressed as Eq. (33), where the effect of renormalization is included on the third term on the r.h.s of this equation. The first and second terms of the r.h.s of Eq. (33) are the same as the conventional linear QME. The time correlation function in the third term of the r.h.s of Eq. (33) is expressed similar to that in the second term (see Appendix). Therefore, timescales $\tau_{B}$ and $\tau_{R}$ in the Born-Markov approximation ($\tau_{B}\ll\tau_{R}$) applied in Eq. (33) are the correlation time of the bath ($\tau_{B}\sim\hbar/\gamma$ in this study) and the relaxation time of the relevant system ($\tau_{R}\sim 1/k_{ad}$ in this study) when approximately derived by the conventional linear QME. The final form of our QME to which the secular approximation ($\tau_{S}\sim\hbar/|E_{d}-E_{a}|\ll\tau_{R}$ in this study) is applied for numerical analysis is expressed by Eqs. (38) and (39). Our QME is expressed as closed formula only by the reduced density matrix of the relevant system. The nonlinear terms that include the coefficient $C$ in our QME are new terms yielded by the time-dependent renormalization procedure in this study, while the new coefficient $C$ does not appear in the methods by Willis and Picard, and by the c-MoP. In the c-MoP method, Born approximation is applied in the derivation of the QME. In the method by Willis and Picard, the Born approximation, the Born-Markov approximation and the approximation of the Markov limit are stepwise applied in the derivation of each parallel QME for the relevant and irrelevant systems. In this study, at each stage of approximations summarized above, approximation is performed so that the contribution of the coefficient $C$ included in the renormalized perturbation term remains in the lowest order. Depending on the value of the coefficient $C$, the non-linear terms in our QME contribute to correct the conventional linear QME in the relevant system through the dynamics in the irrelevant system. This is not included in the methods by the c-MoP, and by Willis and Picard, and is a new effect added by this study.

In this paper, through the renormalization technique in which the dynamics of the irrelevant system are renormalized into the dynamics of the relevant system, the renormalized irrelevant system appears to achieve thermal equilibrium in a mathematical sense. In addition, the strength of the renormalized perturbative terms $H^{\prime r}$ should be less than that of the original perturbation terms $H^{\prime}$. Hence, the applicable timescale can be extended beyond that without the renormalization technique. The conspicuous difference between the numerical analysis of our proposed formalism and that of the conventional analysis appears when the timescale is much larger.

The inset of Fig. 2 shows the numerical results over a short period of time. It can be seen that the initial condition for numerical analysis with our theory is same as that by HEOM. With HEOM, however, the probability immediately following $t=0$ increases slightly, because the process is non-Markovian over a short period of time. In this study, we assumed that the value of $\gamma$ is large. This means that the relaxation time is very fast. When the value of $\gamma$ is small, the reaction process is non-Markovian, and the upper limit of the time integration $t_{max}$ in Eqs. (48) and (49) is small. In such a case, since the values of Eqs. (48) and (49) are large, quantitative discrepancies appear between the numerical results by our theory and those by HEOM, as shown in Figs. 3 and 4.

The numerical results in Figs. 3 and 4 show that the proposed theory is strongly dependent on the size of the energy gap. When the size of the energy gap is smaller than the thermal energy, with $|E_{d}-E_{a}|<k_{B}T$, the state immediately after transition at site $|d\rangle$ $\to$ site $|a\rangle$ in the relevant system is near the thermal equilibrium state for the irrelevant system without renormalization. In such a case, the dynamics for the irrelevant system may be applicable as a thermal bath. However, in the opposite case, where $|E_{d}-E_{a}|>k_{B}T$, the state immediately after transition in the relevant system passes far from the thermal equilibrium state for the irrelevant system without renormalization. In such situations, insofar as there is considerable effective relaxation time in the environment before attaining the steady state, the dynamics of the environment might be far from the equilibrium state immediately after the occurrence of a state transition. Hence, the perturbative approach using the conventional projection operator is inapplicable in such a situation. The effect of renormalization in this study should account for the feedback of such nonequilibrium states on the dynamics of the irrelevant system.

However, there were some quantitative discrepancies between the results of the proposed theory and those of the HEOM. This may be due to the delocalization effect caused by the small energy gap, because the electronic state approaches resonance in such cases. In this paper, our numerical analysis considered the transfer process in terms of the theory for the weak electronic coupling case exclusively. We may need to analyze the dynamics in the opposite limiting case as well, that is, using electron-phonon coupling as the perturbation term. The resulting QME could be applied to the relaxation processes in the delocalized exciton states. To develop our theory more generally, the discrepancies between it and HEOM could be corrected using the perturbative approach in the modified Redfield theory [28, 26], which reproduces Förster theory and Redfield theory for the limiting cases. Recently, a coherent modified Redfield theory [29, 30] was constructed using a pure dephasing reference system master-equation method [31] based on the modified Redfield theory. However, to use this modified Redfield theory, we may need to derive a more exact nonlinear QME in the intermediate coupling case.

After the analysis shown in Fig. 4, it is clear that the theory by nonlinear QME, namely for the case of $C=1$ in Eq. (32), is difficult to apply in all parameter areas, where it holds that the energy gap $E_{d}-E_{a}$ is smaller than the reorganization energy $\lambda$. In this parameter region, it is evident that the linear master equation for the case of $C=0$ in Eq. (32) achieves better results than the nonlinear formalism. This shows that there may be a relationship between the parameter region and the variable $C$. In other words, the degree of the feedback effect from the irrelevant system changes depending on the parameter region. In order to include this effect, it is necessary to introduce a formalism that adjusts the variable $C$ according to the parameter region.

We previously researched the application of a variational principle to QME (which results in a so-called variational master equation) [11, 12]. According to this approach, the free energy function is used to obtain an optimized trial function for the variational method. The variational parameters in the variational master equation should generally vary as functions of time, although the optimized variational parameter is applicable to the thermal equilibrium state. Although some trial functions for the polaron problem have been proposed [13, 14, 15, 16], from the point of view of the polaron transformation, we suppose that the time-dependent renormalization in this paper corresponds to a polaron transformation. Specifically, our numerical analysis was limited exclusively to strong constraints in the case of $C=1$ for the NLQME. If weaker conditions could be imposed by including the variable $C$ as a variational parameter, these problems might be partially or fully overcome.

The numerical results for the proposed theory show that this theoretical approach can be applied over large timescales. For example, in a photosynthetic antenna and its reaction center, a fast light-harvesting process occurs, followed by a slower electron transfer reaction. In the leading studies in this field, attempts have been made to understand the mechanism that produces electron transfer reactions after excitation energy is transferred to the reaction center [32, 33]. Each of these elementary reactions has an intrinsic timescale, and they are hierarchically different from each other. The proposed theory could apply to cases in which the thermalization process produced by the nonequilibrium dynamics of the environment has an important influence on the reactions.

Finally, the only modification to the approximated final form in the proposed theory is the inclusion of the nonlinear term in Eq. (31). Concerning the numerical analysis, such a modification is easy to implement in the source program, and the computational cost is almost the same as that for conventional QME. Hence, the proposed theory offers advantages for the analysis of reaction dynamics over large timescales.

In this paper, a new QME was constructed by applying time-dependent renormalization to the nonlinear QME derived by Shibata and Hashitsume[25]. The general expression of the proposed equation, which might account for the influence of dynamic environmental feedback on the irrelevant system, has nonlinear terms in addition to the conventional linear QME. As an example, the transfer dynamics in relevant two-site systems with a nuclear vibrational environment based on the conventional theory for the weak electronic coupling case were numerically analyzed in a comprehensive parameter region. We found that, when the energy gap is large, the proposed theory reproduces the results given by HEOM, even when the differences between the conventional theory and HEOM are considerable. Consequently, the proposed theory offers significant advantages when studying the dynamics of a system over large timescales, particularly when the thermalization caused by the non-equilibrium dynamics of the environment has an important influence on the reactions.