Nonequilibrium steady-state picture of incoherent light-induced excitation harvesting

The observation of unexpectedly long-lived oscillatory features of ultrafast spectroscopic signals measured on photosynthetic pigment–protein complexes Engel et al. (2007); Panitchayangkoon et al. (2010) has generated much excitement during the last decade. The idea that (dynamic) coherences modulating these signals may be directly relevant to natural light harvesting, Cheng and Fleming (2009); Ishizaki and Fleming (2011) which is triggered by stationary incoherent sunlight, Brumer and Shapiro (2012); Brumer (2018) has been driving intense research activities in the field. Mančal (2020); Cao et al. (2020) Insights from ultrafast spectroscopies are indispensable in determining the underlying Hamiltonian of the system under investigation, from which the dynamics of excitation energy transfer (EET) under any excitation condition may be inferred (as long as the excitation is sufficiently weak). Mančal and Valkunas (2010); Chenu, Malý, and Mančal (2014); Brumer (2018); Scholes (2018); Mančal (2020) A number of experimental Engel et al. (2007); Panitchayangkoon et al. (2010) and theoretical Ishizaki and Fleming (2009a); Sarovar et al. (2010) studies speculating about a positive impact of coherences and entanglement on the biological EET have been performed on the so-called unloaded systems. Such systems feature no transmission of photoinduced excitations to the reaction center (RC or load) from which the excitation energy is eventually harvested. Furthermore, the time scales addressed in these studies are generally much shorter than those representative of exciton recombination, either radiative or nonradiative. van Amerongen, Valkunas, and van Grondelle (2000) Direct experimental insights into the dynamics of molecular aggregates initiated by incoherent light are limited. Turner et al. (2013) Therefore, at present, the possible relevance of some sort of quantum coherence for EET under natural conditions is best examined within appropriate theoretical models.

Such models should contain a realistic description of photoexcitation by natural incoherent light, whose intensity is essentially constant from the molecular viewpoint. Thus, the physically plausible description of natural light harvesting should feature continuous generation of electronic excitations by light, their continuous delivery to the RC, and their continuous loss by recombination. Kassal, Yuen-Zhou, and Rahimi-Keshari (2013); Tscherbul and Brumer (2018); Manzano (2013); Brumer (2018); Jesenko and Žnidarič (2013); Chuang and Brumer (2020) The EET is then studied from the standpoint of nonequilibrium steady states (NESSs), which arise as a result of excitation photogeneration, phonon-induced relaxation, dephasing, trapping at the RC, and recombination. Furthermore, the coupling between the radiation and absorbing pigments is, in general, weak, so that its second-order treatment is reasonable. Then, the only information we need about the radiation is its first-order correlation function, Mančal and Valkunas (2010) which can be either modeled by appropriate expressions Mančal and Valkunas (2010); Sadeq and Brumer (2014) or obtained by a suitable ensemble average. Sadeq and Brumer (2014); Chenu and Brumer (2016) Excitation and deexcitation events can be treated within the Born–Markov approximation Breuer and Petruccione (2002) by Lindblad dissipators, Fassioli, Olaya-Castro, and Scholes (2012); Chan et al. (2018) or by establishing the Bloch–Redfield quantum master equation. Tscherbul and Brumer (2015, 2018) Approaching the problem from the perspective of open quantum systems, one can introduce an appropriate spectral density of light–matter coupling, Pachón and Brumer (2013); Pachón, Botero, and Brumer (2017); Olšina et al. (2014); Chuang and Brumer (2020) and possibly treat it even beyond the second order. Olšina et al. (2014)

Reasonable models of photosynthetic EET should not overlook the non-Markovian interplay between photoinduced electronic excitations and nuclear reorganization processes, Ishizaki and Fleming (2012); Chenu and Scholes (2015) whose relevance is emphasized by ultrafast spectroscopic studies. To that end, a number of studies attempt to combine an explicit treatment of the photoexcitation step with a nonperturbative approach to the excitation–environment coupling. Fassioli, Olaya-Castro, and Scholes (2012); Chan et al. (2018); Maguire, Iles-Smith, and Nazir (2019); Dijkstra and Beige (2019); Janković and Mančal (2020); Olšina et al. (2014) The method of choice for an exhaustive treatment of excitation–environment coupling are hierarchical equations of motion (HEOM). Ishizaki and Fleming (2009b); Tanimura (2020) In the accompanying paper, Janković and Mančal (2020) we combine HEOM with a second-order treatment of light–matter coupling for light of arbitrary properties. Our method correctly captures light-induced reorganization processes and nonequilibrium evolution of the bath between the two interactions with light.

