Abstract

Even photosynthesis—the most basic natural phenomenon underlying life on Earth—involves the nontrivial processing of excitations at the pico- and femtosecond scales during light-harvesting. The desire to understand such natural phenomena, as well as interpret the output from ultrafast experimental probes, creates an urgent need for accurate quantitative theories of open quantum systems. However it is unclear how best to generalize the well-established assumptions of an isolated system, particularly under nonequilibrium conditions. Here we compare two popular approaches: a description in terms of a direct product of the states of each individual system (i.e., a local approach) versus the use of new states resulting from diagonalizing the whole Hamiltonian (i.e., a global approach). The main difference lies in finding suitable operators to derive the Lindbladian and hence the master equation. We show that their equivalence fails when the system is open, in particular under the experimentally ubiquitous condition of a temperature gradient. By solving for the steady state populations and calculating the heat flux as a test observable, we uncover stark differences between the formulations. This divergence highlights the need to establish rigorous ranges of applicability for such methods in modeling nanoscale transfer phenomena—including during the light-harvesting process in photosynthesis.

1. Introduction

The rapid recent development of quantum control, quantum engineering, and information processing at the nanoscale in biological, chemical, and condensed matter systems has led to a crucial need to improve our understanding of open quantum systems. The typical physics assumptions of an isolated system, particularly under nonequilibrium conditions, cannot hold for systems probed on the quantum scale at optical pico- and femtosecond scales [1]. The accurate description of excitation dynamics at these ultrafast scales is essential to understand fundamental processes for Life on Earth such as photosynthesis [24]. For such reasons, theoretical physicists have begun to develop theoretical and experimental tools to study the dynamical behavior of an open system interacting with its environment. The key lies in identifying an accurate way of removing the environmental degrees of freedom and hence obtain a closed equation of motion for the reduced system of interest [59]. Quantum master equations have been tested for exactly solvable interacting harmonic systems in thermal equilibrium and nonequilibrium conditions providing reliable results [10]. However, many systems of interest (e.g., optically probed semiconductor quantum dots or biological photosynthetic processes) show fluctuations which are far from equilibrium. Moreover, recent evidence suggests that excitation energy transfer in biological systems, particularly photosynthetic membranes, might involve some level of quantum coherence [11, 12].

In this paper, we provide analytic results for the steady state of an open quantum system interacting with two reservoirs at different temperatures following local and global approximations. This mimics the situation in photosynthetic membrane reaction centers, and elsewhere in natural and artificial systems, in which transfer occurs between few-level molecular complexes that in turn may be coupled to reservoirs with different effective temperatures. In both cases we use the same assumptions and the same procedure; the only difference is the election of the free Hamiltonian and, therefore, the basis set. The local approach is commonly used to model incoherent transfer phenomena in small quantum systems [4], while the global approach is a more rigorous way to calculate the quantum properties such as coherence and entanglement [13, 14]. For the simple case of an interacting qubit dimer where one of the qubits is coupled with a thermal bath, the approximations are equivalent only for the case of zero bath temperature [15]. Similar studies have been performed for two interacting quantum harmonic oscillators under nonequilibrium thermal conditions, where the local approach is found to violate the second law of thermodynamics for the nonsymmetrical case [16] (nonidentical systems). This discrepancy represents a serious challenge for modeling of quantum systems. We consider an open system comprising interacting subunits (qubits), which could provide a simple representation of interacting two-level systems in a photosynthetic membrane reaction center. The question of how to solve the problem can then be addressed in two ways. (i) Diagonalize the Hamiltonian of the open system and solve the problem in terms of the diagonal basis set (i.e., global approach). (ii) Use the direct product of the individual basis set of the interacting subunits (i.e., local approach). We here apply both these methods in parallel and show explicitly how each formulation leads to a different result when the number of interacting subunits is greater than one.

This paper is divided in four parts. In the second part the methods are presented for an arbitrary system. Analytic expressions for populations and heat current are derived and applied to a two-level system. In the third part, the quantum system is extended to an interacting qubit dimer where the quantities are calculated. The fourth part is devoted to analysis and conclusions.

2. Formalism

