Scattering Theory of Efficient Quantum Transport across Finite Networks

TL;DR

This paper develops a scattering theory demonstrating how designed disorder in finite networks can significantly enhance quantum excitation transfer efficiency, especially under centrosymmetry and spectral doublet conditions.

Contribution

It introduces a novel scattering framework showing how disorder and spectral properties optimize quantum transport in finite networks.

Findings

Disorder can accelerate excitation transfer compared to simple structures.

Centrosymmetry and spectral doublets are key to efficiency enhancement.

Fluctuations in doublet splitting and coupling strength control transfer efficiency.

Abstract

We present a scattering theory for the efficient transmission of an excitation across a finite network with designed disorder. We show that the presence of randomly positioned networks sites allows to significantly accelerate the excitation transfer processes as compared to a dimer structure, if only the disordered Hamiltonians are constrained to be centrosymmetric, and to exhibit a dominant doublet in their spectrum. We identify the cause of this efficiency enhancement in the constructive interplay between disorder-induced fluctuations of the dominant doublet's splitting and the coupling strength between the input and output sites to the scattering channels. We find that the characteristic strength of these fluctuations together with the channel coupling fully control the transfer efficiency.

Figures

Equations

$$S(E)=1-2\pi i{\hat{W}}^{†}\frac{1}{E-H_{\mathrm{eff}}}\hat{W},$$

$$\begin{array}{rl}S(E)& \equiv 1-2i{W}^{†}\frac{1}{E-H_{\mathrm{eff}}}W,\\ & \text{with}{H}_{\mathrm{eff}}=H-iW{W}^{†},\end{array}$$

$$p_{c\rightarrow c^{\prime}}(E)\equiv\mathopen{}\mathclose{{}\left\lvert S_{c,c^{\prime}}(E)}\right\rvert^{2}\,,$$

$$\tau_{c\rightarrow c^{\prime}}(E)\equiv{\rm Im}\mathopen{}\mathclose{{}\left\{{S_{c,c^{\prime}}(E)}^{-1}\frac{\rm d}{{\rm d}E}S_{c,c^{\prime}}(E)}\right\}\,.$$

$$H=\left(\begin{array}{ccc}{E}^{\prime }& v_1\dots v_n& V\\ v_1& & v_n\\ ⋮& H_{\mathrm{int}}& ⋮\\ v_n& & v_1\\ V& v_n\dots v_1& E^{\prime }\end{array}\right),$$

$$W=\sqrt{\frac{\Gamma}{2}}\mathopen{}\mathclose{{}\left|{\rm in}}\right\rangle\mathopen{}\mathclose{{}\left\langle\Psi_{\rm in}}\right|+\sqrt{\frac{\Gamma^{\prime}}{2}}\mathopen{}\mathclose{{}\left|{\rm out}}\right\rangle\mathopen{}\mathclose{{}\left\langle\Psi_{\rm out}}\right|\,,$$

$$H_{\rm eff}=H-i\frac{\Gamma}{2}\mathopen{}\mathclose{{}\left|\rm in}\right\rangle\mathopen{}\mathclose{{}\left\langle\rm in}\right|-i\frac{\Gamma^{\prime}}{2}\mathopen{}\mathclose{{}\left|\rm out}\right\rangle\mathopen{}\mathclose{{}\left\langle\rm out}\right|.$$

$$\mathopen{}\mathclose{{}\left|\pm}\right\rangle\equiv\frac{1}{\sqrt{2}}(\mathopen{}\mathclose{{}\left|\rm in}\right\rangle\pm\mathopen{}\mathclose{{}\left|\rm out}\right\rangle)$$

$$H_{\rm eff}=\begin{pmatrix}E^{\prime}+V-i\frac{\Gamma}{2}&\mathopen{}\mathclose{{}\left\langle\mathcal{V}^{+}}\right|&&\\ \mathopen{}\mathclose{{}\left|\mathcal{V}^{+}}\right\rangle&H^{+}_{sub}&&\\ &&E^{\prime}-V-i\frac{\Gamma}{2}&\mathopen{}\mathclose{{}\left\langle\mathcal{V}^{-}}\right|\\ &&\mathopen{}\mathclose{{}\left|\mathcal{V}^{-}}\right\rangle&H_{sub}^{-}\end{pmatrix}.$$

$$\mathopen{}\mathclose{{}\left\lvert\mathopen{}\mathclose{{}\left\langle\eta^{\pm},\pm}\right\rangle}\right\rvert^{2}=1-\epsilon\,,$$

