Engineering vibrationally-assisted energy transfer in a trapped-ion quantum simulator

Many important chemical and biochemical processes in the condensed phase are notoriously difficult to simulate numerically. Often this difficulty arises from the complexity of simulating dynamics resulting from coupling to structured, mesoscopic baths, for which no separation of time scales exists and statistical treatments fail. A prime example of such a process is vibrationally assisted charge or energy transfer. A quantum simulator, capable of implementing a realistic model of the system of interest, could provide insight into these processes in regimes where numerical treatments fail. We take a first step towards modeling such transfer processes using an ion trap quantum simulator. By implementing a minimal model, we observe vibrationally assisted energy transport between the electronic states of a donor and an acceptor ion augmented by coupling the donor ion to its vibration. We tune our simulator into several parameter regimes and, in particular, investigate the transfer dynamics in the nonperturbative regime often found in biochemical situations.


I. INTRODUCTION
Charge and energy transfer are essential to many important processes in chemistry, biology, and emerging nanotechnologies. Such transfer processes often occur in noisy thermal environments that strongly modify the transfer dynamics and, in some cases, even improve the transport efficiency or robustness [1][2][3][4][5][6].
In these processes, the dominant sources of fluctuations and noise are often intramolecular vibrations and solvent dynamics [4]. Understanding the influence of this molecular environment on transport dynamics requires solving complex and often fully quantized models that become intractable to theoretical treatments even for systems of moderate size. Moreover, direct experimental studies of these vibrationally assisted energy transfer (VAET) phenomena are difficult to interpret since the underlying mechanisms cannot be isolated. A quantum simulation of such models, on the other hand, offers the possibility to both isolate and control the interesting aspects of the underlying mechanisms.
In a simple model featuring vibrationally assisted energy transfer, the environment consists of a thermalized vibrational degree of freedom that can assist the exchange of quantized excitations between a donor and an acceptor site [see Fig. 1(a)]. In general, these sites exhibit different energies such that transfer only occurs if the excess energy is taken up or provided by the vibrationas such, the environment assists in the transfer process. This model captures the important features of vibrationally enhanced phenomena, such as the dependence of transfer efficiency on the spectral properties and temperature of the environment.
Here, we demonstrate VAET in isolation and under fully controlled conditions. We encode the VAET process in a trapped-ion quantum simulator, where energy transfer between the electronic states of two ions is enhanced when coupled to an environment in the form of the thermal vibrational motion of the ion crystal. We observe the hallmark feature of VAET, namely, the strong dependence of the energy transfer efficiency on the temperature and the spectral characteristics of the environment. In addition, we tune our quantum simulator into nonperturbative parameter regimes, similarly to what is encountered in models of biochemical processes.
The model we consider couples two-level energy sites [donor (d) and acceptor (a)] with strength J described by x [see Fig. 1(a)]. In the absence of additional interactions, this coupling will cause a single excitation to oscillate between the sites with frequency J. A relative energy detuning ∆, represented by a term ∆ 2 σ (d) z , imposes an energy cost to move an excitation between the sites, therefore inhibiting complete transfer of excitations. A long-lived vibrational mode in the environment is modeled as a harmonic oscillator with frequency ν eff and is coupled to the sites in the form of κ 2 σ z (a + a † ). Quantum mechanically, the role of the environment may be understood as providing an extra degree of freedom, which helps to satisfy energy conservation in the transfer process. For instance, if ν eff ≈ √ ∆ 2 + J 2 , intersite transfer can occur, provided the environment changes its vibrational quantum number by 1 [see Fig. 1(a)]. Higherorder processes also occur: If ν eff ≈ √ ∆ 2 + J 2 /k for any integer k, intersite transfer occurs and the environment changes by k vibrational quanta. A classical analogue to this process is that the intersite energy difference ∆ Without the presence of an environment (κ = 0, left drawings), the transition probability from donor to acceptor states is attenuated in the presence of an energy barrier ∆. By coupling to an environment (κ > 0, right drawing), an excitation can move between donor and acceptor sites by exchanging energy with a phononic environment. Time traces illustrate the three situations: (1) ∆ = 0: black trace (theory) and data points (•). (2) ∆ > J, without assistance from the environment: blue trace (theory) and data points (×). (3) ∆ > J with assistance from the environment (VAET process): red trace (theory) and data points ( ). (b) Schematics of the ion trap and laser beams generating the simulated Hamiltonian. Internal levels of the ions (blue spheres) serve as energy sites. A laser beam illuminating both ions generates the site-site coupling with strength J. A localized beam generates the coupling to the environment with strength κ and controls the detuning ∆. (c) Laser tones (orange arrows) generating the simulated Hamiltonian. Vertical black lines represent available transitions, with the atomic resonance of the respective ion in the center. The detuning ∆ is introduced via an ac-Stark shift by imbalancing the relative power of the two tones of the local beam, as indicated in the figure. is modulated, parametrically moving the excitation between the sites. The resulting Hamiltonian for the VAET model is then ( = 1) The remainder of this article is organized as follows. After discussing details of the implementation in Sec. II, we present our experimental results in Sec. III. Here, we first study the case where the detuning ∆ is larger than the site-site coupling J such that appreciable energy transfer can only occur with the assistance of the environment. Then, we make ∆ comparable to J and operate the simulator in a regime where all Hamiltonian terms compete with each other. We then initialize the environmental mode in a thermal state with variable temperature; as such, this model corresponds to a biochemical scenario where the coupling of one pigment molecule to its environment is dominated by coupling to one slowly relaxing harmonic mode. We finish by providing an outlook and conclusions in Sec. IV.