Consider a quantum system under nonequilibrium thermal conditions. Each reservoir is modeled as an infinite collection of harmonic oscillators in thermal equilibrium at temperature given by ,  . We assume that the coupling strength between the reservoirs and the central subsystem is weak; hence we can express the total density operator as a direct product of the reduced density operators of the open system and the reservoirs and : . The Hamiltonian of the whole system is then where is the Hamiltonian of the open system and the terms and are the reservoir Hamiltonian and the interaction term associated with the reservoir , respectively. The usual route to this problem is to solve directly the second-order expansion of the Liouville von Neumann equation. This is achieved by obtaining the reduced density operator for the central subsystem through a partial trace over the reservoirs degrees of freedom. Within the interaction picture representation, the dynamics for the whole system is given by The Hamiltonian in the interaction picture representation is defined as , where is the evolution unitary operator of the free system. The free system is considered as a subsystem whose solutions are well known. For the present work, we use the two methods to describe the free system: global and local.

2.1. Global Approach

For this global approach, it is useful to express the interaction term in the form where the operators act over the open system’s degrees of freedom while the operators act on the reservoir . The operators are chosen in such a way that they follow the following commutation relationship: This decomposition is always possible [5, 10, 17]. For instance, for an operator , where and are eigenstates of , it results in . The dynamics, as was mentioned above, is governed by the Liouville von Neumann equation of motion . By using the Born-Markov approximation, the dynamics of the open system in terms of its reduced density operator is given by where is the Lindbland superoperator associated with reservoir [13] where the indices and run over the complete range of operators and is the spectral density is defined by [10, 13, 17] where is the interaction picture representation of reservoir operator . For each selection of the open system, a new set of operators is found.

2.2. Local Approach

Similarly, the problem can be addressed by using a free Hamiltonian which is formed by summing all zero-point Hamiltonians of each subunit. For a simple subunit such as the qubit, the Hilbert space is spanned by two states with energy gap of , and the Hamiltonian can be written in terms of the Pauli matrices as . For instance, for the case where the open system is composed by a linear chain of qubits and considering an XX-like interaction between them, the Hamiltonian of the open system can be written as , where is the contribution of each subsystem and . describes the interqubit interaction. Consequently, the free Hamiltonian and interaction Hamiltonian can be identified as Hence in the interaction picture representation, the Hamiltonian can be written as where creates an excitation in mode of reservoir with a coupling strength of and energy . The subindex labels the subunit interacting with the reservoir ; that is, when and for . Note that for the energy associated with the central system we use only one subindex, for example, , while for the energy associated with one mode of the bath we use two, for example, . We use the commutation relations for Pauli operators and and for bosonic operators . As a result of this transformation, we can use (2) to solve the problem. Specifically, we take the partial trace over the reservoirs and use the Born-Markov approximation. Furthermore, we take the continuous limit for the reservoir frequencies and the wide band limit on the interaction with the central subsystem. This procedure yields to a delta function in the energies that collapses all the energy spectra of reservoir into the energy of the subunit . Under these conditions, the quantum optical master equation for a chain of qubits, whose endpoints ( and ) interact with bosonic reservoirs, is given by where denotes the spectral density function associated with the interaction between the qubit and the reservoir given in terms of the spontaneous emission rate and the Bose-Einstein distribution For simplicity, we consider .

2.3. Observable: Thermal Energy

As a test observable, we use the steady state thermal energy transferred from the reservoir to the quantum system given the nonzero thermal bias. The steady state heat flux is defined as the trace of the product of the Liouvillian superoperator with the subsystem Hamiltonian: The resultant expression for the heat flux in the global approach can be written in the following compact form: In a similar way, the result for the heat flux in the local approach is found to be As expected, formula (14) gives a zero flux value when the reservoirs are set at the same temperature. Furthermore, it can be seen that the global and local approaches lead to the same result when applied to a single qubit under nonequilibrium thermal conditions where the system operators are the same, . The steady state population of the excited state of one qubit system with energy gap is found to be where is a simple universal function defined as Furthermore, the heat flux for this system is found to be As can be seen, formula (18) leads to a zero value when the thermal bias is set as zero. In addition, note the heat flux is always positive for positive bias () and negative, otherwise.

3. Divergence for a Dimer