Recently, a number of groups have suggested that stationary coherences in the energy basis (interexciton coherences) or local basis (intersite coherences) under incoherent illumination may improve the light-harvesting efficiency in photosynthetic systems. Tscherbul and Brumer (2018); Tomasi and Kassal (2020); Yang and Cao (2020); Jung and Brumer (2020) However, the majority of existing theoretical approaches to NESSs in photosynthetic light harvesting typically feature a simplified treatment of the photoexcitation Jesenko and Žnidarič (2013); León-Montiel, Kassal, and Torres (2014); Manzano (2013) or a simplified treatment of excitation relaxation and dephasing. León-Montiel, Kassal, and Torres (2014); Tscherbul and Brumer (2018); Manzano (2013); Yang and Cao (2020); Jung and Brumer (2020) Also, efficient algorithms that avoid the explicit temporal propagation in the computation of NESSs induced by natural incoherent light have just begun to be developed. Axelrod and Brumer (2018) Therefore, there is still an urge to construct theoretical methods that circumvent the disparity between the time scales of EET dynamics and incoherent excitation sources, and yet meet the two requirements outlined in the above text.

Another pertaining issue is the origin of stationary coherences, i.e., whether they are primarily induced by incoherent radiation or by the coupling to the protein environment. The authors of Ref. Olšina et al., 2014 argued that the NESS coherences ultimately stem from the entanglement of electronic excitations with the environment. This entanglement has been systematically studied both analytically Lee, Cao, and Gong (2012); Lee, Moix, and Cao (2012) and numerically Lee, Cao, and Gong (2012); Cerrillo and Cao (2014) within the undriven and unloaded spin–boson model. The two-level system displays noncanonical equilibrium statistics, Strümpfer and Schulten (2009); Moix, Zhao, and Cao (2012) whose deviation from the canonical equilibrium statistics can be conveniently measured by a single parameter. Lee, Cao, and Gong (2012) This parameter can be interpreted as the angle by which the basis in which the system’s Hamiltonian (or the system–bath interaction Hamiltonian) is diagonal should be rotated to obtain the diagonal reduced density matrix (RDM). The basis in which the RDM is diagonal is thus singled out by the environment, and the corresponding basis states are known as the preferred (or pointer) states within the framework of the decoherence theory. Zurek (2003); Schlosshauer (2007) The concept of preferred basis is useful whenever representation-dependent issues, such as the ones we are after in this study, arise.

In this paper, we extend the ideas developed in Ref. Lee, Cao, and Gong, 2012 to examine the properties of the NESS that arises in an incoherently driven and loaded excitonic aggregate. Adapting the algorithm presented in Ref. Zhang et al., 2017, we devise a procedure to find the NESS of our recently proposed HEOM that incorporates incoherent photoexcitation, Janković and Mančal (2020) and to properly define light-harvesting efficiency under incoherent illumination. Our theoretical approach fully respects the continuity equation for excitation fluxes. The NESS is most conveniently described in the so-called preferred basis, in which the steady-state RDM is diagonal. Such a description of a driven and loaded aggregate may be regarded as analogous to the normal-mode description of a system of coupled harmonic oscillators. While our theoretical method is quite general and applicable to arbitrary excitonic networks, we investigate the properties of the NESS using the appropriately parameterized model dimer. For realistic values of the load extraction time, we conclude that light-induced coherences are completely irrelevant in the NESS, which is then close to the noncanonical equilibrium of the undriven and unloaded dimer. We find that the NESS of the driven and loaded dimer is intimately related to the dynamics of the driven, but unloaded dimer, that takes place on the time scale of the excitation trapping at the load. This close connection between the dynamic and stationary picture is correct provided that the load extraction time is longer than the time scale of coherence dephasing, which is in principle accessible in ultrafast spectroscopies.

The paper is structured as follows: The model and method are presented in Secs. II and III, respectively. In Sec. IV, we discuss the relation of the preferred-basis description of the NESS to more standard descriptions conducted in the site or excitonic basis. Section V presents numerical results, with an emphasis on the relation between the time-dependent and stationary picture. In Sec. VI, we outline how the methodology presented in this work could be applied to multichromophoric systems. Section VII concludes the paper by summarizing its principal findings.