II. EXPERIMENTAL IMPLEMENTATION
We implement this model in a trapped-ion quantum simulator consisting of two trapped 40 Ca + ions confined in a radio-frequency Paul trap [ Fig. 1(b)]. Ca + has a ground S 1/2 orbital and a metastable D 5/2 orbital. The magnetic substates |S (m j = 1/2) and |D (m j = 1/2) form a qubit, addressed by an optical transition near 729 nm [15]. The electronic states of the two ions play the role of the energy sites. The state |DS (|SD ) corresponds to a single excitation localized to the donor (acceptor) site. The two-ion crystal has six normal vibrational modes, only two of which are relevant to implement the Hamiltonian in Eq. (1): the axial stretch mode with ω ax ≈ 2π × 1.3 MHz and the radial rocking mode at ω r ≈ 2π × 2.1 MHz. The rocking mode serves as the thermally occupied bosonic environment in the simulation, while the axial stretch mode mediates the coupling between the donor and acceptor sites. The additional vibrational modes of the ion crystal generate Hamiltonian terms that rotate with frequencies of at least 2π × 300 kHz. We operate our simulator in the regime where all Hamiltonian parameters are of order a few kHz and therefore neglect these additional couplings.
The amplitude of each term in the model Hamiltonian Note that the theory estimates include the effect that the detuning ∆ fluctuates around its mean due to relative intensity noise of 0.02 as well as due to variations of the ion-laser coupling strength via finite-temperature effects. In (a), the standard deviation of the corresponding distribution is 2π × 0.23 kHz, whereas in (b) it is 2π × 0.1 kHz. Solid lines in the spectral plots represent a numerical solution where a small frequency offset adjusts for a systematic bias of the vibrational frequency measurements. In all plots, the statistical error is smaller than the markers. The measured data points in the spectral plots are connected with a dashed line to guide the eye.
[Eq. (1)] is controlled by adjusting the strength and frequency of various laser tones as summarized in Fig. 1(c).
In particular, the site-site coupling J 2 σ is implemented by a Mølmer-Sørensen quantum interaction [16] via the axial vibrational mode. This interaction is generated by applying a global laser beam with tones detuned from the qubit transition by ±(ω ax − δ ms ), where [16], where η ax ∼ 0.05 is the Lamb-Dicke parameter for this mode, and Ω G 1,2 are the Rabi frequencies of the two laser tones. For calibrating the simulator, J is measured independently by operating the simulator with ∆, κ = 0. We extract J by fitting the population transfer as a function of τ sim to a sinusoid. Note that in the single excitation manifold, spanned by |SD and |DS , the σ interaction is the same as the energy exchange interaction σ + , which is more commonly used in the context of charge and energy transfer in biochemical systems.
The site-environment coupling is engineered via a bichromatic laser beam localized to the donor ion. The two tones of this beam are detuned from the optical tran-sition by −ω r /2 and ω r /2 + ν eff [17]. Note that ν eff is defined as the difference between the ion-crystal rockingmode frequency ω r and the frequency splitting between the two laser tones. This generates the effective interac- Here, Ω L i is the on-resonance coupling between the |S and |D states generated by the ith tone of the local beam. The parameter η r = 0.039(1) is the Lamb-Dicke parameter for the radial rocking mode. Experimentally, we adjust the laser powers of each of the tones to reach the desired coupling κ.
To calibrate the coupling to the environment, κ, we measure the Rabi frequencies Ω L 1 and Ω L 2 directly in a two-stage procedure. First, both tones in the local beam are blue-shifted in a double-pass acousto-optic modulator (AOM) by ω r /2, such that one tone is resonant with the carrier transition. In that configuration, we measure the frequency of Rabi oscillations and extract Ω L 1 . Then, we shift both tones red by ω r /2 and again measure Rabi oscillations, extracting Ω L 2 . In order to maintain a stable environment frequency ν eff , we actively stabilize the radial trap frequency to within 2π × 250 Hz over the duration of one time scan (several minutes), following a method detailed in Ref. [18]. The rocking-mode frequency ω r is determined via optical spectroscopy on the |S → |D transition. However, ac-Stark shifts from the probe beam make a small correction on the order of a few hundred Hertz necessary. The systematic effect could be avoided by directly exciting the vibration with either a modulated optical radiation pressure force [19] or by measuring all center-of-mass modes using electrically oscillating fields.
We control the detuning ∆ by adjusting the power balance of the two tones on the local beam generating the site-bath coupling. The resulting detuning is measured by setting ν eff = 2π × 30 kHz, sufficiently far offresonant such that the coupling to the vibrational mode can be neglected. In the case ∆ J, the dynamics can be fit to the one given by the simplified Hamilto- z , extracting J and ∆. When ∆ J, the population transfer is too small to be fit for the detuning. In this case, we measure ∆ in a Ramseytype experiment [20]. First, a global π/2 pulse is applied to both ions. After an interrogation time τ , a second π/2 pulse is applied and the parity P is recorded. The operator P is defined as P (|SS ) = P (|DD ) = 1 and P (|SD ) = P (|DS ) = −1. The parity P (τ ) oscillates with frequency ∆. In both cases, an analytical correction is applied to account for the small change in Stark shift (up to 2π × 200 Hz) arising from moving one laser tone by 2π × 30 kHz as necessary to decouple the bath from the dynamics.
The actual experimental procedure for the simulation is as follows: We start by Doppler cooling all the vibrational modes of the ion string to a mean occupation number of 6-12, followed by optical pumping both ions to the state |SS . We further cool the axial center-of-mass and stretch modes to the ground state via resolved sideband cooling [21]. The assisting mode is then prepared via resolved sideband cooling [21] to an adjustable mean thermal occupation 0.04 n 12. The thermal occupation value is chosen by varying the duration of the cooling process. For small temperatures, the average population is extracted by comparing red-and blue-sideband excitation of the mode [22], while for higher temperatures, the sideband strength is compared to the carrier transition strength. The donor is then excited via a local rotation, leading to the combined electronic state |DS . Then, the local and global laser beams generating the model Hamiltonian are applied for a time τ sim . Finally, the combined electronic state of both ions is measured by recording the fluorescence with a charged-coupled device (CCD) camera [23]. For each parameter setting {J, κ, ∆, ν eff }, the simulation is run 100-500 times. The transfer probability is then given by the population in the state |SD . For all data, we report the conditional probability P acc that the system has undergone state transfer, i.e. P acc = P SD /(P SD + P DS ). This corrects for an average population loss from the simulation subspace of 10% arising from imperfect state preparation and incoherent excitation of optical transitions due to spectral impurities in our qubit laser at 729 nm.
In addition to probing the time dynamics of the VAET process by varying the simulation time τ sim , we also investigate its spectroscopic properties. For the spectroscopic measurements, we measure the energy transfer probability for fixed {J, κ, ∆} and simulation time τ sim , varying the frequency ν eff .

