Dynamical transitions in Markovian exciton transport

The large deviation theory has recently been applied to open quantum systems to uncover dynamical crossovers in the space of quantum trajectories associated to Markovian evolutions. Such dynamical crossovers are characterized by qualitative changes in the fluctuations of rare quantum jump trajectories, and have been observed in the statistics of coherently driven quantum systems. In this article we investigate the counting statistics of rare quantum trajectories of an undriven system undergoing a pure relaxation process, namely, Markovian exciton transport. We find that dynamical crossovers occur in systems with a minimum of three interacting molecules, and are strongly activated by exciton delocalization and interference. Our results illustrate how quantum features of the underlying system Hamiltonian can influence the statistical properties of energy transfer processes in a non-trivial manner.


Introduction
Full counting statistics (FCS) is a powerful mathematical tool that enables the characterization of a counting process by systematic computation of the generator of its statistics. The technique was originally developed for theoretical modelling of mesoscopic charge transport [1,2], where current fluctuations can provide vital information on the nature of the transport mechanism [3][4][5][6][7][8]. More generally, FCS is of importance in the theory of open quantum systems as it provides an alternative characterization of the steady-state properties of a quantum system subjected to noise [9][10][11]. The applicability of FCS to open quantum systems stems from the realization that the dynamics of quantum systems in noisy environments can be described in terms of evolutions of ensembles of trajectories [12]. A trajectory, in this context, is a realization of the system's dynamics subjected to particular bath-induced jumps at stochastic intervals, and represents a sequence of configurations in the time domain. It is therefore the discrete nature of the quantum jumps that renders applicability to the tools of FCS.
In a recent contribution Garrahan and Lesanovsky [13] applied the mathematical framework of FCS to investigate the dynamical phases of quantum systems subjected to Markovian noise. FCS states that a knowledge of the generating function enables complete characterization of the statistics. The contribution of Garrahan and Lesanovsky [13] was to reinterpret the counting field of FCS as a means by which sub-ensembles of rare trajectories can be accessed, and to investigate how the statistical generator is changed as a function of the counting field. For sub-ensembles of rare trajectories the frequency of the number of jumps after a long observation time, a quantity known as the activity, deviates from the jump rate in the steady-state, i.e. from the jump rate of the entire ensemble of trajectories. By unravelling the trajectory space of a master equation and identifying the rare sub-ensembles, one can uncover subtle dynamical features, such as dynamical transitions, that cannot be determined from the static properties, or the ensemble averaged dynamics [14]. Dynamical transitions occur if sub-ensembles of trajectories with different activities demonstrate qualitatively different fluctuations. They have so far been predicted in photon counting [13,[15][16][17], electron counting [18], quadrature measurements [19] and dissipative quantum walks [20], suggesting the possibility of their existence in a host of other systems.
The formalism developed by Garrahan and Lesanovsky is known as the s-ensemble theory. The parameter s has the same mathematical structure as the counting field of FCS; it also has the added feature that s = 0 allows access to sub-ensembles of rare trajectories. In some applications of the s-ensemble theory rare trajectories are referred to as the 'non-equilibrium' trajectories to distinguish them from the typical trajectories. This terminology can lead to confusion as in the theory of open quantum systems, a non-equilibrium state is characterized by a finite rate of entropy production [21]. In the absence of external drives, rare trajectories, or typical ones, do not produce entropy. In this paper we are interested in dynamical features of quantum systems at thermal equilibrium and we will therefore avoid this terminology.
In this paper we apply the s-ensemble theory to investigate the statistics of rare trajectories of quantum evolutions describing excitation energy transport in a molecular complex. Assuming strong electronic interactions, and weak system-environment coupling, exciton relaxation can be described by Lindblad type evolution [22]. The exciton relaxation process results in incoherent transfer between electronic eigenstates, reaching a steady-state where the average population of the eigenstates obeys the thermal Boltzmann distribution. At a microscopic level, however, transfer of excitation still occurs in a dynamical manner in the steady-state. For instance for a system of two eigenstates at thermal equilibrium, detailed balance dictates that k 1→2 P 1 = k 2→1 P 2 where k i→ j is the rate of energy transfer from the state i to the state j and P i is the thermally equilibrated population of the state i. We demonstrate that, in this system, despite the classical nature of the dynamics, the statistics of the trajectories are not, in general, Poissonian. We observe temperature-dependent dynamical crossovers for systems containing a minimum of three sites. Statistical crossovers are found to be strongly influenced by exciton interference.
This paper is organized as follows: in section 2 we introduce the conceptual and mathematical framework of the large deviation principle; the mathematical framework behind the s-ensemble theory. In section 3 the Markovian Lindblad master equation and the s-ensemble theory are introduced. In section 4 we apply the s-ensemble theory to investigate the statistics of rare trajectories in systems of two, three, four and seven sites and present numerical results that demonstrate the dynamical transitions in the statistics of trajectories. The electronic and vibrational parameters of the Fenna-Matthews-Olson (FMO) photosynthetic complex are used in the computation of the numerical results, as the Lindblad master equation considered in this paper has been found to make a reasonable prediction of the steady-state of the complex. Section 6 summarizes the main conclusions of the paper.

Large deviation principle and counting statistics
A counting process can be decomposed into a number of trajectories where in each trajectory the process occurs at different stochastic intervals [23]. In the context of quantum jumps, the counting process is the number of jumps observed after a given time t. Let us denote the probability of observing K jumps after a time t by P t (K ), where ∞ K =0 P t (K ) = 1 and {0 P t (K ) 1}. For a counting process that satisfies the large deviation theory, P t (K ) acquires a large deviation form for large t. That is where the function ϕ(k) is known as the rate function [24] and k = K /t is the jump rate. For a simple (unimodal) counting process, ϕ(k) has a local minima corresponding to the most probable jump rate after a long observation time. The generating function of a given probability distribution is defined to be where s is a conjugate variable to K , and is known as the counting field in FCS. If the counting process has a large deviation form, the generating function can be written as Z t (s) is now referred to as the large deviation function, while θ(s) is the cumulant generating function. The cumulant generating function and the rate function are related via a Legendre transform The average and the standard deviation of the number of jumps in the steady-state can be deduced from the cumulant generating function, and are given by Similarly, higher derivatives of θ (s) evaluated at s = 0 are the higher order statistical cumulants. FCS (or the steady-state statistics) is therefore recovered by looking at the derivatives of the cumulant generating function at s = 0. The applicability of the large deviation approach to the identification of dynamical phases becomes evident when the behaviour of the cumulant generating function away from s = 0 is considered. To illustrate how sub-ensembles of rare trajectories can be accessed via the exponential factor e −s K , let us construct an associated stochastic process defined by the probabilities [25] q t (K , s) For s = 0 we recover the original counting process. For s < 0 (s > 0) the exponential factor amplifies (lowers) the probabilities of rare events. The coupling field s, thus provides a systematic means by which rare events can be accessed: unlikely events of the counting process {P t (K )}, correspond to typical events of the associated process {q t (K , s)}. Rare events can thus be accessed without generating a large number of typical events and relying on the statistics of the original counting process. In statistical mechanics the counting process is the number of microstate configurations accessible to the system under certain constraints, making Z t (s) the partition function. Analogously, for an open quantum system, Z t (s) is a measure of the number of trajectories that can participate in a given (rare) sub-ensemble. One can also identify ϕ(k) as the entropy: extrema of ϕ(k) therefore correspond to the typical trajectories. The function θ(s) is the equivalent of a free energy and its derivative with respect to s, is a dynamical order parameter known as the activity. Non-analyticities in free energy lead to abrupt jumps in the order parameter. These jumps correspond to first-order thermodynamic phase transitions, or (firstorder) dynamical transitions in the space of trajectories of a quantum system.
First-order dynamical transitions at s = 0 have important consequences for the steady-state dynamics, and indicate the coexistence of dynamical phases, characterized by distinct activities. The implications of transitions at s = 0 are not always clear and tend to be problem-specific.
Sometimes the original master equation can be mapped into another master equation, with different jump operators, that undergoes dynamical transitions at s = 0 [13]. In this paper we endeavour to explore the possibility of dynamical transitions in electronic energy transfer (EET) as the size of the system is varied, and will not explore such mappings.

s-ensemble theory and electronic energy transfer
In this section we apply the large deviation theory to investigate the trajectories of a particular open quantum system, modelled by the Markovian Lindblad master equation. The system of interest consists of a number of molecules (sites) coupled via electronic couplings and weakly interacting with their environment. Each molecule is modelled as a two-level system in its electronic coordinate, where the two levels correspond to the ground and the first excited electronic states. The electronic degrees of freedom are assumed to be diagonally coupled to a bosonic bath of harmonic oscillators. We assume that there is one excitation in the system that can be coherently shared between multiple molecules (a molecular exciton). The Hamiltonian of this system is given by where J i j is the electronic coupling between the ground state of site i and excited state of site j, the ket |i represents the state in which site i is excited and all other states are in the ground state, h ik is the electron-phonon coupling between site i and mode k of the bath, {b † k , b k } are the raising and lowering operators for mode k of the bath and ω k are the bath frequencies.
Due to the assumption of weak system-bath interaction, the exciton states are delocalized and the influence of the environment is to induce jumps between them at stochastic intervals. In the steady-state, the average population of the exciton states is given by the Boltzmann factor: P i = e −βω i /Z , where P i is the population of the exciton state with energy ω i and Z is the partition function. From a microscopic standpoint, dynamic transfer of excitation still persists in the steady-state. This microscopic transfer consists of incoherent jumps between the excitonic eigenstates. For instance for a system of two sites, steady-state implies that k 1→2 P 1 = k 2→1 P 2 , where k i→ j is the transfer rate from the eigenstate i to the eigenstate j. Moreover, the ratio of the transfer rates satisfies detailed balance: k 1→2 = k 2→1 e −βω 12 , where ω 12 is the energy difference between the two states. k 1→2 P 1 is also the average number of jumps per second between the two states: upward and downward jumps occur with equal frequency in the steady-state. For rare sub-ensembles of trajectories, the frequency of the jumps deviates from the equilibrium rate k 1→2 P 1 .
To compute the statistics of jumps in the steady-state, one can compute the cumulant generating function of this process directly from the Lindblad equation. The Lindblad equation as applied to EET, for weak system-bath coupling with secular approximation, is given by [22] where the commutator describes the unitary evolution and the second term describes relaxation and dephasing in the eigenstates {|α } of the unperturbed system Hamiltonian where {|m } are the localized wavefunctions. For two eigenstates |α and |α , we define the 'intensity factor' to be The intensity factor adds up the squared magnitude of the exciton interference at all sites and determines the contribution of the electronic Hamiltonian to the transfer rate between two exciton states. The factor γ (ω) in equation (9) determines the bath contribution to the exciton transfer rate and is related to the spectral density of the bath J (ω), via the expression γ (ω) = 2π J (|ω|)|n(ω)|, where n(ω) = (e βω − 1) −1 is the thermal occupation number. In the present treatment it is assumed that each molecule interacts with an identical and independent phonon bath. The dissipative term in the Lindblad equation describes the incoherent jumps between the exciton states as induced by the phonon bath. In the quantum jump approach the Hamiltonian evolution of the quantum system is interrupted by jumps at stochastic intervals. A time-record of deterministic (non-unitary) evolution of the wavefunction interrupted by stochastic jumps constitutes a trajectory of the system. An averaging over a large number of such trajectories recovers the reduced density matrix σ (t) [12].
If K jumps are observed after a time t, the (unnormalized) reduced density matrix of the system prior to detection is given by the projection of the total density matrix onto the subspace of K events [9,10]. We denote this projected density matrix by σ (K ) (t). Conversely, the total density matrix can be written as an expansion of all projected density matrices, that is σ (t) = ∞ K =0 σ (K ) (t). The probability of observing K events after a time t is given by P t (K ) =Tr[σ (K ) (t)]. At this point for the sake of clarity we restrict the analysis to the statistics of jumps between two exciton states separated in energy by ω 1 . Note that the total system may still contain an arbitrary number of exciton states. It can be shown that the Laplace transform of the projected density operator, defined by the expression σ s (t) = ∞ 0 σ (K ) e −s K obeys the equationσ Note that σ s (t) has the same structure as the variable q t (K , s) in equation (7). Typical trajectories of σ s (t) thus correspond to rare trajectories of the original master equation. The ensemble of trajectories generated by equation (12) is known as the s-ensemble [13]. Physical time evolution is retrieved at s = 0, whereas s = 0 contains information with regards to rare trajectories of the system. s < 0 (s > 0) generates the active (inactive) trajectories where the number of jumps are significantly higher (lower) than average.
The Markovian nature of the Lindblad equation enables us to compute the statistical moments of the number of jumps in the steady-state via the large deviation principle. Large deviation theory states that for a Markov process, the largest real eigenvalue of the superoperator W s is the cumulant generating function of the statistics [24]. Following the notation of the previous section, we denote this generator by θ(s). In the next section we compute the cumulant generating function and the associated statistical moments for a system of two, three, four and seven sites.

Results and discussion
Adopting the electronic and vibrational parameters of the FMO complex, we first consider two and three sites (chromophores) and study the statistics of trajectories that contain jumps between two given exciton states. We subsequently include further sites in the system and study the statistics of selected jumps. We describe the coupling of the electronic system to the environment via the spectral density [26] where E r = 35 cm −1 is the reorganization energy E r = ∞ 0 J (ω) ω dω, and ω c = 150 cm −1 is the cutoff frequency of the bath [27,28]. This form of the spectral density describes the low frequency contribution of the environmental fluctuation spectra correctly and captures the main features of the dynamics. The sensitivity of the dynamical transitions on the form of the spectral density is an interesting topic that will not be addressed in this work. We compute the matrix W s and solve for the cumulant generating function θ(s) . Figures 1(a) and (b) illustrate the exciton states where two and three sites are included in the dynamics, respectively. The two site model includes sites 1 and 2, and the three site model includes sites 1-3. The site indices are the standard indices assigned to the chromophores in the FMO complex [29]. The exciton states are labelled via a Greek letter and a number, where the Greek letter represents the position of the exciton state in the energy ladder, and the number indicates the site that carries the maximum exciton amplitude. For instance |α : 1 represents the lowest eigenstate, and indicates that this eigenstate is predominately localized on site 1. Due to significant exciton localization in the FMO protein, each exciton can be unambiguously assigned to a unique site index. Figures 1(c) and (d) illustrate the cumulant generating function versus s for different temperatures. The plots suggest that the two cases have similar statistics in the active phase, but begin to diverge close to s = 0. Figure 1(e) illustrates the Mandel parameter for the case of two sites. The Mandel parameter is an indicator of the nature of first-order statistical correlations in the system and is defined to be  [28]. maximum in the Mandel parameter. The maximum is indicative of coexistence of two dynamical phases characterized by distinct activities. In particular, if a maximum occurs at s = 0, there is coexistence of dynamical phases in the steady-state. For the case of two sites, the Mandel parameter is negative for all s. This indicates sub-Poissonian correlations (anti-bunching) where a waiting period for re-excitation is required before the system can relax back to the lower state. As s is increased, the Mandel parameter approaches zero (Poissonian limit). In this limit there are so few jumps that correlations vanish asymptotically.
It is striking that anti-bunching is observed even though the transitions are thermally activated, and the steady-state dynamics obey a classical rate equation. The classical nature of the correlations can be confirmed by considering the rate equation characterized by the following dynamical matrix: where κ and are the relevant equilibrium transfer rates and = κe −βω as required by detailed balance. By diagonalizing W class (s), we arrive at the following expression for Q(s): Indeed, the same expression for Q(s) is obtained if the full Lindblad equation with dephasing terms is instead considered. The corresponding plot of the Mandel parameter for a system of three sites is shown in figure 1(f). The trajectories exhibit a dynamical crossover in their statistics, i.e. the correlations change from sub-Poissonian (Q(s) < 0) to super-Poissonian (Q(s) > 0) as s is varied. The crossover point is temperature dependent and is closer to s = 0 at higher temperatures.
One can continue to add more sites to the model and investigate the statistics of a particular jump, or the collective statistics of a number of jumps. For instance the collective statistics of all jumps to the lowest eigenstate can be investigated. In general, beyond a few sites there are many possible jumps and there is no obvious advantageous collective counting strategy. The statistics of selected transitions for the four site and the seven site models, are explored in the remaining part of this section.
Four interacting sites exhibit the combination of the features of the two site and the three site models, as illustrated in figure 2. Figure 2(b) exhibits the new feature that typical trajectories (s = 0) have different statistics at different temperatures: anti-bunching at higher temperatures and bunching at lower temperatures. In figure 2(c) the trajectories remain anti-bunched over the computed temperature range, similar to the results of the two site model.
We consider all possible pairwise jumps of the four site model and conclude that for sufficiently small s all trajectories become anti-bunched and approach a temperature independent but transition dependent value of the Mandel parameter. For sufficiently large s, on the other hand, all trajectories approach the Poissonian limit. Trajectories that undergo a statistical crossover, approach the Poissonian limit from the positive values of Q(s); those that do not, remain in the Q(s) < 0 portion. The approach to the Poissonian limit from the Q(s) > 0 side is therefore a sufficient indicator of the existence of statistical crossovers. These conclusions are also observed to be valid for larger systems.
Finally, statistics of selected jumps from the seven site system of the FMO complex are illustrated in figure 3. All plots reveal statistical crossovers at certain temperatures. The crossovers are most likely to occur in the unravelling of the jumps associated to fastest exciton transfer rates. The transfer rate between two exciton states separated in energy by ω is given by (ω) = γ (ω) m |c m (α)| 2 |c m (α )| 2 . In this system, the relative magnitude of the transfer rates tends to be dominated by the intensity factor m |c m (α)| 2 |c m (α )| 2 , as the variation of γ (ω) over the energy range of interest is small. In other words, dynamical crossovers are observed in trajectories accounting for jumps between eigenstates producing the largest intensity factor. The statistics of the rare trajectories are therefore sensitively dependent on the distribution of the electronic couplings and the degree of delocalization and interference. For exciton states that  [28].
are close in energy, γ (ω) is large, resulting in fast transfer. Statistical crossovers are therefore also observed in this regime, despite the small intensity factor. An example of this behaviour can be seen in figure 3(a). We conclude by noting that the importance of the Markovian assumption in the present treatment poses a difficulty for many systems exhibiting EET whose environmental interactions are strongly non-Markovian. The formalism can, however, be applied to model certain types of non-Markovianity by extending the Hilbert space of the system Hamiltonian to contain the environmental modes that capture the strong non-Markovian features [30][31][32].

Experimental plausibility
Unlike the familiar problems of photon counting in atomic systems, the transitions discussed in this paper are not accompanied by absorption and emission of real photons. This means that unless the energy exchange with the bath is somehow monitored, the trajectories cannot be directly probed. Monitoring the energetic fluctuations of the bath is certainly not a straightforward task, there is thus an experimental difficulty associated with the direct detection of the dynamical crossovers discussed in this paper. However, it might be possible to devise an experimental scheme that aims to identify indirect experimental signatures of dynamical crossovers. To illustrate this idea consider an external parameter F, that influences the activity of the trajectories and can be varied at will in an experiment (e.g. temperature of the bath, separation between the sites, etc). The task is then to measure the transfer rates between all exciton pairs as this external parameter is varied. From these measurements one can construct Q(F) at s = 0, from the corresponding analytic expression for Q(s = 0). If Q(F) has a local maximum at a given F, we anticipate dynamical phase coexistence in the vicinity of this point. In other words, maxima of Q as a function of an external parameter are likely to indicate the mixing of two dynamical phases [14], and thus the existence of dynamical crossovers. This connection can be theoretically verified by computing the probability distribution of the activity order parameter as discussed in [14].

Conclusion
Markovian master equations are an established method of modelling the relaxation dynamics of many systems exhibiting transport properties. These equations are primarily used to compute the ensemble averages, but in fact contain more information pertaining to higher order statistical moments that cannot be gleaned unless an unravelling formalism is first implemented. By investigating the properties of the statistical generator as a function of the parameter s, the s-ensemble theory extends the unravelling formalism to sub-ensembles of rare trajectories of a quantum process. This enables the identification of dynamical phases, characterized by distinct statistical features. Such an unravelling formalism has been applied in this paper to investigate the statistics of Markovian excitation transport.
The steady-state of a system undergoing EET exhibits surprisingly rich statistics: the features that are often associated with quantum correlations such as anti-bunching can be observed in the steady-state statistics of the thermally activated quantum jump trajectories. System size also plays a key role in the statistics of the dynamical phases as systems of two and three sites reveal qualitatively different correlations: while in the former the jumps are always sub-Poissonian, in the latter they can undergo statistical crossovers. Transitions in larger systems show a combination of these two characteristic behaviours. Importantly, statistical crossovers are found to be activated by electronic interferences: trajectories that account for jumps between two strongly interfering eigenstates (as defined in section 3) are most likely to exhibit statistical crossovers. The structure of the spectral density can also influence the statistics by promoting electronic transitions of particular frequencies.
The physical implications of dynamical transitions away from s = 0 are not always clear and are the subject of debate. For instance, it has been shown that by rescaling time, it is possible to map the dynamics to another system, such that the rare trajectories of the original system are the typical trajectories of the new system [13]. The new system then exhibits dynamical transitions under physical time evolution. It would be interesting to explore whether such mappings can be obtained for the present system.
We leave the matter of the potential implications of the observed statistical transitions to excitation transport as an intriguing open question. As spectroscopy is the established tool of probing EET, one is faced with the question: does a statistical crossover at s = 0 leave its signature on any of the spectral features of a molecular complex? This question is the subject of our current investigation.