We consider the simplest EET system, a molecular aggregate composed of two mutually coupled chromophores (a dimer, a spin–boson-like model Zimanyi and Silbey (2012)). Although a similar model has been repeatedly used by many authors to gain insight into fundamentals of light harvesting under incoherent illumination, Jesenko and Žnidarič (2013); Olšina et al. (2014); Tscherbul and Brumer (2014, 2018); Tomasi et al. (2019); Yang and Cao (2020); Jung and Brumer (2020) our analysis uses a rigorous theoretical approach to investigate in greater detail certain properties of the NESS that have not received enough attention so far. The dimer system can be considered as the minimal model of a photosynthetic antenna with delocalization in which one can study the effects of energy relaxation, dephasing, and (static) disorder in local transition energies (the so-called asymmetric dimer, see Sec. V.1). To be specific, we speak about the model dimer, while we note that the model and the method to be presented are quite general and applicable to multichromophoric situations.

Electronic excitations of the model dimer are modeled within the Frenkel exciton model Valkunas, Abramavicius, and Mančal (2013); May and Kühn (2011) and the corresponding Hamiltonian reads as

In Eq. (1), $|l_{j}\rangle$ is the singly excited state localized on chromophore $j$, $\varepsilon_{j}$ is its vertical excitation energy, while $J_{jk}$ are resonance couplings (we take $J_{kk}\equiv 0$). We limit our discussion to the manifold of singly excited states, which is justified under the assumption that the driving by the radiation is sufficiently weak. The aggregate is in contact with the thermal bath, which represents its protein environment, and which is modeled as a collection of independent oscillators labelled by site index $j$ and mode index $\xi$

The phonon creation and annihilation operators $b_{j\xi}^{\dagger}$ and $b_{j\xi}$ entering Eq. (2) satisfy Bose commutation relations. The aggregate is driven by weak radiation and the generation of excitations is described in the dipole and rotating-wave approximations

In Eq. (3), operators $\mathbf{E}^{(\pm)}$ are the positive- and negative-frequency part of the (time-independent) operator of the (transversal) electric field, while the $eg$ part of the dipole-moment operator reads as

We assume that transition dipole moment $\mathbf{d}_{j}$ of chromophore $j$ does not depend on nuclear coordinates (Condon approximation). The interaction of photoinduced excitations with the environment is taken to be in Holstein form, i.e., it is local and linear in oscillator displacements

We assume there are two possible channels through which photogenerated excitons may decay. The first one is their transfer to the charge-separated state in the RC, in which case they are usefully harvested. On the other hand, exciton recombination, either radiative or nonradiative, is detrimental to the efficiency of EET. While our description of exciton photogeneration and the subsequent phonon-induced relaxation is exact, see Sec. III and Ref. Janković and Mančal, 2020, the description of exciton trapping by the RC and exciton recombination is only effective and relies on the results of more elaborate treatments performed in Refs. Kreisbeck et al., 2011; Shatokhin et al., 2018. There, it is realized that EET from one chromophore to another, as well as the radiative decay to the ground state, are actually mediated by the bath of environmental photons. Performing a second-order treatment of the appropriate interaction Hamiltonians, one ends up with effective Liouville superoperators $\mathcal{L}_{\mathrm{RC}}$ and $\mathcal{L}_{\mathrm{rec}}$ that describe respectively the excitation trapping and recombination on the level of reduced excitonic dynamics. In the following sections, we provide more details on the form of these effective Liouvillians and the manner in which they enter our description.

We use our exact description of weak-light-induced exciton dynamics in molecular aggregates, which is developed in the accompanying paper. Janković and Mančal (2020) The radiation correlation function, which is the only property of the radiation entering our reduced description, is modeled by the following expression Loudon (2000)

where $I_{0}$, $\omega_{c}$, and $\tau_{c}$ are the intensity, central frequency, and coherence time of the radiation, respectively. The radiation is assumed to have well defined directions of propagation and polarization. The weak-light assumption underlying our theoretical approach is well satisfied in a wide variety of photosynthetically relevant situations. For example, for a bacteriochlorophyll molecule irradiated by ambient sunlight, the excitation–light interaction may be estimated to be of the order of $10^{-3}\>\mathrm{cm}^{-1}$, see the accompanying paper Janković and Mančal (2020) and Ref. Chan et al., 2018. This energy scale is much smaller than the typical energy scales ($\sim 10-100\>\mathrm{cm}^{-1}$) of resonance couplings or the excitation–environment coupling. Also, the number of photons per unit time incident on a single photosynthetic complex under sunlight at the surface of the Earth can be estimated to be of the order of $1000\>\mathrm{s}^{-1}$, see the accompanying paper Janković and Mančal (2020) and Ref. Chuang and Brumer, 2020. In other words, the corresponding time scale is orders of magnitude longer than time scales typical for excitation transport, trapping at the load, and recombination.

