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.

Introduction

The understanding of molecular structures provides life sciences with tools to explain complex cell-biology phenomena. Biological complexity of mesoscopic objects along with the quantum behavior of their basic elements leads to interesting unsolved questions awaiting comprehensive answers[1]. The interdisciplinary field of “quantum biology” is the natural area for combining quantum physical methods and tools to investigate, model and simulate biological systems on a mesoscopic level[2-4].

Many biological processes are powered by the energy released from the hydrolysis of adenosine triphosphate (ATP). On a physical level, it can be viewed as a vibrational bound state of the ATP molecule to a protein with energy equal to 0.49 eV. In the 1970’s Davydov proposed a mechanism for the localization and transport of the associated vibrational energy in the α-helix region of a protein by means of a so-called Davydov soliton[5-8]. Although this model has been used for a theoretical description of experimentally observed unconventional absorption bands in proteins[9,10], direct experimental evidence for existence of this soliton is still missing.

The Davydov soliton is a subclass of richer polaron phenomena, i.e., excitations mediated by phonons originally introduced by Landau[11]. Polarons have been broadly studied theoretically as well as experimentally in a condensed matter context[12] and recently in different areas of ultracold physics: ultracold ions[13-19], polar molecules[20-24], ultracold Rydberg gases[25-30], and strongly-interacting ultracold Bose and Fermi gases[31-38]. In this paper we show that the Davydov soliton can be created and observed in a suitably prepared system of ultracold atoms, confined in an optical lattice and off-resonantly coupled to a Rydberg state[39-42]. Such a system can be regarded as a dedicated quantum simulator within the broader class of the Holstein-Su-Schrieffer-Heeger (HSSH) Hamiltonian[43,44]. First, we investigate dynamical properties of the system within the semi-classical Davydov approach assuming an infinite number of phonons. Second, we use an exact evolution approach to study the dynamics of the system in the non-classical few-phonons regime. In a particular, experimentally accessible parameter regime, both approaches confirm the existence of soliton solutions. This also indicates that with this Rydberg quantum simulator it is possible to study the regime in between the so-called small and large polaron regimes of the HSSH model, which cannot be described by perturbative theoretical methods.

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 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 . This motion is however quantized and therefore it is driven by a simple harmonic oscillator-like Hamiltonian:where and are the position and momenta operators related to i-th atom, while operator annihilates vibrational excitation of the i-th atom. Local motion defines a natural scale of length, .

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-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 α = Θ/2Δ (l ∈ {s, p}) are determined by a total Rabi frequency of a driving field Θ and a total laser detuning Δ. In this basis of dressed states the dipole-dipole interaction between neighboring atoms is (R is the spatial distance between the atoms), besides additional contribution to the energy gap between local states |0〉 and |1〉, may induce transitions (excitation hoppings) between internal states of neighboring atoms |0〉|1〉 ↔ |1〉|0〉. 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 can be viewed as a local transition operator |1〉〈0| between dressed Rydberg states. The spatial dependent parameters W and J are related to the dipole-dipole forces induced by Rydberg dressing and are given by[49]:where and Δ = Δ + Δ. In a static situation, when all atom positions are frozen, the energies W and J are site-independent with values controlled by dipole-dipole interactions between neighboring atoms at fixed lattice spacing R0. 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 W = W0 + g(u + u), J = −J0 + g(u − u), where gW and gJ are the appropriate Taylor expansion coefficients of (4) around R0. 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 and J to the Hamiltonian (3) one obtains HHSH Hamiltonianwhere the first term describes the excitation dynamics on the lattice, the second term describes the excitation-vibration 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 = 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 = 0 eigenstates are well described by product states in Fourier space i.e. , etc. These states construct a good basis for perturbation expansions in the gW/J0 parameter for the weak coupling limit; (ii) a second limiting case, called the small polaron, corresponds to a zero hopping energy term J0 = 0 in which the Hamiltonian can be diagonalized in the dressed-polaron picture, i.e.where . The dressed-polaron eigenstates with eigenenergies form a good basis for a perturbative description in the small g/ħω0 parameter. A perturbative calculation shows that the excitation-vibration coupling is given by a dressed hopping amplitude [12], 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 , with R the Rydberg blockade radius[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 J0, and time is measured in units of ħ/J0, i.e. we set J0 = ħ = m = 1.v

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 |vac〉, or slightly delocalized on two neighboring sites |vac〉, where |vac〉 is the vacuum state of the system fulfilling the condition |vac〉 = |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, |vac〉, where time-dependent functions ψ(t) play the role of probability amplitudes for finding an excitation at site i. Consequently ρ(t) = |ψ(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: |vac〉, where amplitudes u(t) and p(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:

These are Davydov type equations[5-8], which describe the dynamics of an excitation ψi coupled to a gradient of a classical phonon field u forming an effective self-trapping potential. An alternative derivation provided by Kerr[51] is based on the Heisenberg equations of phonon position and momentum operators. A complementary approach to the above Davydov equations, which are based on the Lagrangian variational principle, is one based on the Dirac-Frenkel-McLachan (DFM) variational principle. This approach is commonly used in quantum molecular dynamics[52-58], in which equations of motion for the variational parameters are obtained from the minimization of , where δϕ denotes possible variations of ϕ with respect to the variational parameters.

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, g}. 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 Tmax = 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 = 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 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, which we investigate for the the second initial state (Fig. 2). As can be seen, a slight delocalization of the initial state together with non-local coupling gJ dramatically enhance the non-spreading behavior of the wave packet. It is a direct consequence of the non-local terms in (7).

Figure 1

Left panel: Maximal value of an effective width of the spatial density profile max[σ(t)] as a function of ω0 and g for vanishing coupling g = 0. A sharp crossover between non-spreading excitations (blue) and spreading excitation (dark red) is clearly visible. Different regions of the phase-diagram (bordered with white lines) correspond to a distinct nature of the exciton-vibration dynamics. Right panel: Evolution of the excitation width σ(t) for different points on the phase diagram (marked as white squares on the left panel).

Figure 2

Maximal value of the wave packet width max[σ(t)] for different initial states |Ψ0〉 and (top and bottom row, respectively) and different non-local interactions g = {0, 3, 5} (appropriate columns from left to right). Note that strong enhancement of the non-spreading behavior takes place for stronger g and for a smeared out initial state.

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 approach we represent Hamiltonian H as a matrix in the Fock basis spanned by many-body states |vac〉, i.e., states with an excitation located exactly at site i and with selected vibrational states mi for all sites. An arbitrary state of the system can be expressed as an appropriate superposition of the basis states. Since the operator of a total number of vibrations in the system does not commute with the Hamiltonian (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 ≠ 0) plays a dominant role in the excitation stabilization process. Therefore, without loss of generality, we consider now the minimal-coupling scenario with g = 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 = {0.1, 0.75, 1.5}. It is clearly visible that for larger gW 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 = 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.

Figure 3

Exact evolution of the effective width of the spatial density profile governed by the Hamiltonian (5) for the initial state |Ψ0〉. The bottom row shows the total number of vibrations created in the system during the evolution. Consecutive columns correspond to different local couplings g = {0.1, 0.75, 1.5}. All calculations performed for g = 0 and ω0 = 3. Note that for stronger interactions evident stabilization of the excitation density profile, along with increasing number of created vibrations, is observed.

Experimental parameters

The numerical predictions for the model described by the Hamiltonian (5) are quite general. For a quantum simulator we consider 87Rb atoms confined in an optical lattice determined by lattice spacing R0 = 1 μm and V0 = 100 E (recoil energy )[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  = 3.224 GHz μm3 [66]. We choose the dressing parameters as α = 0.015 and Δ/2 = Δ = Δ = 2.5 GHz. With these values, the system mimics the Hamiltonian (5) with dimensionless parameters ω0 = 4.7, g = 5.6, and g = 5.6. These parameters can be easily tuned since they strongly depend on the lattice spacing R0 and on the set of laser detunings. In this way, a large and interesting area of the phase diagrams presented in Fig. 2 can be covered. The estimated lifetime of Rydberg atoms excited to states with n = 50 is τ = 65 μs and τ = 86 μs[67]. The effective lifetime of a Rydberg dressed state is scaled by a factor α−2 and is sufficiently long to observe the non-spreading excitation behavior. It is worth noting that also other experimental realizations, based for example on Rydberg microtrap arrays[68], can be considered as proper candidates for simulating this system.

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 and g. Finally, we note that also disorder effects in the HSSH Hamiltonian[71] can be studied by introducing incommensurate optical lattices.