We consider a dimer composed of two interacting qubits, where each quit is connected to a different reservoir in thermal equilibrium at temperature   . Figure 1 shows schematically the system to be studied. The Hamiltonian of the qubit dimer is where and , with as the identity matrix. In addition, we assume that the qubit labeled as interacts with the bath labeled as only, for . Hence the interaction Hamiltonian between the dimer and the reservoirs can be written as The operators do not commute with , so for the global approach it is necessary to find a transformation that allows us to write the interaction Hamiltonian in the form of (3), so that condition (4) is fulfilled. The eigenstates and eigenenergies of are   ,   , (), and (), where the amplitudes and constants are given by The transformation of the coupling operators from the individual qubits basis set into the dimer diagonal basis set can be calculated as With this transformation, condition (4) is fulfilled . In this way, the operators can be found to be

For simplicity we consider the symmetric case where the energy gap of each qubit is the same: . In particular, we have that , , and the amplitudes are and . The steady state populations can be calculated for each of the qubits: where the matrix elements are calculated in the diagonal basis. For the weak-coupling case (), the populations are found to be () where . For the steady state, the matrix element is zero; therefore the populations take the simple form of . In a similar way the populations for the local approach can be derived by solving (10) for leading to , where the matrix elements are calculated in the local basis; that is, , , and . The following results follow: In Figure 2 we can see how the populations change as the temperature of reservoir changes while the temperature of reservoir is set to be zero. The two approaches for converge when . On the other hand, population decays in the local approach as while the global approach predicts an asymptotic growth to the mixed state of . This divergence in predictions suggests that one of the approaches is not correct.

Another way to see this discrepancy is by looking at the heat flux. Using the system operators (23), the steady state heat flux (13) for reservoir is The heat flux is expressed as a sum over the energy channels and that depend on the interqubit coupling and qubit energy gap . The flux is always positive for positive thermal bias. On the other hand, the steady state heat flux from reservoir to qubit , calculated in the local approach for the symmetric case, can be found to be The local approach for the dimer leads to an expression for the heat flux expression equal to the one for the monomer weighted by a positive function that depends on the interqubit coupling and the reservoirs' temperature. This function ensures that the flux tends to zero as the interqubit coupling decreases and remains finite for large . This is reasonable since the qubits are weakly coupled (i.e., ), and therefore the heat transferred should decrease. However it also shapes the flux in such a way that an optimal value is found; that is, the heat flux exhibits a maximum with the temperature and then decays to zero as the thermal bias grows. This behavior is not found for either the monomer or the dimer in the global approach.

As an illustration, Figure 3 shows the observable for both approximations and three different values of the qubit energy gap , as a function of the temperature difference. To simplify the presentation, we have set the temperature of the second reservoir to be close to zero and have depicted the flux as a function of . We can clearly see the maximum of for a specific value of the energy in the local approach, while the result for the global approach grows asymptotically to as the thermal bias increases.

The dimer system thus provides the simplest nontrivial physical scenario where a quantitative comparison of quantum nonequilibrium thermal quantities can be performed: those obtained within a rigorous global approach and results coming from a less rigorous and hence restricted in validity, local approach. Furthermore, our results as tested in simple systems show that some physical magnitudes calculated within the local framework could be misleading. This work suggests various future avenues of research, one of which concerns the systematic analysis of the scaling of the divergence “distance” between different approach results as the system size (complexity) and the separation from equilibrium increase.

4. Conclusions

We have shown that even in the symmetrical case, the problem of an open quantum system interacting with two reservoirs at different temperatures leads, within local and global approximations, to different results when the number of interacting subunits in the open system is greater than one. The formulations are equivalent and identical for the monomer. By contrast, the results for an interacting dimer open several urgent questions about the range of applicability of the underlying approximations and hence methods. In the global approach, the populations for both qubits are identical in the steady state. On the contrary, the local approach predicts that the population of the qubit interacting with the cold reservoir is always smaller than the population of the other qubit. Second, the local approach predicts a maximum in the heat flux as a function of the temperature gradient, followed by a gradual decay to zero as the thermal bias grows. By contrast, the global approach predicts a saturation of the flux as the bias increases. Finally, we note that the outcome of the local approach in the strong interqubit coupling limit concludes that the dimer can be modeled as a single qubit which resembles the properties of a classical system. Future work with larger numbers of qubits will elucidate the differences in these approaches, as will a careful comparison to future experiments which are able to distinguish between the divergences in their predictions.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.