For the sake of simplicity, we assume that the baths on both sites are identical, but uncorrelated. The bath correlation function, which is the only property of the bath entering the reduced description, can be decomposed into the optimized exponential series Zhang et al. (2018) ($t\geq 0$)

In Eq. (7), the collective bath coordinate $u_{j}$ is defined in Eq. (5), the expansion coefficients $c_{m}$ may be complex, while the corresponding decay rates $\mu_{m}$, as well as the white-noise-residue strength $\Delta$, are assumed to be real and positive. The bath correlation function is usually expressed in terms of the environmental spectral density $J(\omega)$ [$\beta=(k_{B}T)^{-1}$, where $T$ is the temperature]

which conveniently combines information on the density of environmental-mode states and the respective coupling strengths to electronic excitations. May and Kühn (2011); Valkunas, Abramavicius, and Mančal (2013) We explicitly treat only $K=N_{\mathrm{BE}}+N_{J}$ terms in Eq. (7), where $N_{\mathrm{BE}}$ and $N_{J}$ are the numbers of explicitly treated poles of the Bose–Einstein function and the bath spectral density, respectively. We assume the environmental spectral density of the overdamped Brownian oscillator

where $\lambda$ is the reorganization energy, while $\gamma^{-1}$ is the characteristic time scale for the decay of the bath correlation function $C(t)$.

The exponential decompositions embodied in Eqs. (6) and (7) enable us to formulate the problem as HEOM incorporating photoexcitation. As demonstrated in Ref. Janković and Mančal, 2020, the hierarchy consists of two parts, one in the $eg$ sector, and another in the $ee$ sector. Each density matrix $\sigma_{\mathbf{n}}(t)$ is uniquely characterized by vector $\mathbf{n}$ of non-negative integers $n_{j,m}$, where index $j$ enumerates chromophores, while index $m$ counts terms in the decomposition of $C(t)$ [Eq. (7)]. In order to describe excitation harvesting by the RC and recombination, we augment our formalism by effective Liouvillians describing these two processes. As demonstrated in Ref. Kreisbeck et al., 2011, these Liouvillians appear on each level of HEOM. Performing the appropriate rescalings, which ensure that auxiliary density operators are all dimensionless and consistently smaller in deeper levels of the hierarchy, Shi et al. (2009) we obtain the following equations describing the NESS we are interested in ($\gamma_{\mathbf{n}}=\sum_{j,m}n_{j,m}\mu_{m}$)

In the NESS, the continuity equation for exciton currents should be valid, i.e., the number of generated excitons per unit time must balance the sum of the recombination and trapping exciton fluxes. This is physically clear, and it is seen from a more formal perspective by taking the trace (with respect to the electronic system of interest) of Eq. (11) in which $\mathbf{n}=\mathbf{0}$. This results in

In the continuity equation [Eq. (12)], we define all currents to be positive, while the sign is determined by the ”direction” of the current ($+/-$ if it leads to an increase/a decrease in the exciton number). In more detail, the definitions of currents, which are dimensionless in our description, are

Let us immediately note that currents $J_{\mathrm{gen}},J_{\mathrm{RC}},$ and $J_{\mathrm{rec}}$ are written in a basis-invariant manner. This feature is quite appealing, since one can express currents in terms of populations and coherences in any particular basis. The light-harvesting efficiency is defined as

As discussed in the accompanying paper Janković and Mančal (2020) and in Sec. VI, for realistic values of the light coherence time $\tau_{c}\sim 1\>\mathrm{fs}$, the expression for the generation current may be further simplified to

In other words, possible enhancements in $\eta$ due to coherences in any basis ultimately originate from the expression for the trapping Liouvillian $\mathcal{L}_{\mathrm{RC}}$, as detailed in the following section.

