Non-equilibrium quantum systems: Divergence between global and local descriptions

Even photosynthesis -- the most basic natural phenomenon underlying Life on Earth -- involves the non-trivial 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 non-equilibrium 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). 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.


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 non equilibrium conditions, cannot hold for systems probed on the quantum scale at optical pico-and femtosecond scales [1]. Given that the most fundamental process for Life on Earthphotosynthesis -also processes excitations on these ultrafast scales, has raised the question of whether we truly understand its dynamics [2,3,4]. 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 [5,6,7,8,9].
However, many systems of interest (e.g. optically probed semiconductor quantum dots, or biological photosynthetic processes) show fluctuations which are far from equilibrium. For example, recent evidence suggests that excitation energy transfer in biological systems, particularly photosynthetic membranes, might involve some level of quantum coherence [10,11].
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 [12,13]. 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 [14]. Similar studies have been performed for two interacting quantum harmonic oscillators under non equilibrium thermal conditions, where the local approach is found to violate the second law of thermodynamics for the non-symmetrical case [15] (non identical 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.

Formalism
Consider a quantum system under non-equilibrium thermal conditions. Each reservoir is modeled as an infinite collection of harmonic oscillators in thermal equilibrium at temperature given by β j = 1/k B T j , j = 1, 2. We assume that the coupling strength between the reservoirs and the central sub system is weak, hence we can express the total density operator as a direct product of the reduced density operators of the open systemρ s and the reservoirsρ 1 andρ 2 :ρ =ρ s ⊗ρ 1 ⊗ρ 2 . The Hamiltonian of the whole system is thenĤ whereQ is the Hamiltonian of the open system and the termsR j andŜ j are the reservoir Hamiltonian and the interaction term associated with the reservoir j, 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Ĥ I (t) =Û † 0Ĥ IÛ0 , whereÛ 0 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.

Global approach
For this global approach, it is useful to express the interaction term in the form where the operatorsV j,µ act over the open system's degrees of freedom while the operatorsf j,µ act on the reservoir j. The operatorsV j,µ are chosen in such a way that they follow the following commutation relationship: The dynamics, as was mentioned above, is governed by the Liouville von Neumann equation of motion dρ/dt = −i[Ĥ,ρ]. By using the Born-Markov approximation, the dynamics of the open system in terms of its reduced density operatorρ s is given by: where L j is the Lindbland super-operator associated to the reservoir j [12] where the indices µ and ν run over the complete range of operators and J (j) µ,ν (ω ν ) is the spectral density which is given by wheref j,ν (τ ) = e −iRj τf j,ν e iRj τ is the interaction picture representation of the reservoir operatorf j,ν . For each selection of the open system, a new set of operators {V j,µ } is found.

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 2 × 2 pauli matrices asQ = (1/2)σ z . For instance, for the case where the open system is composed by a linear chain of N qubits and considering an XX-like interaction between them, the Hamiltonian of the open system can be written asQ =Q 0 +Q I , whereQ 0 = N q=1 describes the inter-qubit interaction. Consequently, the free HamiltonianĤ 0 and interaction HamiltonianĤ I can be identified aŝ Hence in the interaction picture representation, the Hamiltonian can be written aŝ whereâ † k,j creates an excitation in mode k of reservoir j with a coupling strength of g (j) k , and the subindex λ labels the subunit interacting with the reservoir j, i.e. λ = 1 when j = 1 and λ = N for j = 2. We use the commutation relations for Pauli operators ([σ + , σ − ] = σ z and [σ z , σ ± ] = ±2σ ± ), and for bosonic operators ([â k ,â † k ] = δ k,k ). As a result of this transformation, we can use equation 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 allow us to collapse, through the delta function δ( λ − j,k ), all the frequency spectra of the reservoir j into the frequency of the subunit λ. Under these conditions, the quantum optical Master equation for a chain of qubits whose endpoints (λ = 1 and λ = N ) interact with bosonic reservoirs, is given by where J (j) denotes the spectral density function associated with the interaction between the qubit λ and the reservoir j given in terms of the spontaneous emission rate γ j and the Bose-Einstein distribution For simplicity, we consider γ j = 1.

Observable: Thermal Energy
As a test observable, we use the steady state thermal energy transferred from the reservoir to the quantum system given the non-zero thermal bias. The steady-state heat flux is defined as the trace of the product of the Liouvillian super-operator 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, the formula (17) 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 non-equilibrium thermal conditions where the system operators are the same, V j,µ =σ +(−) . The steady-state population of the excited state n of one qubit system with energy gap is found to be n = σ +σ− = 1 2 e( ) where e j is a simple universal function defined as Furthermore, the heat flux for this system is found to be . (20) As can be seen, the formula 21 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 (N 1 >N 2 ) and negative, otherwise.

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 T j (j = 1, 2). Figure 1 shows schematically the system to be studied. The Hamiltonian of the qubit dimer iŝ whereσ 1 =σ ⊗I andσ 2 = I ⊗σ, with I as the 2×2 identity matrix. In addition, we assume that the qubit labeled as j interacts with the bath labeled as j only, for j = 1, 2. Hence the interaction Hamiltonian between the dimer and the reservoirs can be written aŝ The operatorsσ ± j do not commute withQ, so for the global approach it is necessary to find a transformation that allow us to write the interaction Hamiltonian in the form of Equation (3), so that the condition (4) is fulfilled. The eigenstates and eigenenergies ofQ are 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σ |s p s p |σ j |s q s q | (26) With this transformation, the condition (4) is fulfilled Q , |s p s q | = (E p − E q ) |s p s q |. In this way, the operators can be found to bê For simplicity we consider the symmetric case where the energy gap of each qubit is the same: 1 = 2 = . In particular, we have that ω 1 = ω 4 = | − K|, ω 2 = ω 3 = + K and the amplitudes are the qubits where the matrix elements ρ ij = s i |ρ s |s j are calculated in the diagonal basis. For the weak-coupling case ( > K), the populations are found to be (K = 1): where e j = e(ω j ). For the steady-state, the matrix element ρ 34 is zero, therefore the populations take the simple form of n 1,2 = (e 1 + e 2 )/4. In a similar way the populations for the local approach can be derived by solving equation (13) for N = 2 leading to n 1(2) = ρ 11 + ρ 22(33) , where the matrix elements ρ ij are calculated in the local basis, i.e. ρ 11 = 11|ρ s |11 , ρ 22 = 10|ρ s |10 , and ρ 33 = 01|ρ s |01 . The following results follow: In Figure 2 we can see how the populations change as the temperature of the reservoir 1 changes while the temperature of reservoir 2 is set to be zero. The two approaches for n 1 converge when k B T 1 >> . On the other hand, the population n 2 decays in the local approach as k B T 1 >> while the global approach predicts an asymptotic growth to the mixed state of 1/2. 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 (27), the steady state heat flux (16) for reservoir 1 is: The heat flux is expressed as a sum over the energy channels ω 1 and ω 2 that depend on the inter-qubit coupling K and qubit energy gap . The flux is always positive for positive thermal bias. On the other hand, the steady state heat flux from the reservoir 1 to the qubit 1, 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 K and the reservoirs' temperature. This function ensures that the flux tends to zero as the inter-qubit coupling decreases and remains finite for large K. This is reasonable since the qubits are weakly coupled (i.e. K → 0), and therefore the heat transferred should decrease. However it also shapes the flux in such a way that an optimal value is found, i.e., 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 Q 1 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 k B T 1 . We can clearly see the maximum of Q 1 for a specific value of the energy k B T 1 in the local approach, while the result for the global approach grows asymptotically to as the thermal bias increases.

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 qubit 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 inter-qubit 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 divergence in their predictions.