{'fraction_non_alphanumeric': 0.03273315717241238, 'fraction_numerical': 0.010221430180872201, 'mean_word_length': 4.842508377213978, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Coherent quantum dynamics launched by incoherent relaxation in a quantum circuit simulator of a light-harvesting complex

Engineering and harnessing coherent excitonic transport in organic nanostructures has recently been suggested as a promising way towards improving man-made light harvesting materials. However, realising and testing the dissipative system-environment models underlying these proposals is presently very challenging in supramolecular materials. A promising alternative is to use simpler and highly tunable quantum simulators built from programmable qubits, as recently achieved in a superconducting circuit by Potocnik et al. In this article, we simulate the real-time dynamics of an exciton coupled to a quantum bath as it moves through a network based on the quantum circuit of Potocnik et al. Using the numerically exact hierarchical equations of motion to capture the open quantum system dynamics, we find that an ultrafast but completely incoherent relaxation from a high-lying bright exciton into a doublet of closely spaced dark excitons can spontaneously generate electronic coherences and oscillatory real-space motion across the network (quantum beats). Importantly, we show that this behaviour also survives when the environmental noise is classically stochastic (effectively high temperature), as in present experiments, and also leads to a novel, transient violation of detailed balance in the population relaxation. These predictions highlight the possibilities of designing matched electronic and spectral noise structures for robust coherence generation that does not require coherent excitation or cold environments.


I. INTRODUCTION
Creating and sustaining 'long-lived' electronic coherences in complex, multi-component supramolecular systems has recently been highlighted as an exciting route towards advanced molecular nanodevices with applications ranging from energy harvesting to optomechanics and sensing 2, 3 . In this context, 'long-lived' refers to decoherence times of comparable duration to the 'functional' timescales of the system, which might, for example, correspond to energy transport times or charge generation, in the case of photovoltaic light-harvesting structures [4][5][6][7][8][9][10][11][12] . However, many -if not all -reliable molecular functions are driven in a thermodynamic direction by noisy interactions between electronic degrees of freedom and their thermal environments, so complete suppression of environmental couplings -as is desirable for quantum computation -is not a fruitful strategy for the multitude of applications discussed in Refs 2, 3 .
Instead, there has been an emerging interdisciplinary focus on understanding how it may be possible to exploit the non-perturbative and non-Markovian dynamics of structured system-environment interactions in nanostructured systems 4,10,[13][14][15][16][17][18][19][20][21][22][23][24] , with mounting theoretical evidence that a transient and correlated interplay of dissipative and coherent dynamics may be advantageous for a wide range of ultrafast optoelectronic processes. Indeed, although this essential idea has an origin in studies of photosynthetic pigment-protein complexes (PPCs), it is in rationally designed, organic functional materials, such as DNA origami, polymer-fullerene heterojunctions, carbon nano tubes and molecular dimer systems, that the existence and potentially beneficial impacts of electronic coherence and 'noise-assisted' dynamics on lightharvesting processes have been most cleanly and recently demonstrated [25][26][27][28][29][30][31] . Examples of theoretically proposed 'noise-asisted' quantum phenomena and their potential applications are reviewed in Refs 2, 3,9 . Regardless of whether Nature got there first, or at all 32 , these latter studies underscore the new possibilities arising from exploiting emerging nanofabrication techniques to tune both the properties of photo-excited states (delocalisation, dipoles, energy spectrum) and their environments to obtain novel optoelectronic materials based on taylored system-environment interactions. Recently, Potočnik et al. have demonstrated the first experimental 'quantum simulator' of an open quantum light harvesting model built from transmon qubits in a superconducting circuit (Fig.1) 1 . Using three individually tunable qubits coupled to a transmission line (for photoexcitation) and a resonator (to detect emission), Potočnik et al. demonstrated the formation of robust, delocalised photoexcited states with optical properties analogous to the Davydov-split (Frenkel) excitonic states found in PPCs or J-aggregates 33,34 . These states are engineered, as in photosynthetic antennae complexes, so that energy absorbed by the highest energy state is spatially directed by dissipation towards the lowest energy state 34,35 , which is proximate to a 'reaction centre' that transduces this incoming energy (here, the resonator). Although impos-sible in any real supramolecular structure, this set up also allows controllable application of environmental dephasing noise of arbitrary strength and spectral properties, which in the basis of delocalised states (vide infra) leads to controllable incoherent transitions between the single particle excited states of the network. This feature of the experiment makes it a near-ideal platform for testing theories of open dynamics, and by varying the noise coupling strength, it was demonstrated that the energy transfered across the network was maximised at an optimal value of the dephasing rate, precisely as predicted by recent theories of 'noise-assisted transport' (also known as 'ENAQT') [36][37][38][39] . In further agreement, this optimal dephasing noise strength was found to be of similar magnitude to the smallest coherent coupling between the qubits, occuring at the 'strong-to-weak coupling' transition point where the lowest energy delocalised eigenstates begin to collapse into effectively localised on-site excitations with sequential hopping transport 1 .
Additionally, they also confirmed that energy transfer is considerably more efficient when the spectral function of the environment is strongly peaked around the energy differences between the excitonic excited states. Such structured environmental spectral functions are characterisitic of molecular vibrations, and have become intensively studied in open quantum system theory due to their multi-scale linear response functions (memory effects) and often non-perturbative coupling to the excited states at specific frequencies. These properties allow potentially qualitative and non-stationary modifications of excited state dynamics, c.f. simple heat baths, that have been connected to (transient) phenomena such as violation of detailed balance, extension of electronic coherence times and vibronic mixing of electronic states 4,8,16,18,20,21,23,24,28,[40][41][42] . However, to describe the dissipative dynamics of systems coupled to such environments requires advanced numerical and analytical techniques, and approaches ranging from manybody methods to advanced master equation formulations have recently been applied or developed for this aim 13,19,[43][44][45][46][47][48][49][50][51][52][53][54][55][56] .
In this article, we explore the real-time dynamics of the three-qubit model implemented by Potočnik et al, using the numerically exact hierarchical equations of motion (HEOM) technique 55,57-62 to address a number of theoretical questions that could be verified in a future time-resolved version of the experiment. Specifically we shall show that the set up of Ref. 1 is an ideal platform to demonstrate the generation of coherence by incoherent processes, in this case caused by the noise-induced relaxation of a high energy state into a closely spaced doublet of lower energy states. This is a timely topic, as most observations of coherent optoelectronic phenomena only appear under coherent excitation, whereas functional light harvesting devices are likely to operate under incoherent illumination, i.e. sunlight 63,64 . Without access to excitation sources that can generate electronic coherencs, i.e. laser pulses, future coherent devices must rely on transient internal or inter-component dynamics to induce wave-like phenomena. The present work offers some insight into how this might be obtained from an engineering of electronic eigenstates to match a structured bath spectral density, highlighting the importance of different energy scales in the problem. We also note that another type of quantum simulator for light-harvesting using trapped ions has also recently been demonstrated 65 .
Specifically, we use HEOM to prove that although dissipative 'population-to-coherence' processes are nonsecular in nature, these often neglected transitions in the perturbative Redfield theory can, following decay of the high energy state, generate long-lasting coherence between the lower energy eigenstates. We demonstrate that these quantum beats are a manifestation of real-space coherent motions that could be detectable in a superconducting circuit experiment. To make further connections to general experimental conditions, we also show that interactions with both quantum and classically stochastic fluctuating environments can generate these coherent dynamics, and find that there is an optimal coupling for coherence generation that lies in an intermediate coupling regime. The possible role of non-Markovianity in these phenomena is also studied, as the structured spectral density we consider has a longer correlation time than the incoherent decay dynamics. However, while we explicitly demonstrate that the use of HEOM is essential to account for the strong environmental memory effects (especially for classical noise), the correlation betwen the generation of coherence and the formal measure of non-Markovianity we use appears to be weak.
The paper is organized as follows. Section II describes the three-state model and the system-bath interaction. Qualitative predictions and concepts emerging from Redfield approach are given in section III with a presentation of the operational HEOM equations. Section IV presents our numerical results for a quantum or classical noise and Section V provides some discussion and perspectives for future investigations.

A. Electronic system
A virtue of the model 'excitonic' Hamiltonian implemented in Ref. 1 is that its dynamics, with typical timescales of µs in superconducting circuits, remains unchanged when all energetic parameters are rescaled to optical frequencies (factor ≈ 10 5 ). In light of this, and in the interest of understanding coherence generation in molecular systems, we will retain the relative energy level and coupling structure of the qubits in Ref. 1 , but work at the time/frequency scales of molecular optics and replace qubits with 'sites' representating chromophores.
The "electronic" system consists of three chromophoric sites (two-level systems) but the active Hilbert space is confined to the single excitation sector and therefore is  described by three states |n (n = 1, 2, 3) corresponding to a single localised excitation on site n. Following the qubit topology of Fig. 1, the model Hamiltonian for the chromophore system is given by As in the experiment, the first two states (|1 and |2 ) are tuned to degeneracy ( 1 = 2 = 0) and strongly coupled by a coherent coupling J 12 while state |2 is weakly coupled to a third state |3 , through J 23 = J 12 /10. J 13 = 0, which is an approximation very close to the experimental realisation (qubit 1 and qubit 3 are not physically close to each other -see Fig.1a. The energy gap between the degenerate levels and the lower state energy is equal to the first coupling 3 = −J 12 . Diagonalising this simple Hamiltonian then leads to the eigenstate spectrum shown in Fig. 1c, which is characterised by a single high energy state and a low-lying doublet of states with an energy splitting ≈ 2J 12 , approximately ten times smaller than the mean doubletto-high-energy state gap. Due to this tuning of states, the eigenstates are highly delocalised over the sites. The high energy 'bright' state is approximately given by |B = 1 √ 2 (|1 + |2 ), whereas the lower-energy 'dark' states are given by |D ± = 1 2 (|1 − |2 ) ± 1 √ 2 |3 . We remark that, here, 'bright' and 'dark' refers to the coupling of these states to the transmission line. As seen in Fig.1a, sites 1&2 are close to the waveguide and both couple to the excitation field with the same strength.
Consequently, the transition dipoles of these sites interfere constructively in the |B eigenstate, making this 'bright', while destructive interference decouples the two dark states |D ± from the excitation fields. Experimentally, this configuration is very useful, as it allows optical population of a single, well-defined state from which transport then proceeds, while the non-emissive nature of the dark states prevents radiative losses and noise which might obscure the signatures of energy flow across the site network. Beyond practical considerations, it has also been proposed that using such dark states to reduce radiative losses could boost the efficiency of 'bio-inspired' organic photovolatics devices 20,40,[66][67][68][69][70][71] . Finally, the resonator is only coupled to state |3 , so only |D ± will emit into this channel. Experimentally, it is this resonator emission that is used to quantify the energy transfer from the |B state. For clarity of discussion, we will not explicitly model the coupling of the system to the excitation and read-out fields, but instead consider dynamics beginning with a population prepared in the |B state. In the quantum simulator, noise is generated by applying a stochastic signal along the flux lines that are used to tune the excited state energies of the individual qubits. This effectively introduces site-selective, stochastic (gaussian) noise that is diagonal in the basis of the localised qubit excitations (pure dephasing noise), but this noise is essentially classical in nature (vide infra and see section IV B). Nevertheless, this noise is generated by a signal generator that can produce almost arbitrary stochastic power spectra (see section III A), providing a versatile tool for probing dissipative quantum transport. In order to connect with molecular systems, we will consider both classical and quantum noise within a common framework in which the environment is treated as a continuum of quantum harmonic oscillators coupled linearly to the electronic system. Following Ref. 1 , we will consider the case of noise applied to only one site of the network, site 2, so that the quantum system-bath coupling is given by where S = |2 2| is the system operator and the displacement enviromental operator is defined by where a † k , a k are the bosonic creation and annihilation operators, respectively, of an oscillator of frequency ω k that is coupled to state |2 with amplitude g k . The Hamiltonian of the oscillator bath is k a k (with = 1). The system bath coupling also leads to an energy shift λ = −1/2 k g 2 k /ω k of state |2 (the reorganisation energy) that is added to the system Hamiltonian to define an effective system Hamiltonian H S,ef f = H S + |2 2|λ. The total systemenvironment Hamiltonian becomes H = H S,ef f + H SB + H B . By diagonalising H S,ef f , all the eigen states are coupled through the environment by off diagonal terms of the system coupling operator in the eigenbasis set The reorganisation energy is an indicator of the coupling strength since it modifies the system coupling operator and the energy gap among the eigenstates. We shall use a dimensionless parameter giving the ratio between the renormalisation energy and the dissipation free E BD+ energy gap 56 . Prior to excitation of state |B , we will always assume the environment osillators to be in thermal equilibrium w.r.t their free Hamiltonian H B . With this common asumption, the behaviour of the reduced density matrix of the electronic subsystem is completely characterised by the environment's spectral function J(ω) = k (g 2 k /ω k )δ(ω − ω k ) and its temperature through the Bose function n(ω) = (e βω − 1) −1 where β = 1/k B T and k B is the Boltzmann constant. Both appear in the thermal two-time correlation function of the oscillator displacements, which ultimately determines the dissipative physics of the system (see below) 72 . The correlation function is given by Heisenberg picture w.r.t. the Hamiltonian H B of the environmental displacement operator X. This leads to In the context of open system theory, the difference between quantum and classical stochastic noise is most clearly seen in C(t); the two-time correlation function for non-commuting operators is complex-valued, while for a classical scalar variable it is real. From Eq. 5 and the fact that J(−ω) = −J(ω), it can be seen that C(t) becomes real in the limit of β → 0 (high temperature limit) and we shall later exploit this to extract results about classical noise from our HEOM simulations. This is further illustrated in Fig. 2 for the structured spectral density that we will consider in our numerical results. The spectral density is here a superohmic Lorentzian function where The parameters are chosen to create a sharp spectral density peaked at the energy gap E BD , as shown in the inset of Fig. 2. The numerical values of the parameterisation used are given in Appendix B, with the p parameter being used as a scaling factor that allows us to vary the reorganisation energy of the bath.

A. Redfield Equations
For open quantum systems, the principal object of interest is the reduced density matrix of the system ρ S (t) = Tr B [ρ(t)], where ρ(t) is the joint density matrix of the system-environment state. In general, determining ρ(t) in order to obtain ρ S (t) is highly demanding, however it is possible in many cases to find approximations that greatly simplify this task, and which also provide very useful concepts and intuition for discussing more complex open physics.
For the case of weak coupling (second order perturbation theory w.r.t. system-bath coupling) leading to a broad spectral function characterised by a cutoff frequency ω c that is much larger than any energy difference in the system Hamiltonian H S,ef f , the dynamics of ρ S (t) are often well-described by the Bloch-Redfield theory. Detailed derivations of the Bloch-Redfield master equation can be found in Refs. 72,73 , here we shall simply state the results of immediate consequence for our results and discussion. Following the Born-Markov approximation, the equation of motion for the reduced density matrix in the interaction pictureρ S (t) = e −iH S,ef f t ρ S (t)e +iH S,ef f t is given by the time-local master equatioṅ where ρ nm = E n |ρ|E m and |E m is the m − th eigenstate of H S,ef f with energy E m , which we number in order of ascending energy. This form, containing explicit time-dependence, is often referred to the Non-Secular Bloch Redfield equation, distinguishing it from the Secular Bloch Redfield equations which are obtained from Eq. (7) by only retaining terms in the RHS summation This last approximation is normally justifed when the energy differences between different transition energies E nm − E jk R −1 nnjk , ∀n, m, i, j, so that on the typical timescales on which the density matrix evolves under the action of the Redfield tensor R nmjk , the highly oscillatory terms with E nm − E jk = 0 average to zero. When this is indeed valid, the Secular Bloch Redfield equations have a particularly simple form, as population (diagonal) and coherence (off diagonal) terms of ρ S are completely decoupled. The populations then obey a Pauli (kinetic) master equation, while any coherences present in the intial condition independently and exponentially decay to zero. Crucially, there are no terms in the SBR equations that allow for the ex nihilo generation of coherences.
However, non-secular terms can create dynamical coupling between populations and coherences, as has been widely discussed in the context of ultrafast spectrosopies. The significance of these terms has also recently been highlighted in a number of papers, showing that their inclusion often provides a more accurate description when compared with more advanced numerical treatments of open quantum systems 74 . Interestingly, wellknown problems related to the potential lack of positivity of reduced density matrices under BR evolution have also been shown to arise from problems related to the failure of the Born-Markov asumption, rather than the structure of the master equations when non-secular terms are included 75 .

B. Non-secular generation of spontaneous coherence via incoherent decay
Of particular relevance for our three-level model is the structure of the population-to-coherence elements of the Redfield tensor R nmii = R * mnii , especially the term R D+D−BB which describes the generation of coherence by a population in the initially excited state |B . This is given by Comparing this to the incoherent decay rate from state |B to the lower energy doublet of states (population-topopulation transfer) we can see that, if J(E D±B ), n(E BD± ) and matrix elements V D±B are equal, then R D+D+BB = R D−D−BB = R D+D−BB . For this case, and in the absence of the timedependence associated with the term R D+D−BB in Eq. 7, the incoherent decay of state B would create population and coherence in the lower doublet at an equal rate, i.e would tend to generate a coherent superposition of the lower doublet states. Indeed, if the decay of the |B state is much faster than the oscillation period (set by E D+D− ), one would expect the relaxation to occur into a pure superposition state of |ψ ≈ |D + ± |D − ( with phase depending on the sign of the matrix elements). If the decay is much slower than the oscillatory period, the coherence generation will be negligible, leading to an incoherent mixture of (approximately equal) populations in the doublets. Within the scope of BR theory, the eigenstate level scheme and spectral density that we consider enables the conditions for coherence generation to be almost perfectly met at least for small reorganisation energy when the model parameters remain close to those predicted by the eigenstates of H S . Firstly, the matrix elements for transition between the |B state and the doublets |D ± induced by the spatially local operator V λ (3) are almost equal, due to the real-space delocalisation of the states and the large energy gap between the |B and |D ± manifolds (their fractional difference is no more than ≈ J 23 /J 12 1). Secondly, by applying a spectral function that is symmetric and peaked at an energy 1 2 (E BD+ + E BD− ), we obtain equality of J(E D±B ). Thirdly, by working at a temperature such that βE BD± 1, or βE BD± 1 all relevant elements of the Redfield tensor approach equality. This only leaves the condition that the transition rates for decay of the |B state should be faster than the period of coherent oscillations in |D ± states. This can be controlled by varying the coupling strength to the environment, although we note that obtaining this condition violates the standard application of the Markov approximation, additionally motivating our use of non-perturbative HEOM methods (see below).
Finally, we note that the same qualitative analysis can be applied to coherences generated between the |B and |D ± states due to decay of the |B state. However, the relevant oscillatory time period to compare to the decay rate is now set by the bright-dark state splitting E BD , which is ten times larger than E D+D− . It could therefore be possible to find a parameter space in which large inter-doublet coherence is generated without any significant coherence generation between the bright and dark states. This scenario defines what we mean by coherence generation arising from incoherent relaxation.

C. HEOM
We recall here the derivation already presented in previous works 55,57-62 and give, for a purpose of completeness, the main equations. For the efficiency of the HEOM algorithm, the spectral density is parametrized so that the two-time bath correlation function is expressed as a sum of complex exponential functions 76 : Explicit expressions of the α k and γ k by the analytical integration of Eq.(5) with the superohmic parametrization of the spectral density (6) can be found in the Appendix of ref. 77 . n cor is the sum of the four terms coming from the four simple poles in the upper complex plane and, in principle, an infinite number of terms related to the poles (Matsubara frequencies) of the Bose function on the imaginary axis ν j . In practice, we find that the number of Matsubara terms remains small (about 10) at and above room temperature for the model under study. The complex conjugate of the correlation function can be expressed by keeping the same coefficients γ k in the exponential functions with modified coefficientsα k according to : withα 1 = α * 2 ,α 2 = α * 1 ,α 3 = α * 4 ,α 4 = α * 3 and α j,matsu = α j,matsu where the α m with m = 1, 4 are related to the four poles of the superohmic Lorentzian function and α j,matsu refer to the Matsubara terms 78 .
By assuming an initial factorization of the total density matrix, the time evolution of the reduced density matrix, in interaction representationρ S (t) is given bỹ where ] is the Liouvillian of the system-bath interaction with the system coupling operator in interaction representation S (t) = e iĤ S,ef f tŜ e −iĤ S,ef f t and the bath operator as given above. Expressions (11) and (12) for C (t − τ ) and C * (t − τ ), correspond to the relaxation of n cor effective bath modes. Each set of the corresponding occupation numbers is represented by a collective index n = {n 1 , · · · , n ncor } and is associated to an auxiliary density matrix. The master equation is then written as a time-local hierarchical system of coupled equations among the auxiliary operators. Each matrix can communicate only with the superior and inferior level in the hierarchy for which one occupation number is varied by one unity n ± k = {n 1 , · · · , n k ±, . . . , n ncor } : In this hierarchy of auxiliary density matrices, the system density matrix is given by top row, i.e. for n = {0, · · · , 0} henceρ S (t) = ρ {0,··· ,0} (t). The level of the hierarchy is chosen until convergence is reached for the system density matrix. As previously stated, it can be seen that the equations of motion that determine the reduced density matrix of the system are completely determined by the expansion coefficients of the correlation functions C(t) and C * (t) that appear in Eqs.11 and 12, and which are ultimately determined by the environment spectral function and temperature.
The HEOM equations are efficient to go beyond the second order perturbative regime of the Redfield equations even if the Markovian approximation could still be valid at higher orders. However, stronger couplings are often linked to non-Markovian dynamics. Signature of non-Markovian behaviour for strong system-bath coupling is analyzed in appendix A where we also illustrate convergence of HEOM equations and compare numerically exact HEOM simulations with Redfield results.

A. QUANTUM NOISE
In each simulation, the initial state is the bright eigenstate |B and the bath is at room temperature T = 298K which corresponds to quantum noise for our spectral density, as the peak frequency Ω ≈ 5k B T . We note, however, that this temperature give an energy scale k B T /E D+D− ≈ 3, which would be expected to drive equal (mixed) populations of the doublet population and rapid coherence loss. From the analysis based on the Redfield theory, we first examine the evolution of the main tensor elements to predict the best range of the η parameter (Eq.4) to create the expected long lived superposed state in the doublet. The three main tensor elements of the 'downhill' transitions at T = 298K R D+D+BB , R D−D−BB , R D+D−BB are displayed in Fig. 3 as a function of the η parameter by accounting for the variation of the coupling V λ (Eq.3) and of the eigen energy gap induced by the renormalization energy. The best expected domain for the coherence creation appears to be around η = 0.015 which corresponds here to a renormalisation of about 15 cm −1 . This value provides the largest decay rate while maintaining equality of the three relevant Redfield tensor terms. The reverse rates R BBD−D− , R BBD+D− remain negligible at room temperature due to detailed balance, potentially extending the lifetime of any decay-generated coherent states. Simulation with HEOM will allow us to probe stronger couplings beyond the perturbative regime and examine the stability of the process. Figure 4 shows the population evolution (i.e. the diagonal elements of ρ S (t)) in the |B state (full lines) and in the doublet |D + (dashes) and |D − (dots) for different coupling regimes. In the perturbative regime (Fig. 4a), the decay is monotonous while for the strong coupling  case (Fig. 4b), oscillations occur which are related to quasi-reversible energy exchange between the system and the environment, as the coupling strength exceeds the line width of the spectral function (strong-coupling cavity limit). These latter dynamics lead to features in the measure of non-Markovianity we present in Appendix A, as predicted in ref. 79 . The possible creation of a superposition state is suggested by the simultaneous growing of population in the |D ± doublet states, although this could also arise without generating any coherence in the doublet. Close to the best expected coupling regime for η = 0.01 (red curves), the populations rise monotonically and plateau to equilibrium values expected from the Boltzman distribution at this temperature. The populations of the two doublet states show oscillatory behaviour in the strong regime η = 0.16 (black curves), with a period and duration much longer than the oscillations seen in the decay of the bright state. We will return to this non-perturbative effect in Section IV B. The critical observable, the modulus of the coherence |ρ D+D− (t)| between the doublet states, may be seen in figure 5 as contourplots in a time and η parameter map or in figure 6a for some selected couplings. The amplitude of the created coherence shows a clear dependence on the system-bath coupling. In the optimal situation corresponding to values near η = 0.015, as predicted by the Redfield analysis, the coherence modulus reaches very close to the maximum possible value of 0.5 in about 100 fs and remains stable for more that 1 ps. For very weak coupling, a coherence is observed but its amplitude remains below 0.1. For the strong coupling, the early dynamics leads to a high amplitude around 0.5 but due to the bath interaction the asymptotic value stabilizes below the optimal coherence.
The purity of the system density matrix T r ρ 2 S (t) is given in figure 6b. In the weak coupling cases ( η = 10 −4 or 10 −3 ), the purity is mainly determined by the mixed state with the initial state which is not yet relaxed. The purity confirms that the most favorable situation is the moderate coupling case around η = 0.01, where the asymptotic purity is well above the purity of a Boltzmann mixture at room temperature (≈ 0.5 ) and shows that incoherent relaxation produces a superposition state in the doublet with relatively little entropy generation. Figure 7a illustrates the stability of the created coherence ρ D+D− (t) for a favourable case with η = 0.01 and the difference with the other coherences ρ BD− (t) or ρ BD+ (t). As shown in figure 7b, the latter completely disappears after 500 fs and their amplitudes never exceed 0.015 so they are smaller by more than one order of magnitude. As previously discussed, this establishes the coherence generation arises from an incoherent decay.
These beats could potentially be detected in an experimental set up similar to Ref. 1 . The rapid relaxation effectively prepares a nearly pure superposition state that coherently evolves over a subsequent time t as |ψ(t) ≈  and |D + − |D − = |3 , it can be easily seen that the evolution of the wave function phases leads to oscillatory real-space motion of the excitation between sites 1&2 and 3. This is illustrated in Fig. 8, where a periodic and nearunity population of site three with frequency E D+D− can be seen. As the resonator emission only arises from population of site 3, oscillations in its population should be observable as a periodic modulation in the resonator signal at a frequency E D+D− .
The stability of the coherence creation via the incoherent decay was checked with respect to the shape of the spectral density. We compare the sharp spectral density with a broader one (the parameters are given in the appendix). They are schematized in the inset of Figure 9. The renormalisation energy is calibrated to be nearly equal in both cases. Figure 9a presents the population evolution in the |B and |D ± doublet for η = 0.01. The decay is slower in the broad peak case but the population in the two doublet states is still growing simultaneously. As shown in figure 9b, the coher-  ence ρ D−D+ (t) presents a similar profile in both coupling schemes. We note that this example serves to show that the coherence generation is a result of the coupling matrix elements and transitions rates (thus appearing at the level of the master equation), and do not arise from vibronic mixing effects that require a strong coupling to a resonant and sharp (underdamped) vibrational mode at the Hamitlonian level. This suggests that the conditions for noise engineering, both in quantum simulators and physical realisations, that are required for noise-driven coherence are, in fact, rather lenient. While we have shown in this section that it is possible to find a parameter regime where quantum noise can lead to coherence generation via incoherent relaxation, the longevity of the resulting, near-perfect superposition states is perhaps not so surprising. Due to the absence of incoherent transitions back to the high energy bright state, the only mechanisms of dephasing in the doublet are due to intra-doublet relaxation and/or pure dephasing. For both the peaked and broad spectral densities we have used, the spectral weight at the energy gap E D+D− is extremely small and pure dephasing vanishes at long times 5,8 , making the doublet state effectively decoupled from the environment. As we shall show, this situation changes dramatically in the experimentally relevant case of classical noise.

B. CLASSICAL STOCHASTIC NOISE
At the limit of very high temperature, the bath correlation function becomes real and therefore corresponds to a coloured classically stochastic noise (see Fig. 2). This suggests that it should be possible to simulate the effects of stocastic noise, including any non-markovianity, by making the simulation 'temperature' much larger than the other intrinsic energy scales of the system and bath, while rescaling the coupling to the bath to maintain physically reasonable transition rates. In the Golden Rule approximation, the decay rates depend on J (ω) (n(ω) + 1) which becomes J (ω) k B T /ω at the high temperature limit. In order to keep the transition rates at similar values to those in Section IV A (which also maintains the significance of our parameter η), we set an artifically high temperature (10 3 K) whilst simultaneously dividing the renormalization energy by a factor of k B T . We note that this procedure captures the essential 'infinite temperature' property of classically stochastic noise: the up and down transition rates are now effectively equal so that any coherence created by relaxation is now subject to potentially strong dephasing noise arising from the rapid transitions that drive the system towards the maximally mixed Boltzmann state. Figure 10a presents the population evolution for different η parameters. The case η = 0.01 (red curves) may be compared with the quantum noise case (see Fig. 4) where this coupling range gave optimal generation of the doublet state. With a classical noise case, the super- position is still created on ultrafast timescales but the peak coherence amplitude is smaller (0.3) and coherence is completely destroyed after just 0.1ps, as seen in figure  10b (red curves). This dephasing time is much faster than the period of coherent oscillations in the dark doublet, so no beating can be resolved in the time domain. Reducing the coupling by an order of magnitude leads to large amplitude coherence generation within ∼ 200 fs, but the dephasing rate is now slower than the beating period, allowing about 2 − 3 cycles of beating to be observed over about 1 ps. For very weak coupling η = 10 −4 (blue curves) the population decay is very slow but a coherence of weak amplitude (< 0.1) is created and maintained during the entire decay (about 2ps) of the system to thermal equilibrium. We therefore confirm that the noise-induced generation of coherences survives in the case of classically stochastic noise and, although it is much more fragile, as seen in Fig. 10, there is a parameter regime where it is possible to resolve the quantum oscillations in the temporal domain. Compared to the quantum (cold) case, we note the following differences in coherence generation that may be relevant for future experiments in superconducting quantum cicuits and/or molecular array architectures. Firstly, the near equality of upwards and downward transition rates creates a competition between the fast relaxation needed to generate coherence and the dephasing that arises from the uphill transitions. As shown above, in order to resolve the beats, a compromise must be struck between the maximum possible amplitude of the coherence and the lifetime of the oscillations. The optimum point will depend on the method of detection and the leveraging between acceptable signal-to-noise (favouring large amplitude) and the available time/frequency resolution. Secondly, at larger coupling stengths, oscillations are also seen in the populations of the eigenstates which are damped on the same timescale as the corresponding coherences. Unlike the quantum case, these oscillations are not due to reversible energy exchange with the environment, but result from the coherent real-space motion of the excitation in the doublet states. Again, an approximate but intuitive understanding can be obtained from the structure of the Redfield equations.
Due to the choice of coupling to the environment (local coupling to site 2), the rapid initial non-secular relaxation in the regime of coherence generation (relaxing to a superposition state) can also be seen as a relaxation of the bright state into the non-stationary dark state |D = 1 √ 2 (|D + + |D − ). As all the uphill rates are also the same, subsequent relaxation back to the bright state only arises when this non-stationary state is populated. However, due to the coherent evolution of the |D state, excitations move in real-space to site 3, which is not coupled to the environment, so that the population in the doublet states is temporarily unable to make any uphil transitions. As the oscillatory quantum beats return population to the |D state, the uphill transitions become allowed again, depopulating the dark doublet states in a periodic way. These motion-induced modulations of the uphill transition rate are the origin of the eigenstate population oscillations seen in Fig. 10a.
Interestingly, this novel modulation of the transition rates effectively results in a transient and periodic violation of detailed balance 70,80 , as the suppression of the upwards transitions leads to an 'overshoot' of population transfer from the bright state, as if the bath were (temporally) at a much lower temperture. This effect is most prominent for the intermediate coupling (η = 0.01), where the coherence is both large enough and long-lasting enough to allow a few near-complete oscillations of the excitation between sites 1&2 and 3. This interpretation is confirmed by looking at the populations in the site basis for this coupling, as shown in Fig. 12. Comparing to Fig. 10a, we see that the eigenstate oscillations occur at precisely the times when site 3 is maximally populated and uphill transitions are suppressed. Experimentally, these coherent modulations of the bright state decay would be detectable through the emission in the transmission/excitation waveguide and would be in antiphase with modulations in the resonator signal. Indeed, because these coherent dynamics directly effect the eigenstates populations, it is likely that they will be even easier to detect that the oscillations in the eigenstate coherences that appear in the regime quantum dissipation. We note that this effect is also responsible for the oscillations in the modulus of the coherence that be seen in Fig. 11.
For completeness, we finally comment on the nature of the population oscillations in the quantum case of Fig.  4. These arise from a completely different physics, and the oscillations only occur between the doublet states at strong coupling. This is due to the increasing relevance of the reorganisation energy, compared to the electronic couplings and detunings between the sites of the network. Again, due to the coupling of the bath to just site 2, this has the effect of detuning that site, so that the effective eigenstates in the lower doublet no longer correspond to the fully delocalised states |D + &|D − . The increasing amplitude of the oscillations we observe with stronger coupling are essentially related to the increasing misalignment of the measurement and effective eigenstate basis. It is clear that a non-stationary state of these eigenstates is still prepared by the rapid relaxation of the bright state, but the 'quality' of the super positions that are formed degrades with increasing coupling strength. This degradation is already suggested by the trends in the Redfield rates shown in Fig. 3, and confirmed by the decreasing purity of the system after relaxation in Fig. 6. This strong coupling effect does not occur in the classical case, as the rescaling of the system-bath coupling needed to take the high-temperature limit means that the reorganisation energy is always negligible compared to the system parameters.

V. CONCLUSIONS
In this article we have demonstrated that the system of bright and dark excitons realised in the qubit architecture of Ref. 1 possesses near-ideal properties w.r.t. the non-secular processes that allow incoherent dynamics to generate coherent wave-like motion. By using the nonperturbative HEOM technique, we have verified that under conditions of classical noise, it should be possible to experimentally resolve these coherent dynamics through the related oscillatory motion across the network in realspace, although the effectively high ('infinite') temperature of the classical stochastic bath leads to strong dephasing via incoherent 'uphill' population transitions from the dark manifold. Nevertheless, a previously unanticipated prediction arises from this fact and is seen in our simulations: the periodic violation of detailed balance caused by the suppression of uphill transitions as the excitation moves coherently away from the site-local source of the noise. At the same time, the existence of resolved quantum oscillations also appears to be sensitive to the correlation time of the environment, with nonperturbative theories showing stronger coherent dynamics at the same coupling stength as a simple Redfield approach. The differences persist over the correlation time of the environment, and it may be possible that this is related to the finite 'switch on' time of the rapid transitions that are present from t = 0 in the Markovian theory. Although in general these effects do not correlate with any obvious feature in the measure of non-markovianity, it is clear that 'memory' or temporally non-local effects constitute another handle by which coherent energy transfer might be manipulated. Indeed, the real-time switch from a weaker to stronger dissipative coupling might be more generally important for coherent real-space motion, as suggested for quasi-coherent charge separation in organic bulk heterojunctions 30,[81][82][83] .
In the case of quantum noise, which might be realised by using the multi-level nature of superconducting qubits to simulate a quantum harmonic oscillator 84 , the longevity of the superposition states is only limited by the form of the spectral density and the strength of interactions at the small energy gap between the dark doublet of states, or at zero frequency (pure dephasing). Spectral functions that vanish rapidly at low frequencies while having large amplitudes at the much larger bright-dark energy gap would therefore be advantageous for coherence generation. However, the coherent dynamics are suppressed at very strong coupling -regardless of the shape of the spectral density -by the growing reorganisation energy of the environment, which detunes and localises the low-lying excitations. Finally, our physical understanding of the numerical results has often relied on predictions from Redfield theory that arise from the site-local noise in our model. Given the capabilities of present simulators to apply site-specific noises, it would be very interesting in the future to consider how applying different spatial and spectral correlations to noises impacts coherent dynamics in quantum energy transfer networks.

Appendix A: Non-Markovianity analysis
We briefly discuss the non-Markovianity of the dynamics for some coupling ranges.
Numerous non-Markovianity witnesses have been proposed in the literature 85-87 but we consider here only the volume of accessible states in the generalized Bloch sphere 79 . Eqs. (14) define a time local dynamical map ρ S (t) = φ t [ρ S (0)] which may be expressed in an operator basis set {G m } for the Liouville space of dimension d 2 by a generalization of the Pauli matrices for d = 3. The expansion of the matrix density in this basis leads to the Bloch representation of the system. In matrix form the map reads F m,n (t) = T r (G m φ t [G n ]) and the volume of accessible states in the Bloch sphere can be obtained from the determinant of this matrix V (t) = det(F). (A1) A non-monotonous decrease of this volume is a signature on a non-Markovian back flow from the bath to the system. This is illustrated in figure 13 by the volume of the accessible states in the Bloch sphere for differnt η parameters. For instance, bumps in the volume are obtained for the strong coupling case η = 0.16 (full line). This justifies that dynamics must be treated beyond the Redfield approximation and probably beyond second order regime. The perturbative regime for η < 0.02 leads to a smooth evolution of the population (see Figs. 4 or 10). On the contrary, oscillations are observed during the decay of the bright state for stronger coupling. This behavior may be related to some non Markovian effects characterized by back-flow from the environment to the system. The level of HEOM hierarchy ensuring convergence of the simulation depends on the strength of the systembath coupling. Level L in the hierarchy corresponds to order 2L in perturbation approach. Convergence is checked in figure 14 by analyzing the coherence that always converges more slowly than the populations. The cases with η ≤ 0.02 remain in the perturbative regime. On the contrary, for the case η = 0.16, the regime is obviously non-perturbative and level L = 4 is required.
We do find that although the Markovian 2nd-order perturbative theory (Bloch-Redfield) allows an intuitive, qualitative understanding of the dynamics, it fails to describe the dynamics with quantitative accuracy (c.f. HEOM results). Even at relatively weak coupling, this can lead to significant differences, as shown in Fig. 15. Given the relatively weak coupling involved, we believe that the orgin of these differences are related to the Markov (time local) approximation in Redfield theory, which does not take into account the long correlation time of the environmental spectral functions that we consider. This is supported by the results in Fig. 15, which show that the differences in dynamics become negligible for times longer than the bath correlation time (≈ 200 fs), as shown in Fig. 2. As the oscillatory coherent dynamics induced by the relaxation are much more prominent in the HEOM results, this indicates that a proper treatment of extended bath correlation times can also be an important factor for noise-induced coherence generation in the classical case, further emphasising the need for methods such as HEOM.