Simulating polaron biophysics with Rydberg atoms

Transport of excitations along proteins can be formulated in a quantum physics context, based on the periodicity and vibrational modes of the structures. Numerically exact solutions of the corresponding equations are very challenging to obtain on classical computers. Approximate solutions based on the Davydov ansatz have demonstrated the possibility of stabilized solitonic excitations along the protein, however, experimentally these solutions have never been directly observed. Here we propose an alternative study of biophysical transport phenomena based on a quantum simulator composed of a chain of ultracold dressed Rydberg atoms, which allows for a direct observation of the Davydov phenomena. We show that there is an experimentally accessible range of parameters where the system directly mimics the Davydov equations and their solutions. Moreover, we show that such a quantum simulator has access to the regime in between the small and large polaron regimes, which cannot be described perturbatively.


The simulator
We mimic the behavior of the above-mentioned bio-molecules with a system of ultracold atoms confined in a very deep one-dimensional optical lattice potential V x V xR ( ) sin (2 / ) where each lattice site is occupied exactly by one atom. We assume that the spatial dynamics of atoms is not completely frozen, i.e., atoms may oscillate in vicinities of local minima with frequency Besides spatial motion, each atom may exhibit changes of its internal state due to the long-range interactions between neighboring atoms via Rydberg dressing mechanism 41,[45][46][47][48] . Following the work by Wüster et al. 49 coupling of the internal state of atoms with their spatial motion can be realized by the off-resonant coupling of two different but degenerated internal Zeeman levels in the different hyperfine states of the ground-state manifold of the atoms. The large hyperfine splitting will permit selective addressing of the levels |g〉 and |g′〉 to two precisely selected, highly excited Rydberg states |nS〉 or |nP〉 (via two-and single-photon transition, respectively) with principal quantum number n and angular momentum equal to 0 or ħ, respectively. A perturbation analysis shows that this coupling results in a quite small admixture of a Rydberg state to the atomic ground states and, as a consequence, atom can be found in one of the two dressed states 49 : where amplitudes α l = Θ l /2Δ l (l ∈ {s, p}) are determined by a total Rabi frequency of a driving field Θ l and a total laser detuning Δ l . In this basis of dressed states the dipole-dipole interaction between neighboring atoms is C R / sp 3 3 (R is the spatial distance between the atoms), besides additional contribution to the energy gap between local states |0 i 〉 and |1 i 〉, may induce transitions (excitation hoppings) between internal states of neighboring In consequence, the excitation |1〉 can be effectively transported across the lattice. This effect is driven by the following Hamiltonian of internal motion of all atoms: where an annihilation operator of an excitation a î can be viewed as a local transition operator |1 i 〉〈0 i | between dressed Rydberg states. The spatial dependent parameters W i and J i,i+1 are related to the dipole-dipole forces induced by Rydberg dressing and are given by 49 : where κ Δ =  R C R ( ) ( / )/ sp 3 3 and Δ = Δ s + Δ p . In a static situation, when all atom positions are frozen, the energies W i and J i+1,i are site-independent with values controlled by dipole-dipole interactions between neighboring atoms at fixed lattice spacing R 0 . However, due to the vibrational motion of atoms, these parameters are position dependent and they couple internal states of atoms with their motional degrees of freedom. In the lowest order of approximation they can be written as where g W and g J are the appropriate Taylor expansion coefficients of (4) around R 0 . Moreover, since the vibrational motion is quantized, the parameters have an operator character when acting in the subspace of spatial motion of atoms. By inserting expanded W i and J i+1,i to the Hamiltonian (3) one obtains HHSH Hamiltonianˆˆˆˆˆˆˆ † † † † where the first term describes the excitation dynamics on the lattice, the second term describes the excitationvibration coupling via the on-site energy, and the last term represents the off-site coupling via the excitation-vibration coupling through the hopping energy. Our implementation can be regarded as a dedicated quantum simulator to study excitation stabilization by vibrations related to the α-helix protein.
For the Rydberg parameters that we consider, the second order terms in the Taylor expansion are significantly smaller than the first-order corrections, which justifies a linear approximation.
For the moment we comment on a special case of Hamiltonian (5) with zero off-site coupling g J = 0, i.e. Holstein model. In this model the excitation is dressed by phonons forming a polaron quasiparticle. Two limiting cases have exact solutions: (i) For zero excitation-vibration coupling g W = 0 eigenstates are well described by product states in Fourier space , etc. These states construct a good basis for perturbation expansions in the g W /J 0 parameter for the weak coupling limit; (ii) a second limiting case, called the small polaron, corresponds to a zero hopping energy term J 0 = 0 in which the Hamiltonian can be diagonalized in the dressed-polaron picture, i.e.
form a good basis for a perturbative description in the small g W /ħω 0 parameter. A perturbative calculation shows that the excitation-vibration coupling is given by a dressed hopping amplitude where for ω 0 large enough, the hopping amplitude vanishes. For the HSSH model in the regime of parameters giving rise to a Davydov soliton, i.e. where all parameters are of the same order, a good small parameter for a perturbative description is lacking, and therefore a variational approach is preferred.
The Rydberg dressing is responsible for the coupling between vibrational degrees of freedom of neighboring atoms. The resulting dressed soft-core interaction is proportional to 41 , gives also rise to an additional energy shift. However, this shift is negligible compared to ħω 0 , and therefore we omit it. In the following, all energies are expressed in units of J 0 , and time is measured in units of

