Emulating Molecular Orbitals and Electronic Dynamics with Ultracold Atoms

In recent years, ultracold atoms in optical lattices have proven their great value as quantum simulators for studying strongly correlated phases and complex phenomena in solid-state systems. Here we reveal their potential as quantum simulators for molecular physics and propose a technique to image the three-dimensional molecular orbitals with high resolution. The outstanding tunability of ultracold atoms in terms of potential and interaction offer fully adjustable model systems for gaining deep insight into the electronic structure of molecules. We study the orbitals of an artificial benzene molecule and discuss the effect of tunable interactions in its conjugated pi electron system with special regard to localization and spin order. The dynamical time scales of ultracold atom simulators are on the order of milliseconds, which allows for the time-resolved monitoring of a broad range of dynamical processes. As an example, we compute the hole dynamics in the conjugated pi system of the artificial benzene molecule.

In the recent years, ultracold atoms in optical lattices have proven their great value as quantum simulators for studying strongly-correlated phases and complex phenomena in solid-state systems. Here we reveal their potential as quantum simulators for molecular physics and propose a technique to image the three-dimensional molecular orbitals with high resolution. The outstanding tunability of ultracold atoms in terms of potential and interaction offer fully-adjustable model systems for gaining deep insight into the electronic structure of molecules. We study the orbitals of an artificial benzene molecule and discuss the effect of tunable interactions in its conjugated π electron system with special regard to localization and spin order. The dynamical timescale of ultracold atom simulators are on the order milliseconds which allow for the time-resolved monitoring of a broad range of dynamical processes. As an example, we compute the hole dynamics in the conjugated π system of the artificial benzene molecule.
The structure of molecules is usually determined by x-ray or electron diffraction. Current advances with femtosecond pulses allow for the time-resolved observation of the atomic positions 1 . In general, orbital wave functions are much harder to image. For simple molecules, the reconstruction of the highest occupied molecular orbitals (HOMO) has been recently achieved by electron momentum spectroscopy 2 , by using laser-induced electron diffraction 3 , and by higher-order harmonic generation 4,5 . For hydrogen atoms the nodal lines of Stark states were resolved via photoionization microscopy 6 . State-of-the-art scanning probe microscopy allows resolving the HOMO as well as identifying the chemical bonds 7 for benzene-based compounds attached to surfaces.
Experimental access to the electronic structure of molecules and their dynamics is essential because even for relatively small molecules the full many-particle problem is not computable using classical computers. This has brought forward the idea of quantum computation for molecules 8,9 . The quantum computation of a hydrogen molecule in a minimal basis has already succeeded 10,11 , but estimates for more complex molecules are not promising. While the number of required qubits appears manageable, the estimated number of quantum gate operations (10 18 for Fe 2 S 2 ) is many orders of magnitude larger than currently available 12 . In the last decade, ultracold atoms in optical lattices have established as quantum simulators for condensed matter systems 13,14 . In addition, ultracold atoms were proposed as simulators for different systems such as neutron stars 15 , black holes 16,17 , quarks 18 , or atoms 19 . Recent experimental developments include prospects such as the single-lattice-site imaging and addressing 20,21 as well as the deterministic preparation and detection of few-atom samples 22,23 .
Here, we propose an experimental scheme for the quantum simulation of non-trivial molecules using ultracold atoms based on existing experimental techniques. In this setting, the ultracold atoms mimic the electrons in a molecule, whereas the optical trapping potential takes the role of the nuclei. We show that ultracold atoms can serve as a tunable model system allowing the investigation of open questions in molec-ular physics. As a concrete example, we focus on a model system for benzene and compute its molecular orbitals for vanishing interaction. For the non-solvable interacting problem, ultracold atoms can serve as a quantum simulator for static and dynamical electronic properties in molecules. We demonstrate numerically that the conjugated π electron system shows strong-correlation effects such as localization and spin order, even when neglecting the interaction with inner particles. On the basis of this subsystem, we also reveal that ultracold atom simulators promise unique insight into electronic femtosecond dynamics. We show that momentummapping in combination with phase retrieval allows imaging molecular orbitals with 1-2 orders of magnitude better than the intramolecular distances. We explain how to use the outstanding tunability of interaction and potential for studying electronic interactions and the dynamics in artificial molecules.