In the literature on the physics of photosynthetic light harvesting, two bases play special roles. The first one is basis $\{|l_{j}\rangle|j\}$ of singly excited states localized on single chromophores (local or site basis). While Hamiltonian parameters are usually known in the local basis, the description of, e.g., absorption properties of a photosynthetic aggregate is usually performed in the excitonic basis $\{|x_{j}\rangle|j\}$, i.e., in the basis of stationary states of the isolated-aggregate Hamiltonian $H_{M}$. Claims about possible impact of coherences on the efficiency of photosynthetic light harvesting are usually made with the coherences in the excitonic or local basis in mind.

Tomasi and Kassal have recently classified different types of possible coherent enhancements of light harvesting efficiency according to the basis in which the excitation decay mechanisms are defined. Tomasi and Kassal (2020) Let us now focus on the trapping at the RC. Two forms for the Liouvillian $\mathcal{L}_{\mathrm{RC}}[\rho]$ are widely used in the literature. Jesenko and Žnidarič (2013); León-Montiel, Kassal, and Torres (2014); Tscherbul and Brumer (2018); Yang and Cao (2020); Jung and Brumer (2020) If site $j_{0}$ is closest to the RC, so that it is essentially the sole site coupled to it, one uses the so-called localized-trapping Liouvillian Jesenko and Žnidarič (2013); León-Montiel, Kassal, and Torres (2014); Tscherbul and Brumer (2018)

where $\tau_{\mathrm{RC}}$ is the characteristic time scale on which the populations are delivered from site $j_{0}$ to the RC. However, even in such a situation, the excitation transfer to the RC should be regarded as the multichromophoric Förster transfer, Jang, Newton, and Silbey (2004); Ma and Cao (2015) so that it is more appropriate to employ the so-called delocalized-trapping Liouvillian León-Montiel, Kassal, and Torres (2014); Tscherbul and Brumer (2018)

The recombination occurs due to both radiative and nonradiative processes. In other words, the lifetime of chromophores’ excited states is determined not only by the fluorescent decay, but to a great extent also by the conversion of the singlet states to triplets. While the radiative recombination is most naturally described in the excitonic basis, the nonradiative processes likely occur in a different basis of states, depending on how the triplet states are coupled among themselves. Overall, it is not possible to say that recombination occurs straightforwardly in neither the delocalized, nor local basis. For definiteness, we assume that the (nonradiative) recombination may occur from each site with the same rate constant $\tau_{\mathrm{rec}}$, and the appropriate Liouvillian reads as

The fact that we have chosen the local version of the recombination should not, however, influence our results, because the recombination process is generally much slower than every other process, it occurs on the time scale of nanoseconds, i.e., several orders of magnitude slower than trapping and other processes, see Sec. V. In situations in which the radiative decay is the major loss channel, it may be more appropriate to consider the recombination Liouvillian in the excitonic basis

where $\tau_{j}^{-1}$ is the Weisskopf–Wigner spontaneous emission rate Chan et al. (2018) from excitonic state $|x_{j}\rangle$.

If we assume the trapping at the RC to be governed by Eq. (18), the trapping current $J_{\mathrm{RC}}$ in Eq. (14) depends only on site populations, and on both exciton populations and interexciton coherences. If, on the other hand, we assume that the trapping is governed by Eq. (19), $J_{\mathrm{RC}}$ is expressed in terms of exciton populations only, while its expression in the local basis contains both site populations and intersite coherences. This has been recognized in recent studies, which ascertain that possible coherent enhancements of the efficiency can be achieved only when the coherence occurs in a basis different from that in which the trapping or recombination are modeled. Tscherbul and Brumer (2018); Tomasi and Kassal (2020)

However, we feel that the notion of coherent efficiency enhancements is not defined well enough. The word ”enhancement” would suggest that there is a reference value of the efficiency (which should be smaller than unity), with respect to which we can expect to achieve an enhancement. Moreover, since the enhancement is supposed to be coherent, one could imagine that the aforementioned reference value should depend on populations only, so that the subsequent inclusion of coherences should enhance the efficiency above that reference value. Whatever the form of the effective trapping and recombination Liouvillians is, there will always be a basis in which $J_{\mathrm{RC}}$ is entirely expressed in terms of basis-state populations only. Such a basis will be denoted as $\{|p_{j}\rangle|j\}$ and termed the preferred basis of the NESS under investigation. The trapping current is then expressed only in terms of the RDM diagonal elements, so that the efficiency value calculated in the preferred basis could be regarded as the reference value above which coherent enhancements due to non-zero values of coherences in some other basis (e.g., in the local or excitonic basis) may be possible. However, our definition of the efficiency [Eq. (16)] is basis-independent, and the fact that the coherences in the excitonic or local basis are non-zero should not be expected to bring about any efficiency enhancements.