$$\begin{split}&{\cal E}^{\pm}\approx E^{\prime}\pm V-i\frac{\Gamma}{2}+s^{\pm},\\ &\quad{\rm where}\quad s^{\pm}=\sum_{i}\frac{\mathopen{}\mathclose{{}\left\lvert\mathopen{}\mathclose{{}\left\langle\mathcal{V}^{\pm},\psi^{\pm}_{i}}\right\rangle}\right\rvert^{2}}{E^{\prime}\pm V-e^{\pm}_{i}},\end{split}$$

$$S_{\mathrm{in}\to \mathrm{out}}(E)$$

$$\begin{split}p(E)\approx&\frac{\frac{\Gamma^{2}}{4}(2V+\Delta s)^{2}}{\mathopen{}\mathclose{{}\left((E^{\prime}-V+s^{-}-E)^{2}+\frac{\Gamma^{2}}{4}}\right)}\\ &\qquad\times\frac{1}{\mathopen{}\mathclose{{}\left((E^{\prime}+V+s^{+}-E)^{2}+\frac{\Gamma^{2}}{4}}\right)},\end{split}$$

$$\begin{split}\tau(E)\approx&\Bigg{\{}\frac{\Big{(}\Gamma^{2}+4{E^{\prime}}^{2}+4E^{\prime}\mathopen{}\mathclose{{}\left(s^{-}+s^{+}-2E}\right)\Big{)}}{\mathopen{}\mathclose{{}\left(\Gamma^{2}+4\mathopen{}\mathclose{{}\left(E^{\prime}+s^{-}-V-E}\right)^{2}}\right)}\\ &\qquad\times\frac{4\Gamma}{\mathopen{}\mathclose{{}\left(\Gamma^{2}+4\mathopen{}\mathclose{{}\left(E^{\prime}+s^{+}+V-E}\right)^{2}}\right)}\Bigg{\}}\\ &+\Bigg{\{}\Bigg{(}\frac{-2s^{-}(V+E)+{s^{-}}^{2}+{s^{+}}^{2}}{\mathopen{}\mathclose{{}\left(\Gamma^{2}+4\mathopen{}\mathclose{{}\left(E^{\prime}+s^{-}-V-E}\right)^{2}}\right)}\\ &\qquad+\frac{2\mathopen{}\mathclose{{}\left(V\mathopen{}\mathclose{{}\left(s^{+}+V}\right)-s^{+}E+E^{2}}\right)}{\mathopen{}\mathclose{{}\left(\Gamma^{2}+4\mathopen{}\mathclose{{}\left(E^{\prime}+s^{-}-V-E}\right)^{2}}\right)}\Bigg{)}\\ &\qquad\times\frac{8\Gamma}{\mathopen{}\mathclose{{}\left(\Gamma^{2}+4\mathopen{}\mathclose{{}\left(E^{\prime}+s^{+}+V-E}\right)^{2}}\right)}\Bigg{\}}.\end{split}$$

$$\displaystyle E=E^{\prime}+\overline{s}\,,{\rm for}\ \Gamma\geq\mathopen{}\mathclose{{}\left\lvert 2V+\Delta s}\right\rvert\,,$$

$$\displaystyle E=E^{\prime}+\overline{s}\pm\frac{1}{2}\sqrt{(2V+\Delta s)^{2}-\Gamma^{2}}\,,{\rm for}\ \Gamma<\mathopen{}\mathclose{{}\left\lvert 2V+\Delta s}\right\rvert\,,$$

$$\begin{split}&p(E^{\prime}+\overline{s})\approx\frac{\Gamma^{2}}{V^{2}}\frac{(2V+\Delta s)^{2}}{\mathopen{}\mathclose{{}\left(\Gamma^{2}+(2V+\Delta s)^{2}}\right)^{2}},\\ &p\mathopen{}\mathclose{{}\left(E^{\prime}+\overline{s}\pm\frac{1}{2}\sqrt{(2V+\Delta s)^{2}-\Gamma^{2}}}\right)\approx 1,\end{split}$$

$$\begin{split}&\tau(E^{\prime}+\overline{s})\approx\frac{4\Gamma}{\Gamma^{2}+(2V+\Delta s)^{2}},\\ &\tau\mathopen{}\mathclose{{}\left(E^{\prime}+\overline{s}\pm\frac{1}{2}\sqrt{(2V+\Delta s)^{2}-\Gamma^{2}}}\right)\approx\frac{2}{\Gamma}.\end{split}$$