Creating artificial molecules
As an example, we discuss in the following how to create an artificial benzene molecule and compute the orbital wave functions. In real benzene, the molecular structure is formed by a ring of six carbon atoms and six hydrogen atoms (Fig. 1c). Thereby, the molecular symmetry is essential for the formation of molecular orbitals, where benzene belongs to the point group D 6h . In our case, the idea is to simulate the electrons in a molecular structure via ultracold atoms in a tailored optical dipole potential which mimics the Coulomb interaction with the nuclei. The optical potential is imposed by the AC Stark shift of a laser field and can be created for the example of benzene by superposing a hexagonal optical lattice V 0 V hc (x, y) 24 with a tight dipole trap V dip (x, y) as shown in Fig. 1a (see Methods). A similar potential can, e.g., also be created by adding a triangular superlattice to the honeycomb lattice (a "benzene lattice" 25 ) or by using a spatial light modulator, which renders the possibility of almost arbitrary two-dimensional potentials 26,27 . In the orthogonal direction a one-dimensional lattice V 1D with depth V 0z or a light sheet is applied. A controlled number of fermionic atoms (e.g. 6 Li atoms with spin ±1/2) can be loaded into this optical   The molecular orbitals of the artificial molecule are computed using the plane-wave expansion method for vanishing interaction (see Methods for technical details and parameters). The resulting single-particle orbitals are shown in Fig. 1b. The six lowest s orbitals represent an energetically separated shell within the inner "C" ring, whereas the higher orbitals (7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17)(18)(25)(26)(27)(28)(29)(30) are sp 2 -hybridized. The orbitals (7)(8)(9)(10)(11) and the antibonding orbital (13) form p orbitals pointing along the C ring. The outwards pointing p orbitals (12) and (14-18) lower their energy by maximizing the overlap to the s orbitals on the H atoms. The hybridized orbitals (25-30) show a nodal plane between H and C ring and thus have a higher energy. Energetically in between lie the delocalized p z orbitals (19-24) with a nodal plane at z = 0. Benzene is one of the prime examples for this type of delocalized conjugation of a π electron system stabilizing the planar ring structure. This property is commonly referred to as aromaticity. Accounting for in total 42 electrons, in benzene three of the six p z orbitals are occupied. We can tune the energetic order of the orbitals, by tuning the parameters of the potential, e.g. strength and width of the dipole trap as well as the lattice depths V 0 and V 0z (see Fig. 3a).