While the above discussion about coherent efficiency enhancements is quite general, it will not be the main topic of our numerical investigations, which focus on a model photosynthetic dimer and in which the efficiency is close to unity, see Sec. V.1. Our central question is how the stationary state of the incoherently driven and loaded molecular system looks like, and how it relates to the two standard pictures, the picture of local states and the delocalized-states picture. Nevertheless, to the best of our knowledge, our work is among the first works that properly describe the principal physical features of excitation harvesting under incoherent light and properly define the light-harvesting efficiency. Chuang and Brumer (2020); Yang and Cao (2020); Jung and Brumer (2020) The arguments of the previous paragraph suggest that, once a proper description of the state in which the photosynthetic systems under incoherent illumination and load find themselves, the involvement of the coherences in the description of photosynthetic light harvesting should be critically reassessed. We thus believe that future applications of our NESS methodology to larger systems, see Sec. VI, could significantly contribute to the debate on possible coherent efficiency enhancements.

The excited-state sector of the steady-state RDM $\rho_{ee}^{ss}$ in the preferred basis reads as

The preferred basis is determined by the competition between excitation generation, pure dephasing, energy relaxation, excitation trapping at the RC and recombination.

The excitonic (or site) basis and the preferred basis of the NESS are connected through a unitary transformation. For our model dimer, the most general transformation of that kind can be parameterized by four real parameters, so that the basis vectors in the preferred basis are expressed in terms of the basis vectors in the excitonic basis as follows Murnaghan (1962)

Due to the phase freedom, we can immediately remove parameters $\varphi_{px}$ and $\psi_{px}$ from further discussion, so that we are left with only two parameters, $\theta_{px}$ and $\Delta_{px}$. From our subsequent discussion, it will emerge that the rotation angle $\theta_{px}$ is closely related to the analogous rotation angle which measures the deviation from the non-canonical statistics in the undriven and unloaded system (no generation, trapping, and recombination). Lee, Cao, and Gong (2012) The phase $\Delta_{px}$ is intimately connected to the rates of excitation trapping and recombination, which remove excitations from the system.

The parameters $\theta_{px}$ and $\Delta_{px}$ can be related to the Bloch angles $\theta_{B}^{x}$ and $\phi_{B}^{x}$ that are commonly used to characterize the basis in which the RDM is diagonal. Lee, Cao, and Gong (2012); Cerrillo and Cao (2014) Namely, one can always normalize $\rho_{ee}^{ss}$ to obtain $\widetilde{\rho}_{ee}^{ss}$ whose trace is unity and whose eigenvalues will be denoted as $\widetilde{p}_{j}$. For the model dimer, the operator $\widetilde{\rho}_{ee}^{ss}$ can always be expressed as

where $\mathbb{I}$ is the $2\times 2$ unity matrix, while $\bm{\sigma}=\{\sigma_{1},\sigma_{2},\sigma_{3}\}$ are three Pauli matrices. Let $\sigma_{3}$ be diagonal in the excitonic basis (i.e., $\sigma_{3}=|x_{0}\rangle\langle x_{0}|-|x_{1}\rangle\langle x_{1}|$) and let $\bm{a}^{x}$ be the corresponding vector in Eq. (24). Then it can be shown that the spherical angles $\theta_{B}^{x}$ and $\phi_{B}^{x}$ on the Bloch sphere are related to parameters $\theta_{px}$ and $\Delta_{px}$ as follows

Equations (25) and (26) relate $\theta_{px}$ and $\Delta_{px}$ to the Bloch angles $\theta_{B}^{x}$ and $\phi_{B}^{x}$ that are straightforwardly obtained from vector $\bm{a}^{x}$. In a similar manner, one obtains parameters $\theta_{pl}$ and $\Delta_{pl}$ of the unitary transformation that connects the preferred and the local basis. The rotation angle $\theta_{px}$ can always be chosen in the range $(0,\pi/4)$, while $\Delta_{px}\in(-\pi/4,\pi/4)$. This is discussed in greater detail in Sec. SI of the Supplementary Material.