$$\tilde{\Gamma }\equiv \frac{\Gamma }{2V},\Delta \tilde{s}\equiv \frac{\Delta s}{2V}.$$

$$P(\Delta \tilde{s})=\frac{1}{\pi }\frac{\tilde{\sigma }}{{\tilde{\sigma }}^2+(\Delta \tilde{s}-{\tilde{s}}_0{)}^2},$$

$$\tilde{\Gamma}<\mathopen{}\mathclose{{}\left\lvert 1+\Delta\tilde{s}}\right\rvert\,,$$

$$\mathrm{Prob}$$

$$=\mathrm{Prob}(\Delta \tilde{s}>\tilde{\Gamma }-1)+\mathrm{Prob}(\Delta \tilde{s}<-1-\tilde{\Gamma })$$

$$=1-{\int }_{-\tilde{\Gamma }-1}^{\tilde{\Gamma }-1}\mathrm{d}\Delta \tilde{s}P(\Delta \tilde{s})$$

$$\displaystyle=1-\frac{1}{\pi}\arctan\mathopen{}\mathclose{{}\left(\frac{\tilde{\Gamma}-1-\tilde{s}_{0}}{\tilde{\sigma}}}\right)$$

$$\displaystyle\qquad\quad-\frac{1}{\pi}\arctan\mathopen{}\mathclose{{}\left(\frac{\tilde{\Gamma}+1+\tilde{s}_{0}}{\tilde{\sigma}}}\right).$$

$$p(E^{\prime }\pm V+s^{\pm })\approx \frac{(1+\Delta \tilde{s}{)}^2}{(1+\Delta \tilde{s}{)}^2+{\tilde{\Gamma }}^2/4},$$

$$P(p)=\int_{\mathbb{R}}{\rm d}\Delta\tilde{s}\,P(\Delta\tilde{s})\delta\mathopen{}\mathclose{{}\left(p-\frac{(1+\Delta\tilde{s})^{2}}{(1+\Delta\tilde{s})^{2}+\tilde{\Gamma}^{2}/4}}\right)\,,$$

$$\delta(f(x))=\sum^{m}_{j=1}\frac{\delta(x-x_{j})}{\mathopen{}\mathclose{{}\left\lvert f^{\prime}(x_{j})}\right\rvert}$$

$$\begin{split}P(p)=&\frac{\tilde{\Gamma}}{4\pi\sqrt{p(1-p)^{3}}}\\ &\times\Bigg{(}\frac{\tilde{\sigma}}{\tilde{\sigma}^{2}+\mathopen{}\mathclose{{}\left(1+\tilde{s}_{0}-\frac{\tilde{\Gamma}}{2}\sqrt{\frac{p}{1-p}}}\right)^{2}}\\ &\qquad+\frac{\tilde{\sigma}}{\tilde{\sigma}^{2}+\mathopen{}\mathclose{{}\left(1+\tilde{s}_{0}+\frac{\tilde{\Gamma}}{2}\sqrt{\frac{p}{1-p}}}\right)^{2}}\Bigg{)}.\end{split}$$

$$\displaystyle(H_{\rm int})_{i,j}\sim{\rm Normal}\mathopen{}\mathclose{{}\left(0,(1+\delta_{i,j})\frac{\xi^{2}}{N}}\right),$$

$$\displaystyle v_{i}\sim{\rm Normal}\mathopen{}\mathclose{{}\left(0,\frac{\chi^{2}}{N}}\right),$$

$$\displaystyle E^{\prime}\sim{\rm Normal}\mathopen{}\mathclose{{}\left(0,2\frac{\xi^{2}}{N}}\right)\,,$$

$$\tilde{\sigma }=\frac{{\chi }^2}{V\xi },\text{and}{\tilde{s}}_0=\frac{{\chi }^2}{2{\xi }^2}.$$

Full text

I Introduction

Excitation transport across finite, discrete and disordered networks Scholak et al. (2010); Levi et al. (2015); Walschaers et al. (2016) defines an abstract model for a variety of quantum transport problems, with applications to realistic physical scenarios which reach from natural or artificial light harvesting Spallek et al. (2017) to the physics of cold Rydberg gases Gurian et al. (2012); Scholak et al. (2014); Deng et al. (2016). Beyond fundamental aspects of disorder-induced localisation on discrete networks, this setting also defines an interesting incidence of quantum control in the presence of (and, possibly, through) disorder, which, by the very nature of disordered systems, enforces a statistical approach. Furthermore, when it comes to the specific context of light harvesting, one faces finite networks embedded into hierarchical super-structures Dostál et al. (2016); Kramer and Rodriguez (2017), with hitherto only barely understood interfacing. For a conceptual understanding of the necessary structural elements which ensure the functionality of light harvesting units, it is indispensable to elucidate how the associated, broadly distributed length, energy and time scales are orchestrated, and how the different elementary building blocks are interconnected Ringsmuth et al. (2012).