Imaging of molecular wave functions
Quantum gas microscopes have demonstrated optical imaging with single-site resolution in an optical lattice 20,21 . However, the imaging of molecular wave functions of artificially created molecules requires an optical resolution several times better than the optical lattice spacing. The idea is to overcome the limited spatial resolution by imaging the particles in momentum space. In ultracold atom experiments, this can be achieved by turning off the optical potential, letting the particles freely expand and imaging the particle density ρ TOF (r) after a certain time of flight. In the language of molecular physics, this would be equivalent to a sudden removal of all nuclei. During the expansion, the interactions quickly become negligible, because there are only few particles with high momenta. This allows for a free expansion in two dimensions with an unaltered lattice in z direction avoiding problems with a limited depth of focus of the imaging system. In addition, the interaction among the particles can be tuned to zero via a Feshbach resonance during the expansion. For suitably long expansion times t (see Supplementary Information 25 ), the momentum density can be expressed as ρ p (q) ∝ ρ TOF (pt/m) with momenta p = hq/a, a the lattice constant of the honeycomb lattice, and h the Planck constant.
In the example shown in Fig. 2b, we assume that we can image the momentum density after time-of-flight with P 2 = 49×49 discrete momenta using a so-called pinning lattice for imaging (see Methods). Preparing initially one particle (or two with opposite spin without interaction) only the lowest orbital n = 1 is occupied and per experimental sequence only 1-2 momenta are recorded. To show the feasibility of the proposed imaging method, we use a random number generator with an accumulated number of 1000 momenta. Note that in principle satisfying results with fewer observed momenta are possible ( Supplementary Information 25 , Fig. S4). Figure 2b (first row) shows that the statics is clearly sufficient to map out the momentum distribution of the orbitals. The phase of the distribution can be retrieved by identifying the nodal lines of the momentum density (ρ p (q) = 0) and mapping it piece- x/a n = 1 n = 1 n = 12 , integrated momentum wave function and density in z direction. In combination with the reconstructed wave function in x-y coordinates (see (b)), this allows to retrieve the three-dimensional wave function for the lowest pz orbital (analogous for (b) and (c) using the lowest Wannier orbital in z direction).
wise to the momentum wave function with p(q) = ± ρ p (q) (second row). Since the potential is symmetric, p(q) can be assumed to be real. The components of p(q) are the Fourier coefficients of the wave function and therefore allow the direct transformation to orbital wave function in real space via with plane waves lattice constant a and the number of grid points S per reciprocal lattice vector. The resulting density is plotted in the third row of Fig. 2b. The figure demonstrates that the effective resolution of the imaging method is very high. In principle, the resolution is given by the diameter in momentum space. However, since the momentum wave function drops exponentially to zero for higher values of q, we can assume p(q) ≈ 0 for |q| 3.5. Since the diameter in momentum space is not limiting, the spatial resolution in the experiment is mainly influenced by the statistics of the momentum density (Supplementary Information 25 , Fig. S4). In contrary, the maximal diameter of the artificial molecule to be imaged is given by the resolution in momentum space and estimated as P/(2q max x −2q min x ) using the Nyquist theorem (3.5a for our example).
In analogy, higher single-particle orbital wave functions can be imaged as depicted in Fig. 2c. The preparation of higher orbitals is analogous to the preparation of higher bands in optical lattices via Bragg transfer 28,29 or via amplitude modulation 30 . Alternatively, increasing the number of particles in the artificial molecule automatically occupies higher orbitals. In this case, the recorded images correspond to the sum of the respective momentum densities (for vanishing interactions). Imaging the wave functions in three dimensions requires to measure the momentum distribution also in the z direction, where Fig. 2d depicts the reconstruction of the p z orbital. Here, a statistics of 100 particles is sufficient to resolve the orbital structure. In our case, the separable densities in x-y and z direction can be simply combined to the three-dimensional orbital wave functions shown in Fig. 2 (fourth row).
Quantum simulation of the many-body problem The single-particle energies of the orbitals are plotted in Fig. 3a and b, where the six p z orbitals are indicated in red. Due to the separability of the potential in z direction, these orbitals are analogous to the six lowest s orbitals but with a nodal plane for z = 0. The lowest and highest orbital are symmetric (A 2u ) and antisymmetric (B * 2g ) under rotation (compare Fig. 1b), whereas the second-third and fourth-fifth orbitals are doubly degenerate and labeled as E 1g and E * 2u . The lower three orbitals have a bonding and the higher three orbitals an antibonding character. While molecular potential curves are usually plotted against the nuclei distance R 0 , the distances are fixed in an optical lattice potential. Due to the same reason the repulsion of the nuclei, dominating the potential curves at R 0 → 0, is not present in optical potentials. However, even though the positions are fixed in our case, we can effectively change the intramolecular distances by tuning the depth V 0 of the lattice potential. A deeper lattice potential decreases the density between the sites and is analogue to a larger spacing R 0 between the nuclei coordinates.
A quantum simulator for molecular problems is needed for non-vanishing interactions, since solving the many-particle problem using the full configuration-interaction method is not feasible for complex molecules even accounting for the molecular symmetries. Unlike electrons, neutral atoms do not interact via Coulomb force, but either by a short-range contact potential or a long-range dipolar potential 31 . For ultracold atoms, the interaction strength between ultracold atoms can be tuned in a wide range via the magnetic field due to the existence of Feshbach resonances. As a concrete example, the s-wave scattering length a s for 6 Li atoms can be tuned between −0.2b and 0.4b, where b = λ/2 = 532 nm is the lattice spacing in z direction. Experimentally, we can probe the energies of the interacting problem (benzene has 42 electrons) by exciting a particle to a reference orbital using a microwave transition, which also changes the internal state. This avoids problems with post-interaction, assuming that the particles in the two different internal atomic states do not interact. Furthermore, the population in this internal atomic state serves as an observable measuring the transfer rate as a function of energy.
To demonstrate how interactions change the electronic configuration let us turn to the conjugated π system, i.e. the p z orbitals with delocalized particles, and neglect all interactions with the lower-lying orbitals. When restricting ourselves to this subsystem, the problem can be solved using the full configuration-interaction method. Experimentally, this subsystem is also accessible by populating only the six lowest s orbitals that are equivalent to the p z orbitals in the x and y direction. In the case of p z orbitals, the experimental system would incorporate the interaction with the energetically lowerlying particles constituting a much more complex situation. For the description of strongly correlated electron systems, it is often desirable to switch to Wannier orbitals (see Methods) that are localized at individual sites (potential minima). The Wannier orbital shown in Fig. 3d allows the formulation of a tight-binding Hubbard Hamiltonian taking into account the tunneling element t 1 and the on-site interaction U , where U scales linearly with the scattering length a s . The next-nearestneighbor tunneling t 2 is considerably smaller than t 1 but large enough to cause an asymmetry in the p z band (E 2 =E 0 in Fig. 3e). The fermionic many-particle Hamiltonian readŝ where the operatorĉ i,σ annihilates a particle on site i in spin state σ ∈ {↑, ↓} andn i,σ =ĉ † i,σĉ i,σ corresponds to the respective particle number. Here, the tunneling matrix elements t 1 (d = 1) and t 2 (d = 2) are included and [i] = i mod 6 incorporates the periodic boundary conditions.
For vanishing interaction (a s = 0), the energy spectrum ( Fig. 3c) represents all possible single-particle excitations of the half-filled band with delocalized particles (Fig. 3e). Increasing the interaction causes two drastic changes on the "electronic" structure. First, the particles undergo the transition from delocalized orbitals to the Mott insulator state, where the particles occupy the localized Wannier orbitals. At the same time, an interaction gap is formed separating the many-particle bands by the on-site interaction energy U . Second, the spins of the particles align in an alternating order on the Wannier orbitals (antiferromagnetic alignment), thereby gaining the second order energy 2t 2 1 /U per particle. The lowest Mott band contains possible spin excitations from this antiferromagnetic ground state. For the interaction strength a s = 0.02b, the ratio U/t 1 ≈ 4 is in the vicinity of the proposed values for aromatic hydrocarbons being still under debate 32,33 . The tunability of interaction strength and potential depth allows an independent control of the ratios U/t 1 and t 2 /t 1 (or the next-neighbor interaction for long-rage interactions 31 ).

