Unraveling excitation energy transfer assisted by collective behaviors of vibrations

We investigate how collective behaviors of vibrations such as cooperativity and interference can enhance energy transfer in a nontrivial way, focusing on an example of a donor-bridge-acceptor trimeric chromophore system coupled to two vibrational degrees of freedom. Employing parameters selected to provide an overall uphill energy transfer from donor to acceptor, we use numerical calculations of dynamics in a coupled exciton-vibration basis, together with perturbation-based analytics and calculation of vibronic spectra, we identify clear spectral features of single- and multi-phonon vibrationally-assisted energy transfer (VAET) dynamics, where the latter include up to six-phonon contributions. We identify signatures of vibrational cooperation and interference that provide enhancement of energy transfer relative to that obtained from VAET with a single vibrational mode. We observe a phononic analogue of two-photon absorption, as well as a novel heteroexcitation mechanism in which a single phonon gives rise to simultaneous excitation of both the trimeric system and the vibrational degrees of freedom. The impact of vibrations and of the one- and two-phonon VAET processes on the energy transfer are seen to be quite different in the weak and strong site-vibration coupling regimes. In the weak coupling regime, two-phonon processes dominate, whereas in the strong coupling regime up to six-phonon VAET processes can be induced. The VAET features are seen to be enhanced with increasing temperature and site-vibration coupling strength, and are reduced in the presence of dissipation. We analyze the dependence of these phenomena on the explicit form of the chromophore-vibration couplings, with comparison of VAET spectra for local and non-local couplings.


I. INTRODUCTION
Recent experimental and theoretical studies of the molecular structures present in biological light harvesting complexes have revealed the delicate interplay of electronic and vibrational degrees of freedom, and how these come together to orchestrate efficient transfer of photoexcitations in such systems [1][2][3][4][5][6][7][8]. Coherent beating patterns in nonlinear spectroscopy signals initially ascribed to long-lived electronic coherence in such systems [9] are now generally agreed to be a due to combination of electronic and vibrational coherence, with a key role played by coupling of the relevant electronic degrees of freedom to long-lived, underdamped vibrational modes of molecules [10][11][12][13][14]. This revelation has brought to light the subtle ways in which vibrational dynamics in molecular complexes can influence electronic and excitonic properties [15][16][17][18][19]. In this work we draw inspiration from these studies of molecular systems and ask whether vibrational degrees of freedom can exert other subtle influences on energy transfer in such complexes. In particular, can the presence of multiple underdamped vibrations in a molecular complex influence energy transfer dynamics in non-trivial ways?
With the recent growth of quantum technologies, controllable artificial quantum simulators have been developed to probe the underlying basic mechanism of the observed long-time coherences [20][21][22]. In Ref. [20], an engineered vibrationally assisted energy transfer (VAET) was experimentally demonstrated for an excitonic dimer emulated in a trapped-ion platform. In that work, not only was a one-phonon VAET process signified by a peak at the vibrational frequency being equal to the excitonic transition frequency unambiguously reported, but also unresolved peaks at smaller frequencies were found. It was suggested that the latter were due to multiphonon VAET processes. This provides further motivation for the study of multiphonon VAET processes. In this work we seek to ascertain the extent to which such multiphonon VAET can be resolved, and also the consequences of any collective behavior of the vibrations for excitonic energy transfer processes. Specifically, it is of interest to explore whether cooperative or interference effects might play a role in the vibrationally enhanced energy transfer. In addition, these systems offer the possibility of finding both the phononic analog of the well-known two-photon absorption [23][24][25], and the inverse phenomenon [26], in which one phonon might simultaneously excite two excitonic transitions, where the latter could be of different frequencies. The latter inverse situation offers a new twist with phonons relative to atoms, namely that in the context of VAET under our consideration, the question that we can ask is whether one phonon from a specific vibrational mode can simultaneously cause an excitonic transition and a vibrational transition in a different vibrational mode?
To address these questions, we consider here a donorbridge-acceptor trimeric chromophore system coupled to two vibrational degrees of freedom. We analyze the dynamics within a single electronic excitation subspace with explicit incorporation of the vibrational states, performing full numerical simulations for the excitation energy transfer probability from donor to acceptor via the bridg-ing chromophore under various conditions. We construct a two-dimensional spectral representation of the vibrationally assisted energy transfer (VAET) probabilities by scanning the frequencies of the two vibrations. These 2D VAET spectra allow identification of several mechanisms through which the vibrational modes can influence and/or enhance energy transfer. These mechanisms include both single mode VAET, and multi-mode VAET in which the two vibrational modes can cooperation and/or interfere. Detailed assignment of the VAET features is facilitated by calculation of the vibronic states resulting from the coupling of electronic and vibrational degrees of freedom, which shows that a number of the VAET features are correlated with the presence of avoided crossings in the vibronic energy spectra. In the weak sitevibration coupling regime we analyze the dynamical results with a perturbative analysis and use of double-sided Feynman diagrams [27]. For both weak and strong sitevibration coupling regimes we then investigate the dependence of the VAET features on exciton dissipation and vibrational temperature.
The context of this study is quantum emulation of excitonic energy transfer in ion traps and their use in elucidating the dynamic consequences of exciton-vibration coupling for energy transport in molecular excitonic systems. We focus first on the behavior when the vibrations are coupled locally to individual chromophore sites, as in the case of trapped ions coupled to transverse modes.The model trimeric system of primary interest in this work is generalized from a dimeric system studied previously with an experimental ion trap emulator [20] and starts from a Hamiltonian in which individual vibrational modes are coupled to Frenkel excitons on specific sites. We then extend this to study of the energy transfer dynamics induced by Hamiltonian coupling to vibrations that are correlated between sites, which naturally results from coupling to the longitudinal modes of trapped ions.
The trapped ion internal states in an ion trap emulation of molecular chromophores can be regarded as pseudo-chromophores. The ability to select different states and modulate their energies with external fields allows exploration of the effects of different energetic landscapes on energy transfer. However, the Hamiltonian for interaction of such internal ionic states with the external vibrational modes of the trapped ions has some subtle differences from the form of Hamiltonian relevant to studies of excitonic energy transfer in molecular aggregates such as natural light harvesting systems [28] or Jaggregates. We show here that despite these differences, the effective Hamiltonian that results from projection to the single electronic excitation subspace maps onto the standard form for excitonic energy transfer in the presence of vibrations that are correlated between different chromophores. We study the influence of these correlations on the effectiveness of the energy transfer, focusing in particular on the role of vibrations in enabling uphill energy transport in the excitonic degrees of freedom. The role of coherence in overcoming energy barriers has motivated discussion of rectification [29], while the coherent coupling of vibrational degrees of freedom in a quantum bath to an uphill gradient of excitonic states has been demonstrated to allow quantum ratcheting of energy transfer over long distances [30]. Clearly the phenomenon of VAET [20] is a prime enabler for such uphill energy transport. In this work we shall explore in detail the ways in which VAET processes facilitate and enhance energy transport, with particular emphasis on the role of VAET with correlated vibrational modes. We shall find that not only do resonant single phonon processes play a key role in enhancing energy transport, but that multiphonon processes can also be strong facilitators of energy transport.
The remainder of the paper is organized as follows. Sec. II introduces the model of the trimeric chromophore system with excitonic sites coupling locally to individual vibrational modes, e.g., as with a linear array of trapped ions coupled to transverse vibrational modes. In Sec. III, we summarize the procedure for numerical calculations of the 2D VAET spectra and the perturbation theory for analysis of the excitation energy transfer features. Full details of the latter are presented in the Supplementary Material [31]. Sections IV and V present and analyze results revealing the various VAET signatures in the weak and strong site-vibration coupling regimes, respectively. In Sec. VI, we present the vibronic energy spectrum and discuss the insights this offers for the VAET processes. We also discuss here the effects of the vibrational cross coupling terms induced by restriction to the single-exciton subspace of a Hamiltonian suitable for emulation of excitonic energy transfer with trapped ions. Sec. VII presents 2D VAET spectra and analysis for a trimeric system realizable for linear arrays of trapped ions with coupling only to longitudinal vibrational modes, resulting in a Hamiltonian with explicit correlations in the coupling of individual sites to the vibrational modes. Sec. VIII provides a summary and outlook for observation of the predicted VAET phenomena in trapped ion experiments, together with a discussion of the implications of this VAET study for understanding excitonic energy transfer in molecular systems.