As a first elementary ingredient of a properly equipped toolbox for the modular modelling of such hierarchical structures, our present contribution establishes a scattering theoretical description of statistical control of point-to-point quantum transport across disordered, finite networks. We adopt the perspective that an excitation, collected e.g. by the antenna complex of a photosynthetic functional unit, is injected into the network at a specific input site, and extracted at an output site from where it is channeled towards the reaction centre – the sub-unit where the incoming photon’s energy is used to drive the ATP cycle Dostál et al. (2016). We elucidate the statistics of the resonance structures in the associated transmission cross sections of the network, together with the concomitant excitation transfer times. In particular, we investigate the interplay between network structure and effective coupling to input and output leads, under the specific, coarse grained constraints of centrosymmetry and the presence of dominant doublet states in the networks’ spectra.

II Model

We build our model as a network described by a Hilbert space $\mathbb{C}^{N}$ (given $N$ network sites). To interconnect the network with the structures it is embedded in, it is attached to scattering channels (or leads) which are described by a space ${\cal H}_{c}\subset\mathcal{L}^{2}(\mathbb{R}^{3})$ (because we generally consider wave propagation in three spatial dimensions). We can thus describe the full Hilbert space of the network together with the external channels as ${\cal H}_{\rm total}=\mathbb{C}^{N}\bigoplus_{c}{\cal H}_{c}.$ Physically this means that we consider a set of bound states on the network, which are coupled to a continuum – a setting familiar from nuclear and atomic physics Fano (1961); Feshbach (1958, 1962, 1967). An incoming excitation is described as a wave packet which is coupled into the network via one of the leads, and subsequently creates an excitation on the network. This excitation will subsequently decay, with variable delay, into either one of the leads.

II.1 Transfer Probability and Dwell Times

For a rigorous description of this scenario we introduce the scattering matrix as defined111Direct application of the formalism of Feshbach (1958, 1962, 1967) leads to

$S(E)=\mathbb{1}-2\pi i\hat{W}^{{\dagger}}\frac{1}{E-H_{\rm eff}}\hat{W}\,,$

hence $W=\sqrt{\pi}\hat{W}.$ Both conventions appear in the literature, and we choose (1) for notational convenience. Brouwer et al. (1997); Celardo and Kaplan (2009); Haake et al. (1992); Lewenkopf and Weidenmüller (1991); Šeba et al. (1996); Stöckmann et al. (2002) by:

[TABLE]

where $H$ is the network’s Hamiltonian represented by an $N\times N$ matrix. The latter encodes the relative positions of the network’s nodes (e.g., a set of chlorophyll molecules) and their mutual couplings. The scattering matrix depends on the energy of the incoming particle, $E$, and has a dimension given by the number $n_{c}$ of attached channels. Also the operator $W$ is described by a matrix, but of dimension $n_{c}\times N$, with entries which determine the coupling between leads and network sites. Consequently, $H_{\rm eff}$ is a non-Hermitian Rotter (2009) $N\times N$ matrix, while $S(E)$ has dimension $n_{c}\times n_{c}$.

Note that every scattering channel does itself support a continuum of modes with continuously distributed energies. Thus, there may be an additional energy dependence of $W$, representing how different modes of the same channel couple to the network. However, we here assume that such energy dependence can be omitted – a common approximation both in light-matter interactions Cohen-Tannoudji et al. (1998) and in mesoscopic physics Lewenkopf and Weidenmüller (1991). Furthermore, the present formulation (1) of the scattering matrix assumes the Lamb shift Cohen-Tannoudji et al. (1998) to be negligible.

We now quantify our control target, the excitation transfer efficiency across the network, in the present scattering theoretical context. We define two figures of merit: the transfer probability $p_{c\rightarrow c^{\prime}}(E)$, and the dwell time $\tau_{c\rightarrow c^{\prime}}(E)$ – from channel $c$ to channel $c^{\prime}$ –, which both depend on the injection energy $E$. The transfer probability reads Brouwer et al. (1997)

[TABLE]

and the dwell time Berkolaiko and Kuipers (2010); Brouwer et al. (1997); Kuipers and Sieber (2008); Šeba et al. (1996); Smith (1960) is given by Brouwer et al. (1997); Smith (1960)

[TABLE]

The latter quantity denotes the phase shift imprinted on an incoming plane wave during the scattering process, and can be interpreted as the time needed for the incoming wave packet to be scattered into the output lead. Let us add that the resonance lifetimes Rotter (2013), determined by the imaginary parts $\Gamma_{i}/2$ of the resonance eigenvalues ${\cal E}_{i}=E_{i}-i\Gamma_{i}/2$ of $H_{\rm eff}$, provide another set of time scales which characterize the scattering process. As long as distinct resonances do not overlap, i.e., \Gamma_{i},\Gamma_{j}\ll\mathopen{}\mathclose{{}\left\lvert E_{i}-E_{j}}\right\rvert for all $i$ and $j$, dwell times and resonance lifetimes are intimately related. However, for overlapping resonances, such direct association (see e.g. Lyuboshitz (1977)) breaks down.

II.2 Design Principles for Efficient Transfer

In a next step, we rely on earlier results which identified centrosymmetric Zech et al. (2014); Walschaers et al. (2013, 2015, 2017); Ortega et al. (2015, 2016) random networks as more efficient than unconstrained random assemblies, and specify our scattering model as given by centrosymmetric Hamiltonians of the form

[TABLE]

with $H_{\rm int}$ a centrosymmetric matrix which represents the hardwiring between the bulk sites of the network. $H$ commutes with the exchange operator $J$, where $J_{ij}=\delta_{i,N-j+1}$ in the site basis. $E^{\prime}$ gives the on-site energy of the input and of the output site of the network, coupled with strengths $v_{i}$ to the bulk sites, and with $V$ to each other. Only input and output site be coupled to input and output channel, respectively, and \mathopen{}\mathclose{{}\left|\rm in}\right\rangle and \mathopen{}\mathclose{{}\left|\rm out}\right\rangle be states fully localised on these respective sites, with \mathopen{}\mathclose{{}\left|\rm in}\right\rangle=J\mathopen{}\mathclose{{}\left|\rm out}\right\rangle. If we represent the input and output channel states by \mathopen{}\mathclose{{}\left|\Psi_{\rm in}}\right\rangle and \mathopen{}\mathclose{{}\left|\Psi_{\rm out}}\right\rangle, respectively, we can express the coupling operator $W$ in (1) as

[TABLE]

and $H_{\rm eff}$ in (1) thus takes the form

[TABLE]

With the additional choice $\Gamma=\Gamma^{\prime}$ Walschaers et al. (2017), $H_{\rm eff}$ becomes a non-Hermitian centrosymmetric matrix.

However, centrosymmetry alone is not yet sufficient to guarantee efficient transport features of the networks, as illustrated in Figs. 1 and 2, where the transfer probability $p_{\rm in\rightarrow out}$ and the dwell time $\tau_{\rm in\rightarrow out}$ are plotted for a typical, centrosymmetric random network of $N=8$ sites. In Fig. 1, we observe asymmetric (e.g. at $E\approx 1.5$, due to the interference of overlapping resonances) as well as symmetric resonance structures, with variable widths and strengths, some of which achieve $p_{\rm in\rightarrow out}=1$.

Yet, Fig. 2 highlights that large resonant transfer probabilities may be associated with – here undesirable – very long time scales. This is consistent with the narrow resonance structures in Fig. 1, which imply long resonance lifetimes,222Also note the negative value of $\tau_{\rm in\rightarrow out}$ at $E\approx 1.5$, which hints at the fact that (3) actually quantifies a phase shift rather than a genuine time. as well as with earlier findings Walschaers et al. (2013, 2015).

We therefore need to introduce another design principle, the dominant doublet condition – which requires that the network Hamiltonian $H$ exhibits eigenstates close to

[TABLE]

More formally, given the centrosymmetric matrix (6) in the symmetry eigenbasis (in which $J$ is diagonal),

[TABLE]

the dominant doublet condition imposes the existence of two eigenvectors \mathopen{}\mathclose{{}\left|\eta^{+}}\right\rangle and \mathopen{}\mathclose{{}\left|\eta^{-}}\right\rangle such that

[TABLE]

with $\epsilon>0$ close to zero.

With the help of lowest order perturbation theory Walschaers et al. (2013, 2015), (9) allows to derive analytical approximations for our quantities of interest. To start with, for the eigenvalues associated with \mathopen{}\mathclose{{}\left|\eta^{+}}\right\rangle and \mathopen{}\mathclose{{}\left|\eta^{-}}\right\rangle one finds:

[TABLE]

where $e^{\pm}_{i}$ and \mathopen{}\mathclose{{}\left|\psi^{\pm}_{i}}\right\rangle are the eigenvalues and eigenvectors of $H_{sub}^{\pm}$ in (8), respectively. Setting all $v_{i}=0$ in (4) implies $s^{\pm}=0$, and reproduces the result for a network shrunk to a dimer composed only of input and output site.

Combining Eqs. (9) and (10) we now determine the scattering matrix (1) which maps the input to the output channel as

[TABLE]

from which, with (10), we infer the transfer probability:333We will henceforth omit the index “${\rm in}\rightarrow{\rm out}$”, for ease of notation.

[TABLE]

where $\Delta s=s^{+}-s^{-}$. Similarly, $\tau(E)$ is derived from (10) with the help of (3), leading to the somewhat cumbersome expression:

[TABLE]

Inspection of (12) and (13) shows that transfer probability and dwell time sensitively depend on the energy $E$ of the incoming excitation. From (13), the energies which maximise the transfer probability $p$ are

[TABLE]

with $\overline{s}\equiv(s^{+}+s^{-})/2$ the average shift of the resonance energy. The transfer efficiency at these resonant energies follows as

[TABLE]

and the associated dwell times are

[TABLE]

Eq. (16) implies that, whenever \Gamma<\mathopen{}\mathclose{{}\left\lvert 2V+\Delta s}\right\rvert, there are two energies given by Eq. (15), at which the incoming wave packet is transmitted deterministically. Consequently, the profile of $p$ as given by (12) exhibits two well-separated resonances, as shown in the top panel of Fig. 3. For \Gamma=\mathopen{}\mathclose{{}\left\lvert 2V+\Delta s}\right\rvert, the two resonances merge, and the maximal transfer probability is achieved at $E=E^{\prime}+\overline{s}$ (Fig. 3, middle panel). The transfer probability starts to decrease as the ratio between $\Gamma$ and \mathopen{}\mathclose{{}\left\lvert 2V+\Delta s}\right\rvert increases beyond unity, as illustrated in the bottom panel of Fig. 3.

The on-resonance dwell times $\tau$ in (17) are governed by the parameter $\Gamma$. Indeed, it directly follows from (17) that, on resonance, $2/\Gamma<\tau<4/\Gamma$. Therefore, to obtain fast transport at large transfer probabilities, we must make $\Gamma$ as large as possible, under the constraint that \mathopen{}\mathclose{{}\left\lvert 2V+\Delta s}\right\rvert>\Gamma.

In the above perturbative approach, the impact of the bulk sites of the network is absorbed in the shifts $\Delta s$ and $\overline{s}$ of the dominant doublet. Apart from these shifts, the system then effectively behaves as a two-level system, comprised only of the input and output sites. Of course, there are additional resonances, at other energies than those of the doublet, as clearly visible in the top panel of Fig. 3. These must have very long life times, since the imaginary parts of the associated complex eigenvalues of $H_{\rm eff}$ determine their widths, and only emerge at higher orders of the perturbative expansion (10). While not accounted for in the analytical theory as presented here, such resonances are fully included in our numerical simulations and have no significant impact on the results presented in Sec. IV below.

III Statistical Analysis

Given our above analysis of the resonance structure with associated transmission characteristics of single realisations of disordered, finite networks, we can now address the statistical properties of the transfer probability which results from sampling over a distribution of such networks, under the above constraints of centrosymmetry and a dominant doublet. For this purpose, we assume that the direct input-output coupling $V$, as well as the channel coupling strength $\Gamma$ have identical values for all realisations of the network structure, and introduce the scaled system parameters

[TABLE]

Under this assumption, fluctuations of transfer probability and dwell time have their origin in the fluctuations of the relative level shift $\Delta s$ of the dominant doublet. That correction’s distribution function is given by López et al. (1981); Walschaers et al. (2013, 2015)

[TABLE]

with $\tilde{\sigma}$ and $\tilde{s}_{0}$ determined by the mean level spacing and the average coupling strength between the dominant doublet and the bulk states.

We saw above (recall top panel of Fig. 3) that efficient excitation transfer across the network can be achieved for $\Gamma$ as large as possible, yet under the condition of well-separated resonances associated with the dominant doublet states. In terms of our scaled system parameters, this latter condition reads

[TABLE]

and we can evaluate the probability to fulfill this condition, given (19), as follows:

[TABLE]

The density (23) is plotted in Fig. 4, for different widths $\tilde{\sigma}$ of the distribution of $\Delta\tilde{s}$, while we set $s_{0}=0$ (what is justified since $s_{0}\ll 1$ in general Walschaers et al. (2015), if the dominant doublet condition (9) is fulfilled).

There is a clear transition from deterministically separated dominant doublet resonances for $\tilde{\Gamma}\ll 1$, to overlapping resonances for $\tilde{\Gamma}\gg 1$. For $\tilde{\sigma}\ll 1$, this change occurs drastically at $\tilde{\Gamma}\approx 1$. However, the transition is smoothened, and shifted towards larger values of $\tilde{\Gamma}$, with increasing values of $\tilde{\sigma}$. Furthermore, overlap of the doublet resonances becomes likely for $\tilde{\Gamma}\sim 1+\tilde{\sigma}$, and large $\tilde{\sigma}$ is therefore desirable, since enhanced values of $\tilde{\Gamma}$ which still comply with (20) induce perfect excitation transfer on even faster time scales (recall our discussion in Sec. II.2 above).

After this first coarse-grained assessment of the statistics of the transfer efficiency via the competition between the scaled channel coupling and the distribution of disorder-induced relative level shifts of the doublet states, we now directly inspect the distribution of the transfer probabilities. To sample the latter at the dominant doublet energies $E^{\prime}\pm V+s^{\pm}$ of the closed system Walschaers et al. (2013), which follow from (15) with $\Gamma=0$, and are slightly detuned with respect to (15) for finite $\Gamma$, one needs to infer the statistics of

[TABLE]

as inherited from (19). This leads to

[TABLE]

what can be evaluated with the help of

[TABLE]

where $f(x_{j})=0$ and $f^{\prime}(x_{j})\neq 0$. This in (25) yields

[TABLE]

The condition $f^{\prime}(x_{j})\neq 0$ is violated on the edges of the domain $[0,1]$ of the probability distribution $P(p)$. The resulting distribution (27) is, therefore, only well-defined for $p\in[\delta,1-\delta]$. The divergences on the edges are an artefact of the power-law statistics (19) of $\Delta\tilde{s}$.

In Fig. 5, we compare the prediction of (27) (colour coded) to the transfer probability of a dimer (white line), which is obtained from (4) with all couplings $v_{i}=0$ between input/output and bulk sites, and also implies $s^{+}=s^{-}=0$. For the dimer, it is clear from our previous discussion, as well from its very structure which is determined by only two relevant coupling constants, $V$ an $\Gamma$, that the fastest transfer time scales are achieved for $\tilde{\Gamma}\approx 1$, as clearly displayed by the plot. In contrast, in the presence of bulk sites, and for a broad distribution of the relative doublet shifts as assumed in Fig. 5 by the choice $\tilde{\sigma}=10$, the transition between efficient and inefficient transport is pushed to values of $\tilde{\Gamma}$ enhanced by roughly one order of magnitude, an observation fully consistent with Fig. 4.

IV Comparison to Numerical Results

Let us finally confront our analytical – though perturbative (at lowest order) – prediction (25) with model calculations for open, disordered finite networks. We follow the modelling in Tomsovic and Ullmo (1994); Leyvraz and Ullmo (1996), and set

[TABLE]

where $\chi$ and $\xi$ parametrize the typical (root mean square) values of the model Hamiltonian’s (4) stochastically distributed coupling constants, and (28) implies that the bulk sites’ Hamiltonian $H_{\rm int}$ is sampled from the Gaussian orthogonal ensemble (GOE).444More details on this model can be found in Chapter 4 of Walschaers (2016). We now combine the results of Leyvraz and Ullmo (1996); Ergün and Fyodorov (2003) with the fact that for the intermediate sites the mean-level spacing is given by $\Delta=\pi\xi/(N/2-1)$ Brouwer (1997), and obtain that

[TABLE]

Furthermore, the results in Walschaers et al. (2015); Bohigas et al. (1993), imply that, to fulfil the dominant doublet condition (9), we must have

[TABLE]

what ensures that $\tilde{s}_{0}$ in (31) is negligibly small, as already assumed above.

In our simulations, we explicitly evaluate the scattering matrix element (1) and the transfer efficiency (2) from the input to the output channel, at energy $E=E^{\prime}+V+s^{+}$. $N$, $V$, $\xi$, and $\chi$ are all varied in the three numerical simulations displayed in Fig. 6, to broadly cover parameter space. To generate data, we choose different values for $\Gamma$, adopted to logarithmic scaling of the data represented in Figs. 6 and 7 (for practical purposes we chose $\Gamma_{i+1}=1.2\,\Gamma_{i}$). For each $\Gamma$, $10000$ random Hamiltonians (4) with (28) are diagonalized to extract the shift $s^{+}$. This allows the ultimate evaluation of $S_{{\rm in}\rightarrow{\rm out}}(E^{\prime}+V+s^{+})$, and Fig. 6 is output by the extrapolation method $\mathtt{SmoothKernelDistribution}$ of Mathematica, with the numerically obtained, discrete data points as input.

Fig. 6 shows good qualitative agreement between the theoretically obtained probability distribution (27) and the numerical results, what confirms that $\tilde{\sigma}$ and $\tilde{\Gamma}$ ultimately are the only parameters which control the statistics of transfer efficiencies. This physically implies that it is not the absolute energy scales set by $\chi$ and $\xi$ in (28) and (29), respectively, but rather the ratio of the coupling strengths (i.e., $\tilde{\sigma}$ and $\tilde{\Gamma}$) which govern the efficiency of the transfer process. The only absolute energy scale to affect the physics is the coupling $\Gamma$ to the external channels, because it directly controls the time scale of the transport.

Inspecting Fig. 6 in more detail does also reveal some quantitative differences between simulation data and theory. Several features, even though they emerge at the same values of $\tilde{\Gamma}$, are considerably broader in the numerical data. To verify whether these discrepancies are numerical artefacts induced by the here employed extrapolation method or due to essential physical features which are missed by our analytical treatment, we select three values of $\tilde{\Gamma}$, of different orders of magnitude, and compare analytical prediction (25) and numerical results for the density $P(p)$ in Fig. 6.

Clearly, the comparison is quantitatively excellent in all dynamical regimes (defined by the different values of $\tilde{\Gamma}$), what identifies the above discrepancies as a deficiency of the extrapolation method employed to produce Fig. 6. More importantly, this comparison also shows that our analysis, and (23) and (25) as the central results thereof, provide statistical guidance to optimise the transfer efficiencies of disordered networks, ultimately through $\tilde{\sigma}$ and $\tilde{\Gamma}$ alone.

V Conclusions

Our present contribution provides a scattering theoretical generalisation of earlier results Scholak et al. (2010); Zech et al. (2014); Levi et al. (2015); Scholak et al. (2011); Ortega et al. (2015, 2016) on optimal excitation transfer across finite disordered networks of dipole coupled two-level systems. Given disordered networks which exhibit the design principles already identified for closed systems – centrosymmetry and a dominant doublet – we have shown that, in order to achieve near-to-complete excitation transfer across such a network in minimal time, the coupling strength to the leads needs to be as large as possible, yet smaller than the dominant doublet splitting, as expressed by (20) above. In other words, the period of the coherent oscillation of the excitation between input and output site, as essentially defined by the dominant doublet splitting in the closed model Walschaers et al. (2013, 2015), has to be shorter than the decay time of the resonances associated with the doublet states. Consequently, correlation functions which quantify, e.g., the population of the input or of the output site, must exhibit (rapidly decaying) oscillations, as also familiar from typical 2D data recorded in experiments on photosynthetic light harvesting units Engel et al. (2007). Our results show in a rather transparent way that such separation of time scales, expressed by (20) and recently evidenced in experiments Dostál et al. (2016) on the molecular superstructure the FMO complex is embedded in, is a necessary condition for faithful transfer of the excitation: Increasing the channel coupling beyond the doublet splitting unavoidably reduces the transfer probability.

We have further shown that disorder-induced fluctuations of the doublet splitting as brought about by slow drifts and/or random distributions of the microscopic network hardwiring, and observed e.g. in the distribution of B800-B850 coupling strengths in experiments on the LH2 complex Hildner et al. (2013), allow to accelerate the transfer process under the above optimality condition (20): A broad distribution of relative shifts of the doublet levels allows for enhanced channel couplings, and, hence, for accelerated transmission, which would lead to a net loss of population transfer in the absence of disorder. Note that the same fluctuations also lead to an effective broadening of the condition for resonant excitation transfer when measured on a disordered ensemble.

Acknowledgements

M.W. thanks the German National Academic Foundation for financial support. R.M. is indebted to the Alexander von Humboldt foundation for support through a research fellowship. A.B. acknowledges funding through EU Collaborative project QuProCS (Grant Agreement No. 641277).