III. ENERGY TRANSFER DYNAMICS
VAET is especially well illustrated in the regime where the detuning ∆ is larger than the coupling J. Under these conditions, the energy transfer is suppressed without the assistance of the environment [ Fig. 2(a)]. Significant energy transfer occurs only at appropriate environmental frequencies corresponding to processes |SD, n → |DS, n ± k , k ≥ 1, where k vibrational quanta are either removed from or absorbed by the environment. The most pronounced energy transfer is observed at ν eff = ±2π × 4 kHz, where the occupation of the environment changes by 1. The unresolved peaks at smaller vibrational frequencies correspond to multiphonon processes where the energy gap is bridged by several vibrational quanta. The physical interpretation of ν eff < 0 is a sign change of the detuning ∆ and corresponds to the situation where the environment absorbs the excess energy.
With the environment prepared near the ground state, the vibrations can take up excess energy, but they cannot provide quanta in the transfer process. As a result, the spectral scans in Fig. 2(b), where the environment contains an average of 0.5 quanta, show a suppression for positive ν eff but exhibit strong peaks at negative ν eff . This asymmetry is expected only when the environment is prepared close to the ground state and together with the quantized nature of the signal, it is a signature of the quantum-mechanical nature of the environment.
To study these processes in more detail, we measure the time dynamics of the single-phonon resonances (ν eff = ±2π × 4 kHz). In the low-temperature setting, the time for maximum energy transfer to the acceptor state is around 1 ms, corresponding to a characteristic time scale of Jκ 2∆ given by a perturbative treatment (see Ref. [24]). The dynamics is faster at higher temperatures owing to stronger fluctuations from the environment. The previously mentioned asymmetry between the positive and negative frequency peaks is also present in the low-temperature time dynamics.
For sufficiently large detuning ∆, the multiphonon processes |SD, n → |DS, n + k can each be spectrally resolved up to a maximum integer value. However, when ∆ becomes comparable to J, these processes cannot be distinguished, as shown in Fig. 3(a). In this regime, the energy sites are partially hybridized because of the intersite coupling, and significant energy transfer occurs even without a coupling to the environment. Here, direct intersite coupling, single-phonon, and multiphonon processes all contribute simultaneously to the transfer dynamics. This parameter regime, where J ∼ κ ∼ ∆, is most relevant to energy transfer dynamics in photosynthetic light-harvesting complexes [25]. We note that, exactly in this regime, approximate methods, which typically proceed by perturbation in one of these parameters, are not effective (see Ref. [24]).
When the environment is cooled near the ground state, we observe almost complete energy transfer between donor and acceptor states [see Fig. 3(b)]. This energy transfer is accelerated with increasing site-environment coupling κ. Similarly, we see that the perturbative treatment developed in Ref.
[24] breaks down at earlier times. We note that the non-Markovianity of the environment is particularly evident for the case of κ = 2π × 0.64 kHz, since the energy returns to the acceptor during the measurement time. Finally, in Fig. 3(c), we show time dynamics at higher temperatures [n = 12, with all other parameters similar to the lowest time trace in Fig. 3(b)]. The overall trend in population transfer is preserved; however, the coherent oscillations are damped at the higher temperature.