So far, we have employed our novel theoretical approach to study in great detail the NESS of an incoherently driven and loaded dimer. The dimer represents the only model, in which one can represent the character of the NESS by only two parameters—the angle of the basis rotation $\theta$ with respect to a chosen basis, and the phase $\Delta$ that is closely connected with the excitation decay rates—and the difference between different NESSs is easy to understand. In other words, the dimer is the only model, where relatively simple “understanding” of the role of different parameters in establishing the preferred basis can be derived. Whether our results can be translated to larger systems is, of course, an important question. This section aims at presenting the basic steps that have to be taken in order to apply our NESS formalism to multichromophoric situations.

We start from the fact that the realistic coherence time of light $\tau_{c}\sim 1\>\mathrm{fs}$ is much shorter than any other time scale in the problem. All our results for the dimer are obtained by solving Eqs. (10) and (11) and do not lean on any assumption about the coherence time $\tau_{c}$ entering Eq. (6). In the accompanying paper, Janković and Mančal (2020) we have argued that, in the limit of small $\tau_{c}$, the first-order light correlation function defined in Eq. (6) may be replaced by

which represents the so-called white-noise model of the radiation. Olšina et al. (2014) Equation (10) is then omitted from further discussion, while the source term (the third term) of Eq. (11) reads as Janković and Mančal (2020)

We recall that the operator $\mu_{eg}=\mathbf{e}\cdot\bm{\mu}_{eg}$ is the projection of the dipole-moment operator $\bm{\mu}_{eg}$ onto the radiation polarization vector $\mathbf{e}$. Using the definition of $\bm{\mu}_{eg}$ in Eq. (4), the matrix elements of the source term in the basis $\{|b_{j}\rangle|j\}$ of the single-excitation manifold can be expressed as

The basis $\{|b_{j}\rangle|j\}$ is in principle arbitrary; it can be the local basis ($b=l$), the excitonic basis ($b=x$), the preferred basis ($b=p$), or any other basis in the single-excitation manifold. Equation (39) is in the form in which the rotational average can be straightforwardly performed with the final result Hein et al. (2012)

where $\mathbf{d}_{b_{k}}=\sum_{k^{\prime}}\mathbf{d}_{k^{\prime}}\langle b_{k}|l_{k^{\prime}}\rangle$ is the transition dipole moment of state $|b_{k}\rangle$. In a typical situation, the system of interest consists of many photosynthetic complexes in solution, and the rotational average is performed over random orientation of individual chromophores’ dipole moments with respect to the polarization direction. The source term in Eq. (40) depends only on relative orientations of transition dipole moments, which are known, e.g., for the widely investigated $Q_{y}$-band excitations of the FMO complex. Read et al. (2008) Representing Eq. (11) in basis $\{|b_{j}\rangle|j\}$ and using the rotationally averaged source term given in Eq. (40), we obtain a description of incoherent-light driven EET that exploits both the light incoherence (short $\tau_{c}$) and experimentally available data on relative orientations of transition dipole moments. Such a description features a much more realistic excitation condition than the one we have employed in the study of model dimer (the selective excitation of a local site).

Our NESS approach can also be used to follow the pathways of light-induced excitations, from the point of their generation, through the chromophore network, to the point of their extraction at the load. To elaborate on this, we compute the $|b_{k}\rangle\langle b_{k}|$ matrix element of Eq. (11) for the RDM ($\mathbf{n}=\mathbf{0}$) and obtain

In Eq. (41), $J^{b_{k}}_{\mathrm{gen}}$ is the excitation generation flux into the singly excited state $|b_{k}\rangle$

$J^{b_{k}}_{\mathrm{RC}}$ is the excitation trapping flux from $|b_{k}\rangle$

$J^{b_{k}}_{\mathrm{rec}}$ is the excitation recombination flux from $|b_{k}\rangle$

while $J^{b_{k}b_{k^{\prime}}}$ (for $k^{\prime}\neq k$) is the net flux of excitations that are exchanged between states $|b_{k}\rangle$ and $|b_{k^{\prime}}\rangle$

The last term in Eq. (41), $J^{b_{k}}_{\mathrm{res}}$, stems from the residual term $2\Delta\delta(t)$ in the decomposition of the bath correlation function into the optimized exponential series, see Eq. (7). This term can, therefore, be made arbitrarily small by explicitly treating a sufficient number $K$ of exponentially decaying terms in Eq. (7), and will not be considered in the further discussion. For the sake of completeness, let us also note that, if the light incoherence is exploited on the level of Eqs. (37)–(40), Eq. (42) for the generation flux should be replaced by