Correlated electron dynamics
One of the outstanding advantages of ultracold atomic systems is that the dynamical timescale is on the order of millior microseconds. Therefore, artificial molecules would allow monitoring dynamical processes taking place on femtoor attoseconds timescales in real molecules. As an example, we compute the hole dynamics after the removal of a particle (see Fig. 4). Again, we restrict ourselves to the conjugated π system, for which the dynamics can be computed. The siteselective manipulation has been demonstrated for an optical lattice using a tightly focused addressing beam 26,34 . After an evolution time t, suddenly ramping up the optical lattice allows freezing the particle number distribution at the individual sites of the ring. Using fluorescence imaging, we can count the occupation number n i,↑ + n i,↓ on each site of the artificial molecule 20,21 . The spin state can be identified by selectively removing one of the spin states before imaging via resonant light 34 . In Fig. 4 the site-resolved dynamics is plotted after the removal of a spin up particle at site i = 3. Without interactions (Fig. 4a), the spin down particles are not influenced by the particle removal, whereas the initial hole in the spin up component oscillates between sites i = 3 and 6. This quantum revival occurs at the revival time T = h/E 0 (with E 0 = 0.13E R ). This corresponds to the energy difference between A 2u and E 1g orbitals. For 6 Li, the revival time corresponds to 0.27ms for a recoil energy E R = h×29.4kHz at a lattice wavelength λ = 1064 nm. Switching on the (repulsive) interaction among the particles (Fig. 4b), the dynamical behavior becomes very complex. For t = 0, the antiferromagnetic order is apparent in the spin down component causing the almost vanishing total density on site i = 3. The neighboring sites are highly occupied with the down component which fills the empty site on short timescales. This interferes strongly with the dynamics of the up component leading to complex density modulations. A priori, it is not clear if the spin order persists at the typical temperatures of ultracold experiment. Therefore, the dynamics shown in Fig. 4 incorporates a finte temperature of 0.05E R /k B with k B the Boltzmann, which corresponds to 70 nK for 6 Li, demonstrating that the quantum simulation of electronic dynamics is feasible with state-of-the-art temperatures.
In conclusion, ultracold atoms have a great potential as a quantum simulator for molecular physics, where solving the full many-particle problem is still out of reach using either classical or quantum computers 12 . This includes the formation of molecular orbitals, electronic interactions as well as the time-resolved monitoring of a broad range of dynamical processes. While we have discussed the creation of an artificial benzene molecule using state-of-the-art experimental techniques, the results can be easily generalized covering more complex molecules. The nuclei positions, kept fixed in the present study, can be controlled using slowly-varying time-dependent optical potentials (e.g. using a spatial light modulator or an accordion lattice 35 ) or coupled to phonons using other atomic species such as ions 36 . This allows simulating the electronic response to vibrational modes. In principle, one could also imagine to access regimes where the Born-Oppenheimer approximation breaks down.