IV. CONCLUSIONS AND OUTLOOK
Realistic models of chemical and biological environments require extending the simple model above to incorporate larger numbers of sites and environmental modes. However, as the number of sites and vibrational modes that must be accounted for increases, these models quickly challenge analytical and numerical methods. This is particularly the case in regimes where a separation of time scales does not exist and perturbative approximations are invalid. Moreover, it may be necessary to include higher excitations of the environment. In particular, vibrationally assisted processes are significant when the vibrational modes are almost resonant with the electronic energy differences, typically of order 100-200 cm −1 in photosynthetic systems. At room temperature, the relevant vibrations are each excited with mean phonon numbers on the order of one to two quanta. Capturing the associated Boltzmann distribution via direct numerical simulation may require truncating each harmonic oscillator Hilbert space above five quanta -the computational resource equivalent of 2-3 qubits. Thus, an N -site model with two vibrational modes per site to capture at least some aspects of the local spectral density would require a Hilbert space size comparable to (1 + 2 × 2)N = 5N qubits. To model the local environment more accurately, a Hilbert space with substantially more dimensions would be required. Currently, it is impossible to study the general dynamics in a Hilbert space equivalent to 50 qubits [26]. Extending our experimental platform, it may be possible to encode the 10-site dynamics of our model in a 10-ion crystal, thereby outperforming bruteforce classical computation. More sophisticated numerical techniques for simulating vibrational wave-packet dynamics have been developed recently, e.g., the multilayer multiconfiguration time-dependent Hartree method (ML-MCTDH) [27]. While such methods substantially reduce the computational burden of computing nonequilibrium dynamics substantially, their computational cost is still exponential with the number of vibrational modes tracked; hence, they eventually scale badly.
Even without outperforming classical resources, our platform can be extended to study more qualitatively interesting physics. For example, by applying phase modulation to the site-environment coupling beams, a broadband bath can be implemented in the simulation. Damping can be implemented by adding sideband cooling of particular vibrational modes using auxiliary ions. In addition, moving to a three-site model would allow the investigation of phenomena such as quantum ratcheting [28]. Finally, we note that one could also study steady-state dynamics of VAET by continuously exciting some of the ions on the |S 1/2 ↔ |D 5/2 transition while providing a sink to other ions using light on the |D 5/2 ↔ |P 3/2 transition. Thus, it may be possible to simulate realistic models of energy transfer processes in light harvesting processes or similar transport phenomena such as present in organic electronic devices including solar cells [29][30][31][32] or as discussed in the context of olfaction [33] and neuroreceptor activation [34].
In conclusion, we have implemented an analog quantum simulation of vibrationally assisted energy transfer using trapped ions. We further demonstrated tuning of the simulator from perturbative to nonperturbative regimes (see Ref. [24] for a more detailed discussion). The latter case is particularly interesting, as VAET dynamics with larger structured environments in this regime becomes inaccessible to numerical treatment on current high-performance computers. We expect that our platform will be capable of simulating complex models, including larger structured environments, with various experimental advances.
Recently, we became aware of related work carried out at ETH Zürich [35].