Dynamical properties of the system
An important question related to the dynamics of the HSSH Hamiltonian is whether the lattice vibrations are able to stabilize the excitation that is initially localized on a specific site K ˆ † Ψ | 〉 = a K 0 |vac〉, or slightly delocalized on two neighboring sites ~⟩ˆ † † where |vac〉 is the vacuum state of the system fulfilling the condition a î |vac〉 = b i |vac〉 = 0 for any i. For certain parameters, a system prepared in these initial states evolves in such a way that the excitation does not spread across the protein. This is attributed to a specific ratio of the interactions of excitation and vibrational degrees of freedom, giving rise to a soliton. This spreading or non-spreading behavior can be extracted from information encoded in the time-dependent density profile where the state of the system at given time t can be formally written as where |Ψ ini 〉 is one of the considered initial states. Temporal spreading of the excitation is then given by an effective width of the spatial density profile This quantity takes the value 1/N for an excitation localized at exactly one lattice site and 1 when fully delocalized. In principle, by analyzing the time-dependence of σ(t) one can easily determine whether the excitation remains localized or whether it spreads across the system. Numerically exact solutions of the evolution problem are very challenging due to the strong non-linear quantum-mechanical coupling between excitation and vibrational degrees of freedom. Therefore generally the evolution of the system cannot be found exactly and some approximation methods have to be adopted.

