From Goldilocks to Twin Peaks: multiple optimal regimes for quantum transport in disordered networks

Understanding energy transport in quantum systems is crucial for an understanding of light-harvesting in nature, and for the creation of new quantum technologies. Open quantum systems theory has been successfully applied to predict the existence of environmental noise-assisted quantum transport (ENAQT) as a widespread phenomenon occurring in biological and artificial systems. That work has been primarily focused on several `canonical' structures, from simple chains, rings and crystals of varying dimensions, to well-studied light-harvesting complexes. Studying those particular systems has produced specific assumptions about ENAQT, including the notion of a single, ideal, range of environmental coupling rates that improve energy transport. In this paper we show that a consistent subset of physically modelled transport networks can have at least two ENAQT peaks in their steady state transport efficiency.


I. INTRODUCTION
Energy transport occurs in many contexts: from circuits and molecular junctions to processes like photosynthesis and the electron transport chain in biology [1][2][3][4][5][6][7][8] . In 2008 the specifics of Environmental Noise-Assisted Quantum Transport (ENAQT) were first laid out in both artificial and natural schema [9][10][11][12] . Since then there has been a proliferation of research into the interplay of coherent quantum transport and noise in many systems [13][14][15] , such as the role environmental phonon coupling can have in overcoming the effects of coherent localisation [16][17][18][19] . From these studies and others, an intuition has been established that ENAQT produces a single 'Goldilocks zone' [20][21][22] where dephasing overcomes limits inherent to fully coherent dynamics, but is not sufficiently aggressive to spoil its favourable transport characteristics, such as those brought about by constructive interference.
Recent work introducing the concept of 'population uniformisation' has made the varying transport behaviour between fully coherent and fully classical limits explicit 37 . Population uniformisation states that the vari- * e.gauger@hw.ac.uk ance in on-site populations has a similar qualitative character to the transport efficiency of open quantum systems, and is minimised when transport efficiency is maximised, even in the presence of disorder or repulsive interactions 18,37 . While this framework explicitly frames things in terms of moving from coherent wavefunctions to classical diffusion and Fick's Law, we yet again see the same, singly peaked ENAQT transport efficiency on the standard systems, including the FMO complex 38 .
In this paper we systematically investigate optimal noise rates across randomly generated transport networks, and show that many have at least two ENAQT peaks or 'Goldilocks Zones' where their transport efficiency is maximised. Our networks are made of twolevel systems, which we model as point dipoles. We arrange these sites with realistic spacing and effective dipole moments to ensure the relevance of our results, and consider both uniform and normally distributed energy landscapes.

A. Network Setup
We set up each system as follows: first we define a sphere with a fixed radius of 10 nm, and place two sites at opposite ends of the sphere to act as injection and extraction sites. We then fill the volume with six additional sites with random positions and dipole orientations to produce a disordered but fully connected network, a similar setup and approach has been utilised in prior works 23,30,39 . Figure 1 shows a typical dipole network generated in this way.
Every site is modelled with an identical dipole moment, we use the effective dipole moments from the bacteriochlorophylls in the FMO complex |d| = 0.114 e · nm( √ 30D) 27,[40][41][42] . With the structure established we can consider the Illustration of a typical dipole system: a fixed volume with injection and extraction sites in at the poles, and randomly oriented and positioned sites within the volume. The system is then coupled to phonon environments with varying strengths (Γ), and the transport efficiency (η) measured. On the right we show transport efficiency curves from a selection of randomly generated networks modelled with phenomenological pure dephasing that exhibit multiple maxima in their transport efficiency.
dipole-dipole interactions, where r is the separation vector between the dipole pair, d i,j are the respective dipole moments and 0 is the vacuum permittivity. To ensure the point dipole approximation remains appropriate we enforce a minimum separation of 1 nm between every pair of sites 43 . Throughout this paper we use nanometres, elementary charge and electron volts, meaning a Coulomb constant (4π 0 ) −1 ≈ 1.43996 eV · nm e −2 .
To generate the dipole on-site energies we use an average on-site energy¯ = 1.5498 eV (12,500 cm −1 ), and will later also sample them from a normal distribution with a standard deviation of 1%, so σ = 0.0155eV, very similar to the energies and disorder found in the FMO complex 40 . The default values we use are summarised in Table I.  With the on-site energies and dipole-dipole interactions defined we can construct the Hamiltonian. We assume there is only a single excitation in the system at any time and construct the following excitonic Hamiltonian where i and j are site indices, and the V i,j are evaluated according to equation (1). We use the single excitation approximation for computational efficiency as it gives an N × N Hamiltonian, and is widely used when modelling open quantum systems with low injection rates, such as light harvesting complexes 22,23,30,39,[44][45][46] .

B. Dynamics
To look at the effect of the environment and finite temperatures on transport properties, we use the full nonsecular Bloch-Redfield master equation 47,48 where ρ s is the system Hamiltonian and ω are the eigenenergy splittings 47,48 . The system-environment interactions, A m , are derived by transforming the relevant site basis operators A deph,i = 2 |i i| − I into the Hamiltonian eigenbasis and S mn (ω) defines the noise-power spectrum associated with the system-environment interaction 30,47,48 . The noise-power spectrum function is where N BE (ω) defines Bose-Einstein statistics at a given phonon inverse temperature β, Θ(ω) is the Heaviside function, allowing phonon-assisted transitions from higher to lower eigenenergies (ω > 0), and J (ω) is the phonon spectral density 48 . In this work we use the Drude-Lorentz spectral density, which has previously been used to model excitonic transfer in light harvesting complexes 49,50 , where Γ scales the rate of noise in the system from interactions with the environment, τ is the correlation time, and τ −1 is the spectral density peak frequency. Finally, we have D [A] ρ which is the dissipator superoperator To model extraction and injection, a shelf state is appended to the system. The extraction operator A ext projects population from the extraction site to the shelf state, A ext = σ + shelf σ − ext , and then that population is reinjected from the shelf state back onto the injection site with the injection operator A inj = σ + inj σ − shelf . Injection and extraction are matched, γ ext,inj = 0.1 eV, changing this value generally changes quantitative values but not the qualitative behaviour 37 .
To complement the Redfield calculations, we also carry out phenomenological pure dephasing calculations with the Lindblad master equation 47,51,52 where A deph,i are Lindblad operators describing the environmental influence on each site i, again these are of the form A deph,i = 2 |i i| − I 48 . All other symbols have the same meaning as in equation (3). This approach is equivalent to the nonsecular Bloch-Redfield master equation for an infinite temperature and a flat spectral density 25 .
We focus here on the steady state ρ ss which is found by calculating the null vector of the system evolution Liouvillian. Our figure of merit then is the excited steady state population on the extraction site This is motivated by the strong correspondence found in prior work between dynamical and steady state transport properties 18,53 , as well as further studies suggesting that the steady state is more natural for photosynthetic systems 2,54,55 . The steady state approach also allows us to avoid any confusion that could arise from the influence of transient effects when comparing different networks.
For each network considered in this paper, we were interested in how this transport efficiency η changes with Γ, the noise rate from coupling to the environment. To do this we considered a large range of noise rates Γ → 10 −7 − 10 eV , and for each value recorded η as well as the full steady state population. This range of Γ was chosen as it was broad enough to capture the values where η has maxima for our networks, and additionally show the transport efficiency decreasing outside these peaks as shown in figure 1.
These results were then filtered to ensure validity, all data presented here has passed checks on the unity of the steady state trace, non-negativity of on-site populations and steady state eigenvalues (see appendix K). From that point we could perform simple peak-finding calculations for each network and spectral density to directly identify in which cases there was more than one optimal noise rate or peak in the transport efficiency curves, and how often this occurred.

A. Environmental Effects
We generated several ensembles of 1,000 dipole networks, each comprising two fixed and six randomly located sites. For our Redfield calculations the key spectral density parameters are the temperature, and the peak frequency of the Drude-Lorentz spectral density (τ −1 ). Lindbladian pure dephasing acts as our infinite temperature limit, being equally present at all frequencies.
We start by considering an ensemble of networks with identical splitting between the two levels on each site, and see how many networks have multiple maxima in their transport efficiency. This approach lets us compare our results to prior works that have made the same assumption of uniform on-site energies when modelling disordered molecular networks and other complexes with dipole interactions 24 We see only a slight sensitivity to the peak frequency in this case, and a slight increase in the behaviour at low temperatures. This is likely because the eigenenergy splittings are quite small when the on-site energies are identical, and as such the transport behaviour depends on the near zero behaviour of the Drude-Lorentz spectral density, which is a relatively stable region. Figure 2 shows a surprising result, contrary to prior wisdom we consistently find about 6% of networks have multiple peaks in their transport efficiency, regardless of the spectral density or temperature. This illustrates that these fully connected networks can have multiple maxima in their transport efficiency, but neglects the importance of on-site energies in transport 25 . We consider the effects of varied on-site energies in figure 3, where the energies are normally distributed as described in  show our main results from this paper. We see that in every situation we simulate, a sizeable subset of our networks have two peaks in the their transport efficiency. The Drude-Lorentz peak frequency has a slight effect on how often we observe this behaviour, but there is a more pronounced sensitivity to temperature. The lower the temperature, the more often we see this behaviour. By extension, this is seen least often -but still clearly represented -in the Lindblad pure dephasing limit. We show what proportion of these results occur within measured FMO reorganisation energies in appendix G.
We note that this double peaked phenomenon generally occurs more frequently in the energetically uniform ensemble. We attribute this to energetic disorder producing greater localisation. Meaning that not only are the energetic differences between eigenstates larger, but also those eigenstates are more tightly confined to specific sites. This greater spatial confinement means there are fewer pathways from injection site to extraction site. The greater energetic differences also raise the chance of some eigenstates being so far detuned that they are effectively inaccessible given the finite range of noise rates Γ we consider.

B. Structural Effects
As shown on the right of figure 1, different networks have transport efficiency curves that take a variety of forms. The relative prominence and position of each transport efficiency peak varies. In some cases peaks are very well separated with a pronounced dip between them. In other cases the two are barely distinguishable, almost smearing into each other. This suggests a strong structural dependence between the form of the curves and the Hamiltonians of their respective networks.
We closely inspected networks such as the one in Figure 4, which produce double peaks across a wide range of temperatures and Drude-Lorentz peak frequencies, to identify what key features might be correlated with having multiple peaks in the transport efficiency of a system. See appendix F for the physical properties of the network. We did not find any strong geometric dependence across networks with multiple maxima in their transport efficiency, but we identified consistent features in the system eigenstates. Specifically, these systems often have one or more large gaps in their eigenenergy distributions, as shown in figure 4 (c).

C. Discussion
An intuitive and tempting interpretation of the multiple optimal transport regimes would be that two broad sets of energy scales are important: the width of each band of eigenenergies, and the separation between bands. The intuition from this being that there is an initial path with an ideal environmental coupling for ENAQT. However, as phonon couplings become larger, eigenstates that were previously far-detuned become more accessible. If the new path is more efficient for exciton transport, then a new peak in the transport efficiency can be observed. Both the energetically uniform and non-uniform ensembles consistently have these gaps in their spectra, induced by the dipole-dipole interactions between sites, and if present, differences in on-site energies.
Testing this hypothesis in appendix A we do indeed find a positive correlation with the relative standard deviation of eigenenergy splittings. However the change in the fraction of multiply peaked networks remains modest, suggesting there are other factors at play. In the following we briefly summarise how doubly peaked behaviour correlates with some other aspects.
In appendix B we tested the energetic separation hypothesis in another way, generating a new independent ensemble of 1000 networks where three of the dipoles have a fixed offset added to their on-site energies. This approach encourages more gaps to form in the eigenspectra and we see a modest increase in the frequency of double peaks.
We also have considered the relative energies of the injection and extraction sites in these systems. In appendix C we show that double peaks occur more frequently when injecting at lower energies than the extraction site, suggesting it may occur more often in less efficient networks (in the sense of requiring 'uphill' energy transport), albeit by no means limited to those. In appendix H we directly compare the maximum transport efficiency of single-peaked and double-peaked networks, and show that double-peaked networks have a larger spread in their transport efficiencies, but can be just as efficient as the single-peaked networks. Appendix I is then concerned with how relevant each peak is in doublepeaked systems. We find that that both peaks typically have a similar prominence, though the peak at higher system-environment couplings tends to be more efficient.
Another consideration is the number of potential paths in a system from the injection site to the extraction site. We generated an ensemble of networks made of paired, disordered, nearest-neighbour chains that only connected at shared injection and extraction sites, giving only two paths across the system. In appendix D we show that while this strongly reduces how often a network has multiple transport efficiency maxima, we do still observe it against all spectral densities. We further show in appendix E that double peaks can be observed against Ohmic and superohmic spectral densities as well.
To consider the effect of system density, in appendix J we reduce the minimum separation between sites and the total system volume to better match the chromophoric density seen in light-harvesting complexes. Again we find a similar subset of networks with multiple optimal noise rates, though one that less favours multiple ENAQT peaks at low temperatures.
Overall, our analysis suggests there are a multitude of factors at play which can positively correlate with an increased occurrence of doubly peaked networks. The analysis in this paper has been focused on networks with double peaks as that is what we observe for these sys-tems. We believe that more than two ENAQT peaks are possible, and that networks with more sites and potential paths from source to sink may present such behaviour. Given the large amount of parameters involved in these systems we have presented many conditions that allow for these multiple peaks to occur in a range of system geometries, but do not find any condition that strongly correlates with the multiple ENAQT peaks being present.

IV. CONCLUSION
In this paper, we have shown that transport networks with realistic and microscopically resolved vibrational interactions frequently feature more than one optimal regime for transport efficiency. This runs counter to expectations of there being a single 'Goldilocks zone' in any and all cases. We observe that these multiple optimal transport regimes can occur for energetically ordered or energetically disordered networks, and occur more often in networks with less evenly spaced eigenvalues.

CONFLICTS OF INTEREST
There are no conflicts to declare. positively correlated with an increased fraction of networks displaying two optimal transport regimes. The probability density histograms against the relative standard deviation are shown for the energetically disordered and energetically uniform ensembles in figure  Histogram of probability density against the relative standard deviation of network eigenenergy differences for 1,000 energetically uniform networks, previously described in figure 2. Without an energy landscape, we see a wide spread of relative standard deviations defined by the geometric properties of the networks. We see a relatively sharp maximum relative standard deviation here due to the exclusion volume or minimum distance we enforce between dipoles when generating our systems.

Appendix B: Networks with offset energies
Based on our intuition of double peaks being related to two effective pathways with different optimal conditions for ENAQT, we generate a new, independent ensemble of networks with two different on-site energy scales. We generate the dipole networks exactly as before, including their random on-site energies. Then for each network we randomly pick three of the bulk dipoles and shift their energy upwards by the standard deviation of our on-site energies (+15.5 meV). This adjustment to the networks increases the probability of there being a larger gap in the system eigenenergies, but leaves the geometric properties of these networks unaffected. As shown in figure 7, we do see some increase in double-peaked transport efficiency in most cases. Though as expected when comparing two independent datasets there are fluctuations in the trends.

Appendix C: Injection and extraction energies
We also consider the importance of where the injection and extraction sites sit in the eigenenergy landscape. We index the eigenstates from lowest to highest energy (λ 1−→N ), and record which eigenstate is most present on the injection site (λ inj ) and the extraction site (λ ext ), respectively.
For our networks with no double peaks, there is general symmetry. For the networks with multiple transport efficiency maxima, we see a preference for injecting at lower eigenenergies than they extract at (λ inj − λ ext < 0). This suggests the effect will be more prominent in less efficient systems. Though we note that the effect is present broadly, also clearly occurring in networks where energy transport should be efficient along a downhill gradient. Figure 8 shows these results for energetically disordered networks, and figure 9 shows the same for energetically uniform networks. We see the same trend in both results, however for the energetically uniform case the eigenstates are more delocalised, often being spread over two or more sites. As a result there are multiple cases where λ inj − λ ext = 0 due to the sites sharing a pair of eigenstates which are equally present on the injection and extraction sites. This does not occur in the energetically disordered ensemble because of the additional localisation.

Appendix D: Double peaks in coupled chains
To further examine the importance of different energy scales, we constructed a simple system from 8 sites. We define two nearest neighbour chains of three sites, and connect them to an injection site and extraction site at either end. Both chains have the same NN coupling (2.5 meV) and average on-site energy (1.55 eV), but have different amounts of disorder in on-site energies. We set their standard deviations to σ = 15.5 and 1.55 meV, respectively, such that each arm will typically have a different amount of phonon coupling that is optimal for overcoming localisation.
This gives us a scenario with two well-separated paths between injection and extraction, rather than the many possible traverses in our dipole networks. As such, there is a general reduction in secondary pathways that have an opportunity to improve transport efficiency. The results of these calculations for 1,000 networks are shown in figure 10 10. Double-peak rates for 1000 'two-armed' networks. We see a decrease in double-peak behaviour in every circumstance compared to our main results shown in figure 3. This suggests that the reduction of available paths or long-range coupling is limiting how often double peaks occur.
As figure 10 shows, double peaks are still present, though always to a lesser degree than in our totally random dipole ensembles. A key difference is the large decrease of double peaks with Lindblad pure dephasing, or at high temperatures but low peak frequencies. So just having two possible pathways across a system is not enough to remove the possibility of double peaks, but does lower the chances of observing it.
We observe that when secondary peaks appear in these scenarios, they are often at very high values of Γ, compared to the range we use to see the same behaviour in pure dephasing and Drude-Lorentz models. This is slightly mitigated in the energetically uniform networks where the lack of disorder in on-site energies has the effect of moving these peaks to lower coupling strengths where our standard approach can capture them. We present a clear example in figure 11 of a single energetically uniform network showing double peaks at all temperatures tested for the Ohmic and superohmic distributions.  As such we can state that these double-peaked effects can and do also occur for these power law spectral densities. However, they occur over a much broader range of environmental couplings, and as such, alternative methods suited to intermediate-and strongly-coupled open quantum systems would be needed to provide more robust statistics.

Appendix F: Hamiltonian Details
Here we list out the data for the specific dipole network presented in figure 4.

Appendix G: Reorganisation Energies
The goal of this work has been to demonstrate the existence of a new phenomenon in otherwise very typical transport networks. But answering this question is separate from demonstrating if this occurs at typical amounts of coupling to the environment. Also, answering this question cannot be done for the pure dephasing approach used in many studies, and instead requires microscopic couplings to be considered. As such we consider our Bloch-Redfield results in this section, and we calculate the reorganisation energies of our networks in the standard fashion where λ is the reorganisation energy and the other terms keep their definitions from equation (4). With this, we can constrain our results to only those occurring at physically reasonable reorganisation energies, which we take to be anything below 36 meV (290.4 cm −1 ), covering previously calculated reorganisation energies for FMO 60 . Figure 12 shows what percentage of networks with multiple transport efficiency maxima have both peaks below this upper bound. Here we consider the maximal efficiency η max of the networks in our main data with either a single peak or two peaks. As shown in figure 13 we see a very large overlap in the efficiencies shown by either kind of transport efficiency landscape. Suggesting that doublepeaked systems can be less efficient, but have such a spread of efficiencies that they are often equally as efficient as the singly peaked systems.

Appendix I: Relative efficiency of each peak
We now study the relative efficiency of both peaks in the double peaked case. We define the steady state transport efficiency of the low noise rate and high noise rate peak as η low and η high respectively. In figure 14 we consider the ratio of these two quantities η low η high for our main ensemble of energeticaly disordered networks. As figure 14 shows, the peaks at higher noise rates are typically more efficient than those at lower noise rates. However we also note that for the vast majority of systems the two peaks have efficiencies less than a factor 2 apart. The central two bars of the histogram correspond to the range 0.5 < η low η high < 2 and make up 69.1% of the doubly peaked systems. So in most cases, both peaks have a similar prominence.

Appendix J: Dense Networks
In this work we have considered relatively sparse networks in a volume with a 10 nm radius. This meant there was lots of freedom for the dipoles and their respective exclusion volumes to form many different kinds of structures, while respecting the ideal dipole approximation that we used. However if one wants to understand light-harvesting complexes, these typically occur at higher densities.
Here we briefly consider such a dense system, keeping the same model as before with 8 sites, but reduce the full sphere radius to 2.5 nm, close to prior work 23 and similarly reduce the exclusion volume around each site to 0.5 nm so that all 6 interior sites can still fit in the volume. The results are shown in figure 15, where we see comparable results to those in sparse networks, but with a decrease in the incidence of double peaks at low temperatures. The reorganisation energies of these peaks were also considered in figure 16, where we see that the results at 30K are consistently the least likely to occur below the cutoff value.

Appendix K: Analysing Data
In this work we have used both the Lindblad and non-secular Redfield master equations. By its mathematical construction the Lindblad master equation always produces completely positive density matrices which preserve the trace. Redfield master equations are also linearly trace preserving but can produce unphysical density matrices with negative probabilities 48,61 . As this work is looking for peaks in arrays of values, we need to screen these erroneous points as sudden spikes or dips on the otherwise smoothly varying data produce false peaks.
To ensure the steady state populations were physical we checked if the trace was unitary and if all the on-site populations were positive. To ensure the steady states were valid we recorded the eigenvalues of each Redfield tensor steady state and then checked their eigenvalues were all between 0 and 1. Tolerances of 10 −5 were used for the site checks, and 10 −4 for the eigenvalue checks as these were sufficient to remove erroneous points. The points excluded occurred at higher system-environment couplings, while results at lower couplings were rarely if ever excluded. These checks were also applied to the Lindblad results for consistency. With these points removed, simple peak finding algorithms were applied to the remaining valid points in each array: no requirements were placed on the peak prominence, heights or widths.