Methods
Proposed experimental sequence and parameters. The optical potential for the artificial benzene molecule is generated by three laser beams in the x-y plane, which create a honeycomb hexagonal optical lattice V0Vhc(x, y), and an orthogonal one-dimensional lattice V0z cos 2 (πz/b) (Supplementary Information 25 ). Superposing this lattice potential with a tight dipole trap Vdip(x, y) = − dipV0e −2(x 2 +y 2 )/w 2 dip centered on a plaquette of the lattice as sketched in Fig. 1a causes a finite-size system with a hexagonal inner and outer ring, where V0 is the height of the potential (in units of the recoil energy ER = h 2 /2λ 2 m with the wavelength λ and the atomic mass m). For concreteness, we choose if not stated otherwise the lattice depths V0 = 11ER and V0z = 35ER, the relative strength dip = 10 of the dipole potential and the width wdip = 4a of the dipole trap. Here, a = 2λ/3 is the width of the unit cell in the honeycomb lattice, whereas b = λ/2 is the lattice spacing in z direction. The orbital order can be tuned by adjusting the lattice depth V0 and V0z (see Fig. 3a).
A controlled number of ultracold fermionic atoms can be prepared in the tight two-dimensional potential formed by a dipole trap and a single antinode of the orthogonal lattice. Atoms above a certain energy are expelled from the trap via a strong magnetic field gradient 22 . When adiabatically ramping up the honeycomb lattice, the particles fill up the lowest orbitals of the artificial molecule. To image the momentum distribution a time-of-flight measurement is performed, where the atoms expand for a certain time. Subsequently, a deep so-called pinning lattice is ramped up that holds the atoms in place while they are illuminated by an optical molasses 20,21 . The pinning lattice can be a square lattice of the same wavelength λ containing P 2 central lattice sites. Using a high-resolution imaging system the individual sites of the pinning lattice can be resolved 20,21 . In addition to the observed momenta, also the number of particles is retrieved.
For observing the momentum density in z direction, a threedimensional expansion can be performed. In this case, the imaging (from x direction) needs a large depth of field and is therefore restricted to a larger optical resolution. However, for dilute filling of the optical lattice, the atomic positions can still be identified even if the optical resolution is several times the lattice spacing 37 (see Supplementary Information 25 ). Computation of orbitals. We compute the molecular orbitals by applying the plane-wave expansion method for non-periodic systems. By decomposing the two-dimensional potential in Fourier components V (x, y) = q vq |q , the Hamiltonian matrix can be calculated in the plane wave basis (2). Here, S is the number of unit cells per dimension, q = (qx, qy) and qi = −1/2+ν/S+d are the wave vectors with d = −D, −D+1, ..., D} and ν = 0, 1, ..., S−1 (ν = 1/2, ..., S−1/2 for S odd). The eigenvectors of the Hamiltonian matrix of dimension (2D + 1) 4 S 4 represent the solution for the single-particle orbitals in the plane wave basis |q . The orthogonal optical lattice in the z direction separates. Assuming that only the central lattice plane is filled in the experiment, we can restrict ourselves to the Wannier function w (n) z (z) of band n in z direction. For the orbitals in Fig. 1 and 3  with suitable coefficients cij 25,38 . For the lattice depths V0 = 6ER and V0z = 25ER ( dip = 10, wdip = 4a), the tunneling element is t1 = 0.154ER, the next-nearest-neighbor tunneling t2 = −0.0094ER and the on-site interaction U = 32.9ERas/b (U (pz ) = 22.7ERas/b for the pz and s orbitals), where as the s-wave scattering length and b = λ/2 (cf. Fig. 3c, d and Fig. 4).