The Davydov approach.
We discuss here the two-step Davydov approach 5 , which results in a semiclassical description of the system. In the first step one assumes that the state of the system | Ψ(t)〉 can be well approximated by the product of two independent states |Ψ(t)〉 and |φ(t)〉 for excitation and vibrational degrees of freedom, respectively, |Ψ(t)〉 = |ψ(t)〉|φ(t)〉. Since the system is initially prepared in the state with precisely one excitation and the number of excitations is conserved, the state |Ψ(t)〉 can be decomposed in the single-particle subspace, t ta ( ) ( ) i i i ψ ψ | 〉 = ∑ˆ † |vac〉, where time-dependent functions ψ i (t) play the role of probability amplitudes for finding an excitation at site i. Consequently ρ i (t) = |ψ i (t)| 2 . The second step relays on a semi-classical treatment of the vibrational degrees of the system. In analogy to other quantum field theories, we assume that the state |φ(t)〉 has classical features, i.e., it can be well approximated by the product of independent coherent states: where amplitudes u i (t) and p i (t) are expectation values of appropriate operators in the state |φ(t)〉. Within these approximations, we calculate the expectation value of the many-body Hamiltonian on the system state and approximate the resulting equation of motion by the classical Hamilton equations 50 , we obtain set of coupled differential equations of the form: Phase diagram. We perform a semi-classical evolution of the system governed by eqn. (7), which allows us to observe spreading or non-spreading evolution of an effective width of the spatial density profile as a function of the parameters {ω 0 , g W , g J }. The results can be visualized by the phase diagrams presented in Figs 1 and 2. These diagrams are obtained by plotting the maximal value of the σ(t) reached during the evolution up to maximal time T max = 10. In order to avoid interference effects during the evolution caused by the boundaries, all calculations are performed with a sufficiently large lattice of N = 50 lattice sites and with periodic boundary conditions. We checked that the numerical results are insensitive to enlarging N on the time-scales of study, and therefore the obtained results are also valid for infinitely large systems. Moreover, the chosen lattice size is similar to current experimental efforts in this direction. First, we focus on the case of a completely localized initial state |Ψ 0 〉 for g J = 0 (Fig. 1). We qualitatively indicate five different regions on the phase diagram (left panel): (I) and (II) where the excitation is dressed by a cloud of vibrations and the excitation does spread; (III) where due to an exponential reduction of the hopping amplitude the excitation is localized in its initial position 59 ; (IV) where the sum of vibration energy and exciton-vibration coupling is larger than the hopping energy, giving rise to Davydov-like soliton behavior; (V) where ω g W 0  corresponding to the Discrete Breathers-like behavior 60,61 . Distinct behavior of the system is also visible in these selected areas in the time evolution of σ(t) (right panel of Fig. 1). This picture can be generalized to non-vanishing coupling g J , which we investigate for the the second initial state 0 |Ψ 〉 ∼ (Fig. 2). As can be seen, a slight delocalization of the initial state together with non-local coupling g J dramatically enhance the non-spreading behavior of the wave packet. It is a direct consequence of the non-local terms in (7).
Numerically exact approach. The results obtained in the framework of the semi-classical Davydov approach can be supported by numerically exact dynamics governed by the many-body Hamiltonian (5). In this   (5), therefore an exact evolution is obtained only in the limit where all Fock states are taken into account. In practice, for numerical purposes, we assume that the total number of excitations cannot be larger than some well defined cut-off M. Then the results are treated as exact if increasing M does not change the outcome noticeably 62,63 . Therefore, for a given M, one can perform calculations only for a small range of parameters for which creations of vibrations is limited. It is worth noticing that numerical complexity grows exponentially with the cut-off M. For our parameters, N = 50 and M = 3 the size of the corresponding Hilbert space exceeds 1.1 million. Other approaches based on quazi-exact dynamics are presented in 64,65 . As already expected from Fig. (2), the numerically exact approach confirms that the on-site coupling (g W ≠ 0) plays a dominant role in the excitation stabilization process. Therefore, without loss of generality, we consider now the minimal-coupling scenario with g J = 0 In Fig. 3 we show the time evolution of an initially localized excitation for ω 0 = 3 and for three different values of the local coupling parameter g W = {0.1, 0.75, 1.5}. It is clearly visible that for larger g W the wave packet of the excitation becomes more stable and spreads less. This effect is directly reflected in the number of vibrational modes created, which can be seen in the bottom row of Fig. 3. One can observe that increasing fluctuations of the total vibrations in the system stabilize excitation. Since we reached the limits of our computational method with this size of the Hilbert space, we cannot increase the coupling parameter further. From Fig. 3 it can be seen that the total number of vibrations for g W = 2 is close to the limiting cut-off. At the same time, however, this is a strong argument for employing a quantum simulator, such as proposed in this letter, to validate the predictions of the semi-classical approach.

Experimental parameters
The numerical predictions for the model described by the Hamiltonian (5) are quite general. For a quantum simulator we consider 87 Rb atoms confined in an optical lattice determined by lattice spacing R 0 = 1 μm and V 0 = 100 E R (recoil energy π = E m R 2 / R 2 2 0 2  ) 45 , i.e., the local trap frequency is equal to 6.2 kHz. We assume Rydberg states with principal quantum number n = 50 for which C sp

Summary
We show that a system of ultracold Rydberg atoms confined in a one-dimensional optical lattice may serve as a dedicated quantum simulator for excitation-vibration dynamics, which is a subclass of polaron dynamics. Since effective parameters of the resulting model can be easily tuned, the system can be used to mimic transport of excitation in biologically active proteins and to perform full quantum mechanical tests of the semi-classical predictions. The proposed scheme may serve as a platform to investigate the HSSH bi-and many-polaron system 69,70 . In particular, the character of the bi-polaron interactions can be tuned from repulsive to attractive by the experimental control parameters g W and g J . Finally, we note that also disorder effects in the HSSH Hamiltonian 71 can be studied by introducing incommensurate optical lattices.