B. Perturbation in κ
The most common approximation in open quantum system models is that the system-environment interaction κ is small. In this regime, master equations, perturbative in κ, are often used to describe the system dynamics [4,36]. However, we cannot resort to master equations in this case because the vibrational mode is not actively damped in the VAET model described by Eq. 1 in the main text. The absence of damping renders the Markovian approximation of a memoryless environment, and therefore master equation approaches, invalid.
An approximation of the dynamics can still be obtained in the small κ regime by performing a perturbation expansion of the interaction picture propagator While we can calculate this quantity to any order in the perturbation expansion, the increasing computational expense renders this approach impractical. Fig. 4 shows various perturbative approximations to the dynamics shown in Fig. 3(c) of the main text. As can be seen from this figure, this perturbative approximation always breaks down at some timescale, and this breakdown happens earlier as κ increases. Our experimental simulations extend into a regime not described by this perturbative treatment. A more sophisticated perturbation expansion in the system-environment interaction is possible by first performing a polaron transformation of the VAET Hamiltonian (the effective system-environment coupling is reduced in the polaron frame) [37]. However, such expansions will also fail to desribe sufficiently long-time dynamics.

C. Large ∆ regime
Another regime where one can develop a simpler picture of the VAET dynamics is when ∆ dominates over the other parameters (J, κ). To see this, we write out the explicit form of the the interaction picture representation of σ z :σ z (t) = f x (t)σ x + f y (t)σ y + f z (t)σ z , with where Ω ≡ √ ∆ 2 + J 2 and in the second equality of each term we have simply rewritten the expression in terms of (J/∆).
Then the propagator in the interaction frame can be explicitly written as where The key observation now is that amplitude of vibrationassisted transitions from donor to acceptor will be proportional to powers of F ± x and F ± y . We will now examine one of these terms, F − x (t), but the following arguments apply to all of these terms that quantify the magnitude of vibration-assisted transitions. Expanding F − x (t), There have been no approximations till this point. Now consider the case where ∆ dominates over J and κ (e.g., Fig.2 in the main text). Then the prefactor κ(J/∆) 2(1+(J/∆) 2 ) ≈ κ(J/∆) 2 is small and due to the oscillatory nature of the integral, |F − x (t)| is small for t 1/Ω, 1/ν eff , except for when ν eff ≈ ±Ω. This explains the resolved peaks in Fig. 2 of the main text; i.e., for large ∆, vibration-assisted transitions are only significant when ν eff ≈ ±nΩ since all powers of F − x (t) (and the other coefficients) are generated by the exponential in Eq. 8.
In this regime, due to this strong resonance condition, one can develop simple few-level models that approximate VAET dynamics.

D. The non-perturbative regime of VAET
Finally, consider the deeply non-perturbative regime of the VAET model where J ∼ ∆ ∼ κ. This is the situation depicted in most of the panels in Fig. 3(a) and 3(b) of the main text. In this case, the perturbative expansion of the propagator in κ is not valid and the approximations based on small J/∆ do not hold. Consequently, the prefactor κ(J/∆) 2(1+(J/∆) 2 ) of the expression for F − x (t) in Eq. 10 is not small. The oscillating integrals still decay for values of ν eff far away from ±Ω, but the larger prefac-tor means that the vibrational frequency must be further away from resonance before |F − x (t)| becomes negligible. This explains the broader peaks in the spectrum shown in Fig. 3(a) of the main text.
For these parameters, we must simulate the full model for accurate predictions at timescales longer than 1/Ω and 1/ν eff .
We note that this regime, where the resonance condition is violated, resembles the breakdown of the rotatingwave-approximation in the Rabi model (which makes the simplification to the Jaynes-Cummings model invalid) [38].