II. THE EFFECTIVE MODEL OF A TRIMERIC CHROMOPHORE SYSTEM
We consider a donor-bridge-acceptor trimeric chromophore system coupled to two undamped vibrations, shown schematically in Fig. 1(a). The model can be described by the Hamiltonian where (settingh = 1) In the electronic Hamiltonian H s , the three sites correspond to the donor, bridge, and acceptor, respectively. Each site is modelled by a two-level system, with 2ω i the transition frequency between its ground and excited states, |g i and |e i , respectively. The Pauli operators σ (i) x,z are given by σ (i) x = |g i e| + |e i g| and σ (i) z = |e i e| − |g i g|. J ij is the coupling strength between adjacent ith-and jth-sites. We assume that the coupling between the first and third sites is vanishingly small. The vibrational mode with creation operator a † (b † ) and frequency ν a (ν b ) is coupled to the bridge (acceptor) site with coupling strength κ a (κ b ). In the quantum simulation context, such local coupling of sites to individual vibrations can be realized by coupling to transverse modes of a linear chain of trapped ions, so we denote this Hamiltonian by H tr . In addition to these undamped single-mode vibrations coupled to the electronic sites included in the above Hamiltonian, we will also incorporate dissipation effects by use of non-Hermitian terms in the Hamiltonian. Since the immediate context of this study is the emulation of excitonic energy transfer by trapped ions, these may be due to coupling to optical dephasing or to spontaneous emission, in addition to an additional damped vibrational environment as is usual in natural light harvesting systems.
We focus on the excitation energy transfer in a single electronic excitation subspace. An effective Hamiltonian within this subspace is obtained using the subspace pro-jection operator Ξ ≡ |egg egg| + |geg geg| + |gge gge|, usingH tr = ΞH tr Ξ. Explicitly, we find where |1 ≡ |egg , |2 ≡ |geg , |3 ≡ |gge , The couplings in the full Hamiltonian are depicted in Fig. 1(a) and resulting couplings of the effective Hamiltonian in the single excitation subspace are illustrated in Fig. 1(b). We note that, just as for a dimer [20], assistance from vibrations for energy transport becomes unnecessary with the transition frequencies of three sites are identical, since there are then no energetic difference between spatially separated sites. Our numerical calculations in this paper will focus on electronically uphill processes, as indicated schematically in Fig. 1(b). Two important remarks on the effective model in Eq. (5) are in order here. First, the counter-rotating terms in the XX-type interaction of Eq. (2) do not conserve the number or excitations and therefore do not survive the projection into the single electronic excitation manifold (since Ξσ − Ξ = 0). Second, it is evident that, in addition to the coupling between an electronic site and its directly connected vibrational mode [see Eq. (4)], the single-excitation effective model of Eq. (5) contains terms that couple a vibration to excited states of its unconnected sites, as indicated by dashed arrows in light green and blue in Fig. 1(b). These terms derive from the projection of the σ z operators through which sites 2 and 3 couple to their individual vibrations, i.e., Ξσ The Hamiltonian Eq. (1 differs formally from the common modelling of vibrational coupling of molecular chromophores, for which only coupling to excited state electonic states is included, i.e., κ a σ [28]. In the Born-Oppenheimer approximation, the ground electronic state of a molecular system is defined to be at the minimum of all relevant vibrational degrees of freedom and therefore there is no linear coupling to the electronic ground state. For a single chromophore, this coupling will just shift the overall energy of the system and both pictures are used for single chromophores in the literature. However, for two or more chromophores, we find that on projection to the single excitation subspace the σ z coupling introduces crosscoupling terms between excitations on different sites that mimic correlations between vibrations at different sites in a molecular system [10,28,[32][33][34][35]. In the following sections we shall see that these cross-couplings can give rise to enhanced collective phenomena in VAET. We note that the cross-couplings in Eq. (5) show both positive correlations between some sites, eg., between donor and acceptor sites for mode ν a and negative correlations between other sites, e.g., between donor and bridge sites for mode ν a . Consequently the overall effect of the correlations is not easily rationalized in terms of arguments for correlated dimers [28].
To illustrate the relevance of this study for understanding of natural photosynthetic systems, Table I shows typical values for the parameters considered in this work and compares them to the corresponding typical values for the Fenna-Matthews-Olson (FMO) light harvesting complex unit of green sulphur bacteria [36,37]. The first row shows typical values from the parameter ranges employed in this work and the third row shows corresponding typical values for the FMO system. The second row shows the result of setting our site energiesω 1 ,ω 2 , andω 3 to the natural value and scaling up the other parameters accordingly. It is evident that the range of parameters available in ion trap emulations scales consistently to the natural system, suggesting that analogs of the phenomena observed here might be present also in some natural light harvesting complexes (for more discussion of this, see the Supplementary Material [31]). Of particular relevance here are systems that have uphill regions in their landscape of Frenkel (site) exciton energies. This includes the FMO monomer complex [36], as well as the purple bacterium Rhodopseudomonas viridis [38] and the CP43 core antenna of photosystem II [39].

III. EXCITATION ENERGY TRANSFER PROBABILITY
We calculate the probability of finding an excitation in the acceptor, i.e., the third site, given an initial excitation localized on the donor, i.e., on the first site, as a function time, P 3 (t). Since there is no direct excitonic transfer from donor to acceptor, in the absence of coupling to vibrations the excitation is transferred via the bridging site. P 3 (t) is given by Here the unitary time evolution is for the whole system, i.e., U |1 1|ρ a ρ b U † , where |1 1|, ρ a , and ρ b are initial states of the chromophoric trimer and two vibrational  modes, respectively, and U = e −iHt is an evolution operator withH given by Eq. (5). Tracing over the vibrational degrees of freedom leads to the reduced excitonic system dynamics and then further taking the quantum average of the site number operator |3 3| gives the population at the acceptor site, P 3 (t). Numerical calculations of the transfer probability in Eq. (6) are performed assuming thermal initial states for the two vibrations, i.e., ρ a = e −νaa † a/k B Ta /T r a [e −νaa † a/k B Ta ] and , which measure the maximum and accumulated population during a given time period t f , respectively, as quantitative measures of the excitonic energy transfer efficiency.
Transforming to an the interaction picture with respect to free Hamiltonian [including H (e) 0 and the vibrational part (i.e., H where ∆ jk = λ j − λ k is the transition frequency between eigenstates and the forms for the coefficients A jk , B jk are given in Appendix A. We use fourth-order perturbation theory with respect to the site-vibration coupling to expand the evolution operator as and correspondingly obtain the transition probability with where the transition amplitudes are given by Unless otherwise stated, in the numerical calculations and perturbative analysis of these in the following two sections we shall study the maximum probability Max[P 3 (t)] during a given time period (0 ≤ t ≤ t f = 400ms), for a symmetric trimeric system with J 12 = J 23 = J andω 3 −ω 2 =ω 2 −ω 1 = ∆ in Eq. (5), coupled to two vibrations at a temperature larger (by a factor of 10) than these energies. This choice of parameters represents an energetically uphill process in the single excitation subspace. We shall present twodimensional (2D) VAET spectra obtained by evaluating Max[P 3 (t)] over the given time period as a function of the two vibrational frequencies ν a and ν b . These frequencies are given in units of the energy difference ∆ 31 between the eigenstates |e 3 and |e 1 of the electronic Hamiltonian H

IV. VAET SIGNATURES IN THE WEAK SITE-VIBRATION COUPLING REGIME
When the site-vibration coupling is weak, e.g., κ a , κ b < ∆, J, the electronic states of the trimeric chromophore system are only weakly perturbed and do not gain substantial vibronic character. In this situation the energy transfer processes are primarily excitonic in origin but assistance by vibrations that are coupled to the excitonic states can be still expected. Fig. 2(a), presents a 2D VAET spectrum in this regime and Fig. 2(b) presents a corresponding schematic diagram summarizing the energy transfer processes responsible for each of the main features of the 2D VAET spectrum. The perturbative analysis of these features is then summarized in Fig. 3. We now discuss the 2D VAET spectral features, starting with those due to single-mode VAET processes and then proceeding to the multi-mode VAET processes.

A. Single-mode VAET
The main features due to single modes in the weak sitevibration coupling regime are represented by the three vertical lines at ν a /∆ 31 = 1, 0.5, 0.25 and the horizontal lines at ν b /∆ 31 = 1, 0.5, 0.25. These signify resonant one-, two-, and four-phonon VAET processes, respectively, assisted by either the vibration ν a which is coupled to the bridging site 2 (vertical lines) or the vibration ν b which is coupled to the terminal, acceptor site (horizontal lines).
For a given number of phonons, e.g., one, two or four, the vertical line is more intense than the corresponding horizontal line, indicating a stronger impact of the bridging vibration ν a in assisting the energy transfer. This can be confirmed by our perturbative analysis as follows. For the resonant one-phonon VAET, perturbative expansion of the transfer probability (see Ref. [31] for full details) shows that the first order term, which is second order in the interaction Hamiltonian, is proportional to α 4 (W −a 31 ) * W −a 31 ∼ t 2 κ 2 a A 2 13 for transitions along the line ν a /∆ 31 = 1, and to along the line ν b /∆ 31 = 1. Here α ∼ 1 for the weak coupling regime and the coefficients A 13 , B 13 are given by Eqs. (A5) and (A6) of Appendix A, respectively. This VAET process is illustrated schematically by the Feynman diagrams in Fig. 3(a). Averaging this transfer probability over the thermal distributions of vibrational states for modes ν a and ν b [31], shows that even under conditions of identical site-vibration couplings (κ a = κ b ) and identical vibrational temperatures (k B T a = k B T b ), there will nevertheless be a higher probability for excitations along the vertical line than along the horizontal line, confirming the stronger impact of the vibration ν a that is coupled to the bridge site. This effect is clearly visible in the 1-dimensional top and right slices of Fig. 2(a), where it can be seen that the probability at the one-phonon VAET peak ν a /∆ 31 = 1 (top slice) is almost four times as large as that of the one-phonon VAET ν b /∆ 31 = 1 (right slice). Closer examination of the perturbative couplings shows that difference arises from the factor of 2 in A 13 relative to B 13 (see Eqs. (A5) and (A6)). A 13 is the matrix element of the electronic coupling to the bridge vibrational mode ν a , between the lowest excitonic state |e 1 and the highest excitonic state For the symmetric Hamiltonian we have e 1 |1 1|e 3 = − e 1 |2 2|e 3 = J 2 /2Ω 2 and e 1 |3 3|e 3 = J 2 /Ω 2 , resulting in a factor of 2 larger coupling for the electronic coupling to the bridge vibrational mode. When the vibrational frequency, either ν a or ν b , is equal to half of an excitonic transition frequency, the transition from |e 1 to |e 3 is still accessible but only via an intermediate state |e 2 . The excitonic transition is then accompanied by an absorption of two phonons. This single-mode two-phonon VAET appears in Fig. 2(a) at the vertical line ν a /∆ 31 = 0.5 and at the horizontal line ν b /∆ 31 = 0.5. We denote these processes as single mode two-phonon absorption (TPhonA) in Fig. 3(b). The corresponding transfer probabilities in the weak site-site coupling limit (J < ∆) are given by with α ∼ 1 for weak coupling as before, and A ij and B ij given in Eqs. (A1)-(A4). Further perturbative analysis with respect to J/∆ reveals that the transition |e 1 → |e 2 → |e 3 is a second-order process when the absorbed phonons are both borrowed from the bridging vibration ν a (since A 12 ∝ J ∆ and A 23 ∝ J ∆ , see Appendix A and Ref. 31), but becomes a third-order process when the two phonons are provided by the terminal mode ν b (since B 21 ∝ ( J ∆ ) 3 and B 23 ∝ J ∆ ). Therefore, by a similar argument to the onephonon VAET above, i.e., multiplying the above transfer probability by a prefactor 2n 2 a or 2n 2 b [31] to account for thermal averaging over the vibrational modes, it is then evident that just as in the one-phonon VAET process, the vibration (ν a ) coupled to the bridge site has a stronger impact on the two-phonon VAET processes than the vibration ν b connected to the terminal, i.e., acceptor site.
Another observable feature of the single-mode VAET in Fig. 2(a) is the four-phonon process at ν a /∆ 31 = 0.25 or ν b /∆ 31 = 0.25, which is described by the Feynman diagrams in Fig. 3(c). Comparison between the corresponding vertical and horizontal lines in Fig. 2(a) supports the conclusion above that the vibration ν a which is coupled to the bridge site has a stronger impact on the energy transfer than the vibration coupled to the terminal site. Since this four-phonon VAET is a higher order process than that considered in our perturbation theory, we do not provide an analytical expression for the transfer probability here.
We also find that the two-phonon VAET at ν a(b) /∆ 31 = 0.5 is dominant over the one-phonon VAET at ν a(b) /∆ 31 = 1). This is particularly marked for the single-mode VAET enabled by the vibration ν a that is coupled to the bridge site. This is evident from the vertical lines in Fig. 2(a) and the 1-dimensional slice located above this. While the dominance of TPhonA VAET is clearly visible for the bridging site vibration, it is also manifested to a lesser degree for the terminal site vibration (see 1-dimensional slice to the right of Fig. 2(a) where the integrated strength of TPhonA VAET is clearly stronger [31], despite a slightly lower maximal value). This dominance of the two-phonon VAET over one-phonon VAET is particularly marked in the weak coupling regime (recall κ a = κ b = 0.01kHz for Fig. 2), but will also be evident in the strong coupling regime results presented in Sec. V below. The reason for this dominance is the relatively high temperature considered here, i.e., k B T a(b) = 1.5kHz, which ensures availability of the required number of phonons for both one-and two-phonon VAET, while the probability of energy transfer increases with average phonon number. Thus, for the one-phonon VAET, the average phonon numbers are n a(b) ∼ 1 when ν a(b) /∆ 31 = 1, while for the two-phonon VAET, the average phonon numbers are n a(b) ∼ 2.4 when ν a(b) /∆ 31 = 0.5. In the latter situation there is a higher than required average phonon number for ν a (b), implying that the two-phonon processes are more likely (1 < 2.4/2 = 120%). This analysis also holds when the temperature becomes so low that the necessary number of phonons cannot be taken from the thermal state, e.g., for k B T a(b) = 0.5kHz, where the two-phonon VAET is still dominant relative to one-phonon VAET.
Finally, we note that the single-mode two-and fourphonon VAET features seen in Fig. 2 are in good agreement with the additional partially-resolved peaks observed in the recent experimental study of single-mode VAET in a dimer system emulated with trapped ions [20].

B. Multimode VAET
The collective behaviors of multiple vibrations enable unique signatures of VAET arising in our trimeric system, relative to those due to individual ones of single vibrations presented above. Such signatures are represented by the diagonal and anti-diagonal lines in Figs. 2(a) and (b). They include cooperativity and interference of the two vibrations. Of particular interest for the former is the manifestation of the phononic analog of Mayer's twophoton absorption [23], constituting the anti-diagonal lines, and the inverse of this that combines vibrational and excitonic transitions, which we refer to as heteroex-citation, and which constitute the diagonal lines.

Cooperative Two-Phonon Absorption VAET
The anti-diagonal lines in Fig. 2(a) signify cooperative processes in which phonons from both modes are involved in a VAET process. For example, the line ν a /∆ 31 +ν b /∆ 31 = 1 signifies a double-mode two-phonon cooperative VAET process. We designate this as a cooperative TPhonA process (c-TPhonA) in Fig. 2 This cooperative process with a simultaneous absorption of phonons from two vibrations in assisting the energy transfer constitutes a phononic analogue of two-photon absorption [23].
At the symmetric lattice point {ν a /∆ 31 , ν b /∆ 31 } = {0.5, 0.5} satisfying the resonance condition ν a = ν b = ∆ 21 = ∆ 32 , the transition from |e 1 to |e 3 proceeds via the bridge state |e 2 and is assisted by two strongly cooperative processes consisting of absorption of a single phonon from one mode, followed by a second phonon from the other mode, with perturbative transfer prob- , as illustrated in Fig. 3(d). The two transfer processes having absorption of two phonons from distinct vibrational modes in different orders can interfere with each other, giving rise to a double-mode two-phonon interference VAET. The Feynman pathways for this interference are shown in the right hand (red) subpanel of Fig. 3(d) and the associated perturbative expression for the probability is Fig. 4 shows the time dependence of P 3 (t) at (0.5, 0.5) (solid green line) together with corresponding time traces for nearby points along the (0.5, ν b /∆ 31 ) (dashed blue and yellow lines). This shows that the symmetric point (green line) is a point of destructive interference along the ν b axis, since it has a smaller maximal probability than the time traces of the nearby points. Similarly comparing the time traces of P 3 (t) for points along the (ν a /∆ 31 , 0.5) line shows that the point (0.5, 0.5) is a point of constructive interference along the ν a axis, since here the green curve has a higher maximal value than those of nearby points (red and cyan dot-dashed curves).
At all points away from the symmetric lattice point, i.e., along the rest of the anti-diagonal line ν a /∆ 31 + ν b /∆ 31 = 1 with ν a = ν b , the transition from |e 1 to |e 3 involves a virtual intermediate state, i.e., |e 2 rather than |e 2 as illustrated in the lower left subpanel of Fig. 3(d). These points therefore show relatively small transfer probability in Fig. 2. We point out that the upper portion of the anti-diagonal with ν b < ν a is no-ticeably more intense than the lower part where ν b < ν a . This is once again a consequence of the stronger impact of the vibration ν a coupled to the bridge site in assisting the energy transfer, relative to that of the vibration coupled to the terminal site. In situations for which the impact of the vibration ν b becomes dominant in the energy transfer, for example the Hamiltonian with site-correlated vibrational modes analyzed in Sec. VII, this effect is reversed.
Similarly, the second anti-diagonal line visible in Fig.  2(a), i.e., ν a /∆ 31 + ν b /∆ 31 = 0.5 at the bottom left, signifies a double-mode, four-phonon cooperative VAET involving two processes of two-phonon absorption each. Here the lattice point, i.e., {ν a /∆ 31 , ν b /∆ 31 } = {0.25, 0.25}, can host an interference between two transfer pathways involving different ordering of the twophonon absorptions, while the rest of this anti-diagonal line, i.e., ν a = ν b , shows a relatively low transfer probability due to the off-resonant nature of the intermediate virtual state. This constitutes a higher-order phononic analogue of two-photon absorption, which is illustrated in Fig. 3(e).

Heteroexcitation VAET
In addition to these cooperative VAET processes for which both vibrations contribute phonons to assist excitonic energy transfer, we also find evidence of VAET processes in which a phonon from one vibrational mode can simultaneously excite not only the electronic system but also the other vibrational mode. This is a new kind of cooperative mechanism of VAET, which we shall refer to as "heteroexcitation". It is evidenced by the diagonal lines in Fig. 2(a), which are also shown with their assignments as solid diagonal lines in Fig. 2(b). For example, the line ν a /∆ 31 − ν b /∆ 31 = 1 indicates that one phonon of vibrational mode ν a generates an electronic transition from site 1 to 3, together with absorption of a single phonon in mode ν b . The line ν a /∆ 31 − ν b /∆ 31 = 0.5 represents processes in which two phonons of vibrational mode ν a generate the same electronic transition, but now together with absorption of two phonons in mode ν b . An analogous interpretation applies to the line We note that the alternative heteroexcitation associated with the diagonal line ν b /∆ 31 − ν a /∆ 31 = 1, shown as the dotted diagonal line at the top left of Fig. 2(b), in which one phonon from vibrational mode ν b generates an excitonic transition together with absorption of a phonon of vibrational mode ν a is not observable in Fig. 2(a). This is because of the weaker impact of the vibration coupled to the terminal site. This vibration now has to provide one phonon to be absorbed by both the timer and the other vibration, which a significantly weaker process at this temperature. However this transition would appear on further increasing the temperature of mode b (k B T b ), so we also show the relevant Feynman diagrams for this process in Fig. 3(e). We note that while the reverse of two photon absorption, namely one photon exciting two atoms, has been discussed previously, it was assumed there that atoms of identical frequency are excited by a single photon [26], consistent with the larger wavelength of optical photons relative to atoms. The heteroexcitations seen here constitute a generalized phononic analog of that optical phenomenon, where now not only the energies, but also the nature of the two degrees of freedom being excited can be different. Another interesting interpretation of this process is that of redistribution of energy from one phonon reservoir to another, mediated by the electronic degrees of freedom.

C. Vibrational temperature effects
One interesting capability of artificial energy transport as studied with emulations using e.g., trapped ions, that is not possible in real molecular systems, is the ability to individually vary the effective temperature of different vibrational modes. Here we assess the effects of these temperatures on VAET features. Fig. 5 presents 2D VAET spectra at three different temperatures from that in Fig. 2, including also the presence of a temper-ature bias between the two vibrations in panels (b) and (c).
In the absence of a temperature bias, i.e, when the vibrational temperatures are equal (panel (a)), comparison with the higher vibrational temperature spectrum of Fig. 2 (k B T a = k B T b = 1.5kHz) shows that collective VAET features such as the cooperative behavior evidenced by the anti-diagonal and diagonal lines in the 2D spectrum become weaker as the vibrational temperature decreases, indicating a suppression of vibrationally assisted energy transfer processes. However the two-phonon VAET is still dominant over the one-phonon VAET, as discussed in Section IV A.
When a temperature bias between the two vibrations is present, as in panels (b) and (c), we find that increasing either k B T a [ Fig. 5 Fig. 5(c)] will enhance the transfer processes assisted by either of the vibrations coupled to the bridge or to the acceptor. This suggests that the weaker impact of the vibration coupled to the acceptor seen above could be enhanced by selectively raising the temperature of this mode in an emulation experiment.
Panels (b) and (c) of Fig. 5 show that the presence of a temperature bias across the two vibrations can also enhance heteroexcitations at |ν a − ν b | = ∆ 31 , ∆ 31 /2, relative to that seen for equal temperatures in Fig. 5(a). Such enhancement of heteroexcitations would increase with further increase of the temperature bias.
The time dependence of energy transfer probability in VAET processes is also strongly dependent on the vibrational temperatures.

D. Dissipative effects
In contrast to the assistance provided by undamped vibrations for excitation energy transfer discussed above, the decay of an electronic excitation at each site, resulting from e.g., spontaneous emission or coupling to damped vibrational environments, is expected to suppress energy transfer processes. Here we study the effect of such relaxation processes, using a non-Hermitian approach that in the single excitation subspace is equivalent to use of the Lindblad master equation with a relaxation operator [40] (see detailed analysis in the Supplementary Material [31]).
Including a non-Hermitian Hamiltonian term, the system under such dissipation is described by H eff = H − by Eq. (5) [31]. The average effect of dissipation on the excitation energy transfer is then obtained by repeating the 2D VAET spectral calculations withH replaced bỹ H eff . Fig. 6 shows the 2D VAET spectrum with dissipation given by parameters γ 1 = γ 2 = γ 3 = γ = 0.001kHz. We see suppression of all energy transfer processes, particularly those along the anti-diagonal and diagonal lines, relative to the no dissipation results in Fig. 2(a). As expected, the single-mode two-phonon VAET is the most pronounced VAET process in Fig. 6. Fig. 7(a) shows that for a specific VAET transfer process, e.g., the single mode TPhonA VAET at the resonant position ν a /∆ 31 = ν b /∆31 = 0.5, the time-dependent probability of finding an excitation at the acceptor, P 3 (t), is increasingly suppressed for all t as γ increases. To analyze which sites contribute to this suppression, Fig. 7(b) shows calculations with different dissipative parameters γ i at each site. Only small variations are seen, within the general trend that a strong dissipation at the donor site provides the greatest suppression (red dotted line), followed by having the strongest dissipation at the acceptor site (blue solid line). Interestingly, when the strongest dissipation is at the bridge site, the energy transport is most robust to the dissipation (yellow and green dashed lines).

V. VAET SIGNATURES IN THE PRESENCE OF STRONG SITE-VIBRATION COUPLING
As the site-vibration coupling strength increases, different VAET features emerge and the balance between single-and multi-phonon VAET processes changes. We explore these changes by considering larger coupling strengths κ a = κ b = 0.03kHz and 0.1kHz, summarized in Figs. 8(a) and 8(b), respectively. The top and right side slices in each of these plots are taken at ν a /∆ 31 = 0.746 (right slice) and ν b /∆ 31 = 0.746 (top slice).
In addition to the basic VAET features from Fig. 2(a) (where κ a = κ b = 0.01kHz), we now observe additional multiphonon VAET processes in Fig. 8(a) that involve three, five, and six phonons, indicated by vertical lines at ν a /∆ 31 = 1/3, 1/5, and 1/6, respectively. For the larger site-vibration coupling strength κ a = κ b = 0.1kHz shown Fig. 8(b), these vertical lines become more distinct and also start to shift noticeably away from the excitonic resonant transition frequencies. We also see that in the strong coupling regime, not only do the one-phonon (e.g., ν a /∆ 31 = 1) and two-phonon (e.g., ν a /∆ 31 = 0.5) VAET processes become more comparable in intensity, but also the impact of the vibration coupled to the acceptor, ν b , becomes comparable to that of the bridging vibration ν a . Thus, we now see local maxima at ν a /∆ 31 , ν b /∆ 31 = 1 and 0.25, in both the right and top slices of Fig. 8(a) and (b). The greater structure in these intensity patterns contrasts with the simpler structure obtained for weak site-vibration coupling in Fig. 2(a).
To further demonstrate the effects of strong sitevibration coupling on the VAET, we plot in Fig. 8(c) the time evolution of the P 3 (t) at the resonance point (ν a /∆ 31 , ν b /∆ 31 ) = (0.5, 0.5) for a two-phonon VAET process at various values of the coupling strength κ = κ a = κ b . When the electronic sites are decoupled from the vibrations, i.e., κ a = κ b = 0, the blue reference curve in Fig. 8(c) shows Rabi oscillations characterized by the transition frequency ν a = Ω = √ ∆ 2 + 2J 2 = 0.52kHz, with corresponding oscillatory period 2π/ν a ∼ 12ms (approx. eight cycles in each period of 100ms). When the  coupling is nonzero, we observe modulated Rabi-like oscillations that show slow oscillations superimposed on the fast oscillations with the frequency Ω. See, for example, the solid orange curve in Fig. 8(c), for which κ a = κ b = 0.01kHz. As κ increases to 0.05kHz (green curve) and beyond to 0.1kHz (red curve), the initial rise of P 3 (t) is faster and the first maximum higher. However, further increase of the site-vibration coupling strength beyond 0.1kHz reverses this trend. In the next section we shall see that this is a result of the formation of vibronic states with strong mixing of excitonic and vibrational degrees of freedom, giving rise to very different transition frequencies.

A. Vibronic states
In order to better understand the origin of the VAET features, we have calculated the energy spectrum for the trimer excitonic system coupled to the two vibrational modes including three vibrational levels [Eq. (5)] as a function of the scaled frequencies (ν a /∆ 31 , ν b /∆ 31 ), for specific values of the coupling strengths. This reveals the energies of the vibronic states formed as a consequence of the two exciton-vibration couplings. Fig. 9 shows the corresponding two dimensional vibronic spectrum for the case of coupling strengths κ a = κ b = 0.03kHz. The figure clearly shows the presence of avoided crossings that derive from the exciton-vibration coupling. For example, along the horizontal line ν b /∆ 31 = 0.5, whenever ν a /∆ 31 approaches a resonant transition frequency of the trimer excitonic system (i.e., zero detuning at ν a /∆ 31 = 0.25, 0.5, 1 as shown in Fig. 2), gives rise to an avoided crossing due to the site-vibration coupling. Each avoided crossing in the spectrum shown in Fig. 9 indicates a vibronic state, i.e., a mixing of the electronic and vibrational degrees of freedom [2].
Some of the avoided crossings in the vibronic energy spectrum correspond to VAET features discussed above. A perturbative analysis of the vibronic energies predicts the presence of avoided crossings at the degenerate states. Thus the specific avoided crossings magnified in Figs. 9(b) and (c) indicate the hybridized vibronic states ( ∼ = |110 ± |300 and |111 ± |301 (i, j, k in |ijk represent the excitonic state |e i , and vibrational occupation states |j and |k , respectively) that give rise to the one-phonon VAET feature along the vertical line ν a /∆ 31 = 1) in Fig. 2(a). Similarly, the states ∼ = |121 ± |301 at the avoided crossing in Fig. 9(d) are associated with the single-mode two-phonon VAET indicated by the vertical line, i.e., ν a /∆ 31 = 0.5, in Fig. 2. We also see vibronic states associated with the cooperative VAET features. In Fig. 9(e), the avoided crossing of states ∼ = |120 ± |301 corresponds to an intersection of the horizontal line ν b /∆ 31 = 0.5 and the off-diagonal line 2ν a − ν b = ∆ 31 in Fig. 2 that indicates the double-mode cooperative VAET.
We note that the avoided crossings in the vibronic spectrum become more pronounced as the coupling strength κ a or κ b increases, consistent with the perturbative analysis. This means that not only does the gap between the two adjacent levels increase, but also the shift from the excitonic resonant transition frequencies (e.g., ν a /∆ 31 = 1, 0.5, 0.25) will be larger. This trend is also visible in the cross-sectional slices in Fig. 8(a) and (b). Consequently, for given frequencies ν a , ν b , increasing either κ a or κ b to values so large that they are comparable with the excitonic energy differences will be expected to suppress energy transfer processes below values seen for smaller coupling. Indeed, this is consistent with the decrease in P 3 (t) seen for large κ a = κ b values in Fig. 8(c).

B. Effect of cross couplings in single excitation subspace
The VAET features presented above are based on the consideration of the effective Hamiltonian Eq. (5) derived as the single electronic excitation restriction of the model in Eq. (1) for the trimeric chromophore system. This trimeric model, generalized from an experimentally investigated dimer for an artifical excitonic system realized in a trapped ion system [20], contains interaction of the vibrations with both excited and the ground states, i.e., κ a σ . We saw that the resulting effective model in the single electronic excitation manifold has cross coupling terms, i.e., an interaction of a vibration with the excited states of unconnected sites.
Here we analyze the effects of the cross coupling terms on the excitation energy transfer. To isolate the effects resulting from these terms, we rewrite the effective Hamiltonian in Eq. (5) as with variable parameter ζ which interpolates between Eq. (5) for ζ = 1 and the usual single excitation manifold effective Hamiltonian for molecular excitons without cross coupling terms for ζ = 0. The latter case corresponds to the full Hamiltonian, Eq. (1), with the sitevibration coupling in Eq. (4) replaced by κ a σ Fig. 10 shows the two-dimensional VAET spectra for two symmetric trimeric systems with identical energy gaps (i.e., ∆ 21 = ∆ 32 ) that allow interference VAET to appear. Comparison of either the two left panels (a) and (c) with ∆ 21 = ∆ 32 = 0.332kHz, or the two right panels with ∆ 21 = ∆ 32 = 0.52kHz, shows that the cross coupling terms inH =H tr (ζ = 1) significantly enhance the energy transfer. For example, the one-phonon VAET at ν b = ∆ 31 = 1.04kHz, which shows high intensity forH tr (ζ = 1) (panel (b)), is considerably less intensē H tr (ζ = 0) (panel (d)).
Comparing now the left and right panels of either the upper (ζ = 1) or lower (ζ = 0) row shows the effect of modifying the energy barrier for both Hamiltonians. Thus the higher probabilities for the two-phonon VAET processes seen in panel (a) are due to the lower excitonic energy barrierω 3 −ω 2 =ω 2 −ω 1 = 0.3kHz which is more similar to the excitonic coupling J = 0.1kHz, than that of panel (b) for whichω 3 −ω 2 =ω 2 −ω 1 = 0.5kHz. The appearance of an interference VAET requires a specific condition, i.e., ∆ 21 = ∆ 32 . To isolate the interference VAET features we therefore present in Fig. 11 2D VAET spectra for two asymmetric systems that do not host any interferences. The upper row of  Fig. 12 where it is evident that they have a maximal values intermediate between those of symmetrically distributed neighboring points, implying an absence of interference at the symmetric points ν a = ν b = ∆ 21 = ∆ 32 . This is in contrast to the 2D VAET spectra for systems with identical energy gaps ∆ 21 = ∆ 32 in Fig. 10, where the interferences at the crossing point of vertical (ν a = 0.332kHz) and horizontal (ν b = 0.332kHz) lines in panel (a) and at the crossing point of ν a = 0.52kHz and ν b = 0.52kHz in panel (b) are clearly visible. The corresponding time traces (not shown), show destructive interference along the vertical lines [ν a = 0.332kHz (Fig.10  (a)) and ν a = 0.52kHz ( Fig.10 (b)], and constructive interference along the horizontal lines [ν b = 0.332kHz ( Fig.10 (a)) and ν b = 0.52kHz ( Fig.10 (b)]. We also present the corresponding results for the effective Hamil-tonianH tr (ζ = 0) in the lower row of Fig. 11, to emphasize the key role of the cross coupling terms in amplifying these cooperative VAET processes. The interference fea-  (color online) Time evolution of the probability P3 of the acceptor at the asymmetric point {νa, ν b } = {0.508, 0.343}kHz in Fig. 11(a) and {νa, ν b } = {0.343, 0.508}kHz in Fig. 11(b), together with nearby symmetric points {νa, ν b } = {0.508, 0.508}kHz,{0.343, 0.343}kHz.
All other parameters are the same as in Fig. 11.
tures are no longer visible here, confirming the critical role of the cross-correlated vibrations in enabling these quantum features.  Fig. 2(a) or Fig. 7.

VII. VAET WITH EXPLICITLY CORRELATED VIBRATIONAL MODES
The energy transfer in the trimeric chromophore system discussed in Sections II-VI relies on the assistance of two independent vibrational modes that are coupled to the bridge and acceptor sites, respectively [see Eq. (4)].
Experimental realization of this ideal Hamiltonian with the local and independent control of the site-vibration interaction for a trapped-ion quantum simulator requires that the ionic states be coupled to transverse vibrational modes. This is more challenging than coupling to the longitudinal vibrations, requiring stabilization with regard to the trapping frequency. Coupling to longitudinal modes allows instead access to normal modes that are linear combinations of local vibrations, exemplified for a three-ion system by the following Hamiltonian [41]: Note that here we have also included a direct excitonic coupling J 13 between donor and acceptor sites. This Hamiltonian describes the coupling of the trimeric system to the symmetric and asymmetric normal modes of vibration along the longitudinal axis of a linear chain of three ions. In contrast to the transverse coupling Hamiltonian of Eq. (1), we now have one vibrational mode, ν c (the symmetric stretch) that shows anti-correlated coupling to the donor and the acceptor and a second vibrational mode ν d (the asymmetric stretch) that shows a more complex correlated coupling with all three sites. Anticorrelated vibrations have been claimed to drive nonadiabatic electronic energy transfer in photosynthetic light-harvesting systems [10]. It is thus of interest to analyze the possibility of VAET processes in such a Hamiltonian possessing both correlated and anti-correlated sitevibration couplings. Following the procedure outlined in Section II to project the Hamiltonian onto the single electronic excitation subspace, we obtain the following effective Hamiltonian: We then perform the numerical simulations to evaluate the maximum energy transfer probability Max[P 3 (t)] for the effective modelH eff =H lng − iγ 2 (|1 1|+|2 2|+|3 3|), where we also include the effects of dissipation via non-Hermitian decay of chromophore excitations with parameter γ. The resulting 2D VAET spectra for the dynamics with and without dissipation in the weak coupling regime are presented in Fig. 13(a) and (b), respectively. Fig. 13(a) shows that this system with explicitly correlated couplings to the longitudinal vibrational modes generates a very different relative impact of the symmetric stretch ν c and the asymmetric stretch vibration ν d on the VAET spectra from that seen for the local transverse couplings in Sections II-VI. This is evident both in the ratios of single mode two-phonon to one-phonon VAET processes (horizontal and vertical lines), and in the variable intensities of the two-mode two-phonon VAET processes (antidiagonal lines). As we explain in detail below, this different impact reflects the fact that with the longitudinal modes, the asymmetric stretch ν d couples more strongly to the bridge site than to the donor and acceptor sites, while the symmetric stretch ν c couples only to the latter. This is true both in the full HamiltonianH lng and in the corresponding effective Hamiltonian, Eq. (12).
The single mode one-phonon VAET line for the asymmetric stretch (vertical line ν c /∆ 31 = 1) is more intense than that for the symmetric stretch (horizontal line ν d /∆ 31 = 1). The two-phonon lines ν d /∆ 31 = 0.5 and ν c /∆ 31 = 0.5 show an even greater disparity, to the extent that the two-phonon VAET for the symmetric stretch (ν c /∆ 31 = 0.5) is not dominant over the corresponding one-phonon process. The strongest feature in the VAET spectrum is now the single mode two-phonon absorption in ν d , i.e., the horizontal line ν d /∆ 31 = 0.5. This parallels the analogous dominance of the single mode two-phonon absorption for the mode coupled to the bridge site in the VAET spectrum of Fig. 2(a) (vertical line at ν a /∆31 = 0.5). The anomalous observation of the one-phonon symmetric stretch VAET for ν c being more intense than the two-phonon process results from the fact that this mode is completely decoupled from the excited state of the bridge site in both Eq. 11 and Eq. 12. This is quite different from not only the interaction of mode ν d in Eq. (12), but also that of modes ν a and ν b in Eq. 5, all of which include some coupling to the excited state of the bridge site in the effective Hamiltonian for the single excitation subspace.
More marked is the behavior in the anti-diagonal lines representing multi-mode VAET (corresponding to TPhonA in Fig. 2(b)). For example, the anti-diagonal line ν c /∆ 31 + ν d /∆ 31 = 1 shows significant intensity in the sector ν d > ν c but negligible intensity in the sector ν d < ν c . When dissipation is taken into account, Fig. 13(b) shows that all transfer processes are suppressed, similar to what is seen for the uncorrelated trimeric system above (see Fig. 6).

VIII. DISCUSSION AND CONCLUSIONS
We have systematically studied the phenomenon of vibrationally assisted energy transfer, VAET, in a donorbridge-acceptor trimeric chromophore system coupled to two vibrations over a range of coupling strengths. In this work we focused on two types of systems. The first derives from a Hamiltonian with uncorrelated local coupling, as would be obtained by coupling to transverse modes in a trapped ion quantum emulator. The second derives from explicitly correlated non-local coupling, as would be obtained by coupling to normal modes of longi-tudinal motion in a trapped ion quantum emulator. The parameters considered in this work are within the regime of current trapped-ion experiments [20] and in all cases we considered a parameter set ensuring energetically uphill transitions from both donor to bridge chromophores and bridge to acceptor chromophores.
In the case of local site-vibration couplings, we found a rich array of VAET phenomena going beyond the onephonon VAET observed previously with a trapped ion quantum emulator [20]. In particular, we also find clearly resolved signatures of two-and even four-phonon absorption processes in the 2D VAET spectrum at weak sitevibration coupling strength, while increasing the coupling strength introduces up to to six-phonon VAET processes. The two-phonon VAET processes constitute a phononic analogue of the well-known two-photon absorption [23] and we refer to them as TPhonA. They are found to be dominant for all coupling strengths, although the relative contributions of both one-and greater than two-phonon VAET processes do increase with coupling strength, gaining intensity from off-resonant contributions in the strong coupling regime. At all values of coupling strength, we find that for every VAET process the vibration coupled to the bridge has a significantly stronger impact than the terminal vibration, consistent with its central spatial location for energy transport across the chain of sites.
We also found that the two vibrations can give rise to multi-mode VAET processes in which they behave collectively, specifically via cooperation and interference that enhance the efficiency of energy transfer relative to that obtained from VAET with a single vibrational mode. This includes cooperative TPhonA in which the two phonons derive from different modes, possibly with different frequencies. We also observe an interesting phenomenon that is formally related to the reverse of this, namely processes in which a phonon from one vibrational mode simultaneously excites both the excitonic states of the trimeric chromophore system and the other vibrational mode. We term this process "hetero-excitation". A vibronic spectral analysis of the VAET features allowed detailed assignment and rationalization of the spectra, revealing the constructive effects of cross coupling terms in our derived effective model. Detailed analysis of transfer processes showing quantum interference was arrived at by considering a generalized asymmetric trimeric analog of the symmetric model for which the bulk of the numerical calculations were made.
The collective VAET features were found to be reduced but not completely suppressed by dissipative effects. We showed that they can however be enhanced by raising the temperature of the vibrational modes, as well as by increasing the strength of the site-vibration coupling.
In the case of explicitly correlated non-local sitevibration couplings, as would be obtained by coupling ions to longitudinal modes, we found generically similar VAET features but with quite different relative strengths. The most important parameter determining the integrated strength of excitonic energy transfer was seen to be the vibrational coupling to the bridge site of the chromophore system.
Our projection of the full Hamiltonian onto a single excitation subspace generates an effective Hamiltonian with induced cross-correlations in the effective site-vibration coupling that can be mapped onto excitonic energy transport for molecular chromophores coupled to correlated vibrational modes. The richness of the VAET spectra found here raise the intriguing question as to whether some of these VAET processes may be operating in natural systems. In particular, the results show the important role played by resonant vibrations in enhancing uphill energy transport. Such modes provide sharp features in the 2D VAET spectra, indicating a significant enhancement of the generic quantum ratcheting of energy transport that is derived from coherent coupling to a quantum vibrational bath [30]. As illustrated in Table I, the parameters considered in our study are scaled versions of parameters found in natural photosynthetic systems. This motivates further analysis of whether examples of the more complex VAET phenomena such as two-phonon absorption and heteroexcitations are present in any natural systems.
Finally, we emphasize that the cooperative behavior of multiple vibrations seen in this work act not only to enhance the excitation energy transfer but also demonstrate a rich set of VAET phenomena. Following the experimental observation of single-mode one-phonon VAET for a dimeric system in a trapped-ion quantum emulator [20], generalization of such experiments to three and more ions [22], as well as to other emulation platforms [21] appears feasible. We look forward to experimental verification of the predictions of VAET signatures for twophonon absorption and for heteroexcitations in emulations of a trimeric chromophore system.

ACKNOWLEDGMENTS
We thank Joseph Broz and Hartmut Häffner for helpful discussions.
Work at the University of California, Berkeley was supported by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences (BES) under Award # DE-SC0019376. Work at Sandia National Laboratories was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division.
Sandia National Laboratories is a multimission laboratory managed and operated by NTESS, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. DOE's NNSA under contract DE-NA-0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government. Here we summarize some key results from the analytical perturbative calculations provided in the Supplementary Material 31 that are used in the main text. In the interaction picture, the coefficients A jk , B jk that determine the perturbative factors discussed in Section IV are given explicitly as functions of J, ∆ and Ω as with ∆ =ω 2 −ω 1 =ω 3 −ω 2 and Ω = √ ∆ 2 + 2J 2 . The energy transfer probability at the acceptor can be written as . Note that in the weak site-site coupling limit, i.e., J ∆, we have α → 1 and β, γ → 0 and the excitonic eigenstates |e 3 and |e 1 can then be approximated by |3 and |1 , respectively [31].

Appendix B: Convergence
As demonstrated in the main text, a large transfer probability can be realized by increasing either the vibrational temperatures k B T a and k B T b , or the sitevibration coupling strengths κ a , κ b . However higher values of temperature or site-vibration coupling might cause convergence issues if the truncation number N of each vibrational Fock space is not large enough. To address this issue and confirm convergence of the numerical results, we show in Fig. 14 the transfer probability for various values of N , at vibrational temperatures k B T a = k B T b = 1.5kHz for several values of the sitevibration coupling. When the coupling strength is weak, e.g., κ a = κ b = 0.01kHz in Fig. 14(a), and when it is increased to 0.05kHz and 0.1kHz as in Fig. 14(b) and (c), respectively, the results are already convergent for N = 15 . This provides evidence for the accuracy of our results in Fig. 2(a). Increasing the coupling strength to an ultra-strong regime [42], e.g., to κ a = κ b = 0.5kHz, which is comparable to the excitonic transition frequency, requires a value as large as N = 30 to achieve a good convergence, as shown in Fig. 14(d). Similarly, convergence of the spectra in the very low vibrational frequency regime (i.e., smaller ν a(b) /∆ 31 ) might be expected to require a larger value of N . However, the general form of the spectra in this region is already converged at N = 15, as shown explicitly in the Supplemental Material [31].

SUPPORTING MATERIAL FOR "UNRAVELING EXCITATION ENERGY TRANSFER ASSISTED BY COLLECTIVE BEHAVIORS OF VIBRATIONS"
In this supporting material we provide details of our perturbation-based analysis that helps to interpret our findings in the main text. We also present additional numerical results that demonstrate features of vibrationally assisted energy transfer (VAET) in the trimeric chromophore system, including plots of integrated probability, t f 0 P 3 (t)dt, as an alternative measure of the transfer efficiency. Finally, we also illustrate the equivalence between the non-Hermitian approach that we employ to study dissipative effects, and a Lindblad master equation approach.

PERTURBATION-BASED ANALYTICAL RESULTS
In the main text we focus on the symmetric version of the effective model because analytical results for this model can be developed. In this section we provide details of the perturbation-based analysis of this model.

The symmetric donor-bridge-acceptor trimeric chromophore model
Relative to the effective three-level model derived in the main text, its symmetric version under our consideration means identical site-site couplings and effective energy gaps, i.e., J 12 = J 23 = J andω 3 −ω 2 =ω 2 −ω 1 = ∆. This constraint on the effective level energies also implies ω 3 − ω 2 = ω 2 − ω 1 = ∆/2, namely, identical difference between frequencies of two adjacent sites in the original model. Under these conditions, the effective Hamiltonian given by Eq. (5) in the main text becomes where the unperturbed part is and the interaction Hamiltonian is given by In an interaction picture with respect to the unperturbed Hamiltonian H 0 given by Eq. (B2), the whole Hamiltonian in Eq. (B1) becomes where U † 0 = e iH0t and An alternative way of deriving ξ(t) and χ(t), in addition to the Baker-Campbell-Hausdorff (BCH) formula, is to exploit the operator decomposition in the eigenstate basis {|e j }, e.g., where eigenstate states |e j satisfy H The eigenstates |e j are given explicitly by and the eigenvalues λ j are with Ω = √ 2J 2 + ∆ 2 . The transition frequencies, i.e., ∆ jk = λ j − λ k , between eigenstates are We then obtain the following expression for the electronic operator ξ(t) [defined in Eq. (B7)], in the interaction picture, where A jk e i∆ jk t are the matrix elements of the electronic operator in the interaction Hamiltonian, i.e., (|2 2|−|1 1|− |3 3|), in the eigenstate basis. The A jk are explicitly, Similarly, the electronic operator χ [defined in Eq. (B8)], in the interaction picture, becomes where Therefore we have the full expression for the Hamiltonian in the interaction picture:

Perturbed evolution operator
In the regime of weak site-vibration coupling, e.g., {κ a , κ b } < max{J, ∆}, the evolution operator in the interaction picture reads with T being time-ordering operator and U (0) I , U I , and U I ] already allows us to uncover collective behavior, e.g., not only cooperative behavior but also interference of two vibrational modes, higher-order terms are neglected in the following.

Transfer probability
The probability of finding an excitation at the acceptor (i.e., the third electronic site) of the trimeric chromophore system is where transition amplitudes at different orders are defined as In the first line of the equation that defines P 3 (t), U † 0 |3 3|U 0 is the projection operator |3 3| in the interaction picture. The initial state of the whole system is |1 1|ρ a ρ b , independent of the chosen picture. We mention that in the last line of Eq. (B26) there are no cross terms like (A (0) ) † A (1) that contain an odd number of vibrational operators and therefore cannot survive the trace (e.g., T r a,b [· · · ρ a ρ b ]).
The transfer probability can be written as a sum of probabilities at different orders, namely, 3 , where These expressions can be further simplified when the weak site-site coupling regime (J < ∆) is considered, as shown below.

Interaction amplitudes
To gain insight from the perturbative expressions of transfer probability above, we define the nth-order interaction amplitude (due to the coupling between the system and the vibrational modes), as W q1x1,··· ,qnxn jk1,··· ,kn−1kn = κ x1 · · · κ xn (X 1 ) jk1 · · · (X n ) kn−1kn where x i ∈ {a, b} and X i ∈ {A, B}. q i ∈ {+, −} represents an emission or absorption of a phonon and j, k i ∈ {1, 2, 3} denote electronic eigenstates (|e 1 , |e 2 , or |e 3 ), respsectively. In addition to the site-vibration coupling strength κ a , κ b for vibrational mode ν a and ν b respectively, the interaction amplitude further depends on A jk , B jk , which are matrix elements of the electronic interaction operator in the eigenstate basis and given above.
In the framework of the fourth-order perturbation theory considered in our work, the first four interaction amplitudes are needed for a full expression of the transfer probability. Here we explicitly write down the first three of them that are relevant in this supplemental material, namely

Weak site-site coupling
In this subsection we specialize to a weak site-site coupling J < ∆ (or J < Ω) regime. This assumption allows us to approximate initial and final states in the site basis roughly by the eigenstates, i.e., and the transition amplitude appearing in Eq. (B26) becomes Therefore the probability in Eq. (B32) becomes The amplitude A (0) is zero and does not contribute to the probability, which is however, not true when J is not small as shown in the next subsection. This implies that, without assistance from the vibrational modes, it is impossible for an excitation to transfer from the initial donor state |1 to the target acceptor state |3 .
a. One-phonon processes By using the amplitude A (1) in Eq. (B43), we obtain the first-order probability as Here we used (W +a ij ) * = W −a ji and (W +b ij ) * = W −b ji and the average phonon numbers of two vibrational modes ν a and ν b are given by respectively. We mention that terms like (n a + 1)(W +a 31 ) * W +a 31 that do not conserve energy have been neglected in P (1) 3 above. This direct transition from |1 to |3 results from the assumption of weak J under which we have used |1 ≈ |e 1 and |3 ≈ |e 3 . It is obvious the transition becomes forbidden when J → 0 due to . When the transition becomes resonant, i.e., ν a = ∆ 31 or ν b = ∆ 31 , the terms in P 3 , corresponding to a one-phonon process assisted by either mode ν a or ν b , can be further simplified as b. Two-phonon processes The two-phonon processes refer to transitions from the donor |1 to the acceptor |3 via an intermediate state |2 . The corresponding second-order probability is where The coefficients in front of interaction amplitudes accounting for the effect of temperature come from averages of products of four vibrational operators, evaluated via Wick's theorem [43], aaa † a † = 2(n a + 1) 2 , a † a † aa = 2n 2 a , a † aa † a = n a (2n a + 1), (B54) aa † aa † = (n a + 1)(2n a + 1), a † aaa † = aa † a † a = 2n a (n a + 1), aa † bb † = (n a + 1)(n b + 1), a † ab † b = n a n b , aa † b † b = (n a + 1)n b , a † abb † = n a (n b + 1).
Substituting the expressions for W ik,kl from Eqs. (B38) into the expressions for P , we find that the resonance conditions for the two probabilities are different. For P (2,2) 3 the resonance condition is ν a = ∆ 31 and/or ν b = ∆ 31 , whereas for P is much more relevant to the two-phonon VAET processes shown in the main text, we examine it further. We simplify each term in P (2,1) 3 by considering resonant transitions, i.e., ν a , ν b → ∆ 21 = ∆ 32 , in which case, 6. Strong site-site coupling The weak site-site coupling J < ∆ considered above does simply the problem, but it is, however, actually not necessary for obtaining perturbation-based analytical results, which only require a small κ a and κ b . In this section we therefore present results that go beyond weak site-site coupling.
To calculate the transition amplitudes, e.g., A (i) , and therefore probability P 3 (t) = P 3 + P (see Sec. B 3), we need to first transform the initial and target states from the site basis into the eigenstate basis. By using eigenstates defined in Eqs. (B11), (B12), and (B13), we have where We then obtain an analytical expression of the probability of finding an excitation in the final state |3 (acceptor) transferred from the initial state |1 (donor).
a. P (0) 3 By first deriving the transition amplitude A (0) = 3|U 0 |1 in Eq. (B27), we obtain the zeroth-order probability As expected, the transition amplitude is independent of interaction with vibrational modes and therefore is consistent with unitary evolution.
b. P (1) 3 The first-order probability P . (B69) The first term on the right-hand-side is obtained as where |A (1)

A remark on the FMO parameters
The parameters we use for simulations in the main text are those relevant to trapped-ion simulations. We show that if these are scaled, they coincide reasonably well with the relevant parameters in typical photosynthetic complexes, specifically the FMO complex. However, in Table 1 in the main text, it can be seen that the site-site couplings (J) in FMO are stronger than the scaled values by a factor of 2 to 5. Therefore, to verify that our observations are valid even in these regimes of stronger site-site couplings, we perform the simulations with the FMO parameters in the line 3 of Table 1 and present the result in Fig. 19. Indeed, the results look similar to Fig. 2 of the main text, with prominent VAET and multi-phonon VAET features.

NON-HERMITIAN APPROACH FOR DISSIPATION
In this section we justify the non-Hermitian approach that we use to study effects of the dissipation. jump term, we have with H eff = H − ih 2 3 j=1 γ j |e j e|. In the main text, we restrict attention to the single excitation subspace, and so projecting this equation into the single excitation yields: