From Classical to Quantum Glasses with Ultracold Polar Molecules
FFrom Classical to Quantum Glasses with Ultracold Polar Molecules
Wolfgang Lechner ∗ and Peter Zoller Institute for Quantum Optics and Quantum Information,Austrian Academy of Sciences, 6020 Innsbruck, Austria andInstitute for Theoretical Physics, University of Innsbruck, 6020 Innsbruck, Austria (Dated: September 28, 2018)We study the dynamics of a bilayer system of ultracold polar molecules, which exhibits classicaland quantum glassy behavior, characterized by long tails in the relaxation time and dynamicalheterogeneity. In the proposed setup, quantum fluctuations are of the order of thermal fluctuationsand the degree of frustration can be tuned by the interlayer distance. We discuss the possibleobservation of a glassy anomalous diffusion and dynamical heterogeneity in experiment using internaldegrees of freedom of the molecules in combination with optical detection.
The recent experimental realization of cold ensemblesof polar molecules [1] has opened a new pathway to ex-plore the dynamics of quantum many body systems withstrong, long-range and anisotropic polar interactions[2–9]. In combination with low dimensional trapping ge-ometries, this allows the realization of stable many bodyphases, for example in the form of 1D wires or 2D pan-cakes [10–13], or as stacked pancakes representing cou-pled multilayer systems [14–18]. Most of the experimen-tal and theoretical studies of polar molecular gases havefocused so far on the quantum degenerate regime, and onequilibrium phenomena, including superfluid and crystalphases, quantum magnetism, and topological phases (fora review see [8, 9]). Instead, our interest below will beon non-equilibrium many-body dynamics. We will showthat a bilayer setup of ultracold polar bosonic moleculescan feature a glass phase , and we present methods toprepare this phase and measure the relevant order pa-rameters with tools available with present experimentalsetups. The unique feature and the theoretical challengeof glass physics with cold molecular ensembles is the pos-sibility to study the crossover from classical to quantumglasses.A glass phase is characterized by a plateau in the re-laxation time-scale, known as aging , with exponentiallyincreasing tails that prevent the system from reaching itsequilibrium state [13, 19–21]. In a structural glass, globalreorganization to the equilibrium is prevented from geo-metric frustration as a result of the dynamics [22]. Whilelarge reorientation is slow, on a local scale, relaxation canbe fast, a phenomenon known as dynamical heterogeneity [23]. In a classical glass this relaxation dynamics is domi-nated by thermal fluctuations, and a considerable under-standing of the glass transition has been gained from var-ious theoretical methods [6, 19, 25–28] and experimentalmodel systems [29–31]. In contrast, the question of theinfluence of quantum fluctuations on the glass relaxationdynamics is far less well understood. Recent theoreti-cal studies indicate, based on analytical [28, 32, 33] andnumerical [4, 5, 36] methods, that quantum fluctuations ∗ Electronic address: [email protected] s = 0.5as = as = 0.25 a s = 0.5as = as = 0.25 a ss (b)(c) (d) (e) (a) BABA
FIG. 1: (Color online) Bilayer setup of trapped polarmolecules with parallel dipole moments (a) and anti-paralleldipole moments (b). The lower layer represents a high-densitydipolar crystal, the upper layer is loaded with low density (seetext). A typical route to enter the glass phase is to quench thesystem from the liquid to a supercooled phase. In contrast,here, the glass phase is reached by decreasing the layer sepa-ration s . (c) Sketch of the phase diagram of a single layer ofbosonic polar molecules with the ratio between kinetic energyand interaction K and r d from Ref. [10]. Note, that the limit m → ∞ corresponds to the classical limit as the extension ofthe wavefunction vanishes with λ dB = h/ (2 πmk B T ) / . Thearrow indicates the parameter region of interest. (d) For par-allel dipoles [panel (a)], the effective volume and symmetry ofthe defects (dimers, trimers) depends on the interlayer sepa-ration. The ground state from a classical simulated annealingcalculation at T = 0 results in patterns include triangular,cubic and asymmetric structure where particles in the crystallayer (red) are displaced by the molecules in the defect layer(blue). (e) For anti-parallel dipoles [panel (b)] defect patternscan range from triangular, cubic to fivefold symmetries. can enhance but also inhibit the glass transition.In this Letter, we propose and analyze a protocol toprepare and measure a glass phase in a bilayer setup ofcold polar molecules (see Figs. 1a and b) in the regimeof the cross over from a classical to a quantum glass. Weassume ultracold molecules prepared in their electronic a r X i v : . [ c ond - m a t . qu a n t - g a s ] S e p and rovibrational ground states, where a static electricfield E ≡ E e z oriented perpendicular to the trapping lay-ers in the xy -plane induces an electric dipole moment d ≡ d e z . Thus the molecules will interact accordingstrong, long-range and anisotropic dipole-dipole interac-tion. The stacked pancake potential of Figs. 1a andb can be realized with a 1D optical lattice with layerseparation s controllable by laser parameters. For strongconfinement the motion of the molecule within each layercan be described as an effective 2D dynamics, while thetunneling between adjacent layers will be suppressed by asufficiently high barrier. The Hamiltonian for our bilayersystem thus has the form H = H A + H B + H AB . Here H A ( H B ) is the Hamiltonian for the intralayer dynamics, H X = (cid:88) i ∈ X (cid:18) p i m + V T ( r i ) (cid:19) + (cid:88) i,j ∈ X d X π(cid:15) r ij ( X = A, B )(1)as sum of kinetic energy and a tranverse trapping poten-tial, and the (purely repulsive) dipolar interaction within-plane distance r . The interlayer interactions are H AB = (cid:88) i ∈ A,j ∈ B d A d B π(cid:15) R ij (cid:32) − s R ij (cid:33) . (2)with s the layer separation and R the distance betweenthe molecules with R = r + s . Note, that our modelallows for different dipole moments d A and d B in layers A and B , respectively. Different effective dipole momentsin the two layers can be achieved e.g. by a gradient inthe electric field E z or by engineering of the interactionsusing internal rotational degrees of freedom [8]. As wewill see below, this feature is essential to observe clearsignatures of a glass phase.The basic steps of our protocol to study glassy dynam-ics are as follows. (i) We start with two uncoupled layers( s large), where layer A is in a high density 2D crystalphase and layer B in a low density gas phase. To preparethe initial state for the dynamics, the layer separation isquenched to a small distance s , which leads to the for-mation of defects with various symmetries and patternsdepending on the value of s (see Fig. 1). (ii) We evolvethe system to find glassy dynamics, identified by dynam-ical heterogeneity, the deviation from the linear diffusionlaw and a plateau in the relaxation (see Fig. 2). (iii)In order to measure these features in a possible experi-ment with cold molecules we introduce marker molecules ,i.e. molecules prepared in a different internal state, whichallow tracking the time evolution of the particles (see Fig.3). Preparation of the initial state - We first preparemolecules in two uncoupled layers: a dipolar crystal ofpolar molecules in layer A and a low density phase ofdefects in layer B . For uncoupled layers, the phase oflayer A is described by parameters r d = Dm/ ¯ h a and K = k B T a /D , where D = d A / (4 π(cid:15) ). The first pa-rameter is the ratio of the dipolar interaction D/a fora given mean intralayer distance between the particles a , and the kinetic energy ¯ h /ma , which in the dipolar crys-tal phase is r d (cid:29)
1. The second parameter measures thetemperature in units of the interaction. For bosons therelevant phase diagram is sketched in Fig. 1(c) [10], whichat low temperatures shows a low density (2D) superfluidphase, a high-density crystal phase and a high tempera-ture fluid phase. The theoretically predicted conditionsfor the formation of a self-assembled crystal are, in thecase of bosons, r d >
20 and
K < . r d and K are a requirement for tem-perature, density and dipole moment. While in presentexperiments with KRb the crystalline phase has not beenrealized, the ongoing effort in the laboratory to preparecold ensembles of LiCs molecules [37] with its large elec-tric dipole moment µ = 5 . m = 92uprovides a promising candidate to obtain dipolar crys-tals. In this case for a given temperature T = 0 . µ Kthe requirement on r d and K corresponds to an effec-tive (induced) dipole moment of d A = 2 . a = 0 . µ m. We assume that the density in layer B islower than that in layer A while all other parameters arethe same. In the following we find it convenient to definelength, energy and time in the reduced units a , D/a , a (cid:112) ma /D , respectively.To prepare the initial state for the glassy dynamics, thetwo layers are brought together by a quench in the layerseparation. The speed of the quench is chosen such, thatnon-adiabatic effects are negligible [38]. We consider twocases: the dipoles can be aligned parallel [see Fig. 1(a)] oranti-parallel [see Fig. 1(b)]. Fig. 1(d) depicts the result-ing defect patterns for parallel dipoles d B = 4 d A wheremolecules belonging to different layers attract each other.This attraction leads to stable dimers or trimers [8] if thetwo layers are separated by s ≤ a . These composite de-fects are reminiscent of interstitials in a crystal. Here,the effective volume of the interstitials can be tuned bythe ratio of the dipole moments d A /d B and the inter-layer separation s . For various layer separations rangingfrom s = 0 .
25 to s = a the resulting patterns includehexagonal, cubic and asymmetric defects. We note thatthe binding energy of similar impurities has been recentlystudied in a fermionic system [39].In the anti-parallel case, d B = − d A , the interlayer in-teraction is repulsive. This can be achieved using internalrotational degrees of freedom [8]. The defect patternsdepicted in Fig. 1(e) include cubic, five-fold symmet-ric and triangular symmetries reminiscent of vacancy de-fects. Again, the effective volume of the defects can betuned by the layer separation.In both cases, parallel and anti-parallel bilayer setup,the combined system is a mixture of effective defects andmolecules. This is an experimental realization of a binarymixture of dipoles which is a well known glass-formingliquid in the classical regime [30, 31]. We note, thatthe anti-parallel setup allows one to induce defects withfive-fold symmetry which may allow one to implement aclassical spin liquid model, where fivefold symmetries arestudied [21], in the quantum regime. Glass Dynamics - Below we study the relevant glass or-der parameters from the dynamics of the proposed setupafter the quench. In the classical regime we use moleculardynamics simulations. In the quantum regime we employthe recently developed dynamical path integral methods,which has been applied to the glass transition in Ref. [4](for details on the numerical algorithms see SM [38]).We note that these methods include quantum effects butignore the exchange statistics between the particles andthat we are interested mainly in the transition from theclassical phase to a phase where quantum fluctuationsbecome important.A glass phase is identified by a dramatic increase in therelaxation time which divergences with τ r ∝ exp[ − / ( T − T g )] when approaching the ideal glass temperature T g [19]. This corresponds to an extended plateau in theself-intermediate scattering function [19], a two-time cor-relation function defined as F ( k ∗ , ∆ t ) = 1 N (cid:88) j (cid:104) e ik ∗ [ r j ( t ) − r j ( t +∆ t )] (cid:105) . (3)Here, angular brackets (cid:104)·(cid:105) denote ensemble averages overmany realizations of the experiment after the quench.The sum runs over all particles in both layers N = N A + N B and k ∗ = | k ∗ | is the absolute value of the char-acteristic k -vector of the coupled system correspondingto the first peak in the static structure factor S ( k ) = N (cid:80) i (cid:80) j (cid:104) e − i k [ r i − r j ] (cid:105) where r is the distance in the xy -plane. Note that Eq. (3) bears some similarity with theFourier transform of the density-density correlation func-tion, however here, the sum runs over individual particlesat different times. In the path integral picture, the anal-ogous order parameter to the classical case Eq. (3) readsas [4]Φ( k ∗ , ∆ t ) = 1 N ¯ hβ (cid:90) ¯ hβ dλ (cid:104) ρ † ( t + ∆ t + iλ ) ρ ( t ) (cid:105) . (4)Here, ρ ( k , ∆ t ) = (cid:80) Nj e i kr j and the integral runs over theimaginary path integral time. The semi-classical approx-imation allows one to follow the real time propagation ofthe path integral [2] neglecting exchange (for details seeSM [38]).A crucial dynamical characteristic of the glass phase isits deviation from the linear diffusion laws. In particu-lar, the local diffusion of individual particles is spatiallyheterogeneous, a characteristic of the glass phase knownas dynamical heterogeneity . The local diffusion of parti-cle j is ∆ r j (∆ t ) = (cid:104) r j ( t + ∆ t ) − r j ( t ) (cid:105) and the meansquared displacements (cid:104) ∆ r (∆ t ) (cid:105) = 1 /N (cid:80) i ∆ r i (∆ t ).In the presence of quantum fluctuations the analogousmeasure in the path integral picture is the root meansquared displacement of the projections of the path in-tegral with ˆ r j ( t ) = N ¯ hβ (cid:82) ¯ hβ dλr j ( t + iλ ) and ∆ˆ r j (∆ t ) = (cid:104) ˆr j ( t + ∆ t ) − ˆr j ( t ) (cid:105) . (c) (d) -2 -1 (b) -2 -1 -3 -2 -1 (a) -1 s=0.2as=0.3as=0.4as=0.5a -1 classical limit parallelanti-parallel FIG. 2: (Color online) Glass order parameters in the classical(top) and semi-classical regime from path integral calculations(bottom). (a) In the glass phase, the average mean squareddisplacement of the particles deviates from the constant dif-fusion. In the anti-parallel setup the plateau extends over thewhole sampling time below
K < .
02 of the combined sys-tem. (b) The relaxation of the time dependent structure fac-tor F ( k ∗ , t ) for s = 0 . a with d B = 4 d A (black) and d B = − d A (red) diverges when approaching the glass phase [symbols asin (a)]. (c) Relaxation dynamics with initial temperature K = 0 .
01 in the anti-parallel setup for various choices of effec-tive r d of the combined system. Parameters are chosen suchthat the classical limit m → ∞ enters a glass phase and forsmaller r d the system melts due to large quantum fluctuations(blue). This allows one to study the cross-over from a clas-sical glass to a quantum phase. (d) Relaxation dynamics asa function of the layer separation with effective r d ≈ − K = 0 .
02 in the anti-parallel setup. While for distance s = 0 . a the system is a liquid due to large quantum fluctu-ations (blue) one reaches the glass phase when approaching s = 0 . a (black). The number of particles used was N A = 780and N B = 200 in the classical simulations and N A = 168 and N B = 50 in the path integral simulations using 50 time-slicesper particle and periodic boundary conditions. Averages aretaken from 100 independent runs. For LiCs a time of ∆ t = 10 in reduced units corresponds to ∆ t = 3 . The unique feature of the bilayer system of ultracoldmolecules is the possibility to study glassy dynamics inboth the classical and quantum regime, which are charac-terized by the dominance of thermal vs. quantum fluctu-ations, respectively. This transition can be controlled bythree tunable parameters: s representing the layer sepa-ration, the dimensionless temperature K , and the dipolarcrystal parameter r d . The case of large r d ≡ Dm/ ¯ h a corresponds to the classical limit [compare Fig. 1(c)].The associated classical glass transition can be studiedby varying K with fixed r d , while varying the initial de-fects with s [see Fig. 1(d) and (e)]. Below we will discussthis limit in Fig. 2(a,b). The transition from the clas-sical to the quantum regime is achieved by decreasing r d . This increasing role of quantum fluctuations on theglassy dynamics is summarized in Fig. 2(c) for various r d , while Fig. 2(d) discusses the effect to defect patterns,as controlled by s [compare Fig. 1(d) and (e)].The averaged root mean square displacement as a func-tion of time is shown in Fig. 2(a). The dynamics does notfollow the Einstein diffusion law (cid:104) ∆ r (∆ t ) (cid:105) (cid:54) = 2 Dd ∆ t ,where d is the dimension of the system, but developsa plateau larger than the sampling time, which cor-responds to a divergence in the thermodynamic limit,for an effective K of the combined system smaller than K < .
02. Fig. 2(b) shows the according relaxation dy-namics F ( k ∗ , ∆ t ) for a final layer separation s = 0 . a inthe classical regime. Both setups with parallel and anti-parallel dipoles show a growing plateau and enter a glassphase.Including quantum fluctuations the relaxation dynam-ics of the system changes considerably. Fig. 2(c) de-picts the relaxation dynamics Φ( k ∗ , ∆ t ) calculated froma path integral simulation as a function of r d at a fixed K . With increasing r d , the dynamics approaches that ofthe classical molecular dynamics simulation correspond-ing to r d → ∞ where the system is glassy. By lowering r d , the system melts due to large quantum fluctuations.In the bilayer system, one can also tune the defect sizewith the layer separation (see Fig. 1). Fig. 2(d) shows aset of parameters, where one can reach the glass phaseby variation of the layer separation. In contrast to theclassical model system of colloids [31], the actual dynam-ics after the initial thermalization corresponds to that ofan isolated system (micro-canonical). However, this doesnot influence the long tail part in the relaxation which isindependent of the microscopic dynamics [41]. Measurement of Dynamical Heterogenity - One way tomeasure the anomalous diffusion as a glass feature is toconsider a subensemble of polar molecules in a small spa-tial region ( marker molecules ), which are addressed witha laser and transferred to another internal (e.g. hyper-fine) state [4]. The position of this molecular cloud ata later time can be measured with optical techniques.The ultimate tool for such position measurements will beprovided by the quantum gas microscope for molecules,which - following the atomic case - will provide single site/ single molecule detection with a resolution of the singlesite of an optical lattice [42]. Dynamical heterogeneitycan then be measured as follows: We mark molecules inthe crystal layer A in a small region with laser waist-size R > a . After the system evolves in time, the positionsof the marker molecules is measured, and the extensionof the cloud of marker molecules is given by the radius ofgyration R g = 1 / (2 N ) (cid:80) i,j ( r i − r j ) . In a liquid a cloudof tagged particles spreads out linearly with time whilein an amorphous solid it is constant R g /R = 1. Mostimportant, in a liquid and crystal, the dynamics is inde-pendent of the initial position of the makers molecules asthe dynamics is spatially homogeneous. However, in theglass phase due to dynamical heterogeneity the scalingdiffers significantly depending on the initial position (seeFig. 3). As the initial size is R from the laser waist, (a) xy slow fast (c) (d) (b) FIG. 3: (Color online) Dynamical heterogeneity of a quan-tum glass in an AMO experiment: Molecules in an excitedinternal state act as marker molecules to identify glassy dy-namics. (a) Typical initial configuration of layer A (red) andlayer B (blue) depicted as the projection of the paths ontoreal space. (b) The parameter ψ i measures the relative or-der in the vicinity of particle i (for definition and details seeSM [38]). Here, regions of large order ( ψ i ≈
1) are regionsof small classical diffusion and a small extension of the parti-cle wavefunction. Regions of low order ( ψ i ≈
0) feature largediffusion and larger quantum fluctuations. Note, that dynam-ical heterogeneity is a more general feature of glasses, which isalso present in systems without apparent partial crystallinity(see SM [38]). (c) Snapshots of two possible initial choicesof marker molecule positions (green and orange) at ∆ t = 0with R g = R . The particle positions indicated by spheresrepresent the maximum of the path probability in layer A (gray) and layer B (lightblue). (d) After evolving the sys-tem in time, the radius of gyration of marker molecules aftertime ∆ t = 1 . × (5ms) differs depending on the initialposition. In the glass phase, due to dynamical heterogeneity,the growth of R g strongly depends on the choice of the initialposition. The cloud of maker molecules placed in an inactiveregion (green) of the glass follows the dynamics of an amor-phous crystal R g = R while molecules in a mobile region(orange) show a fast increase in size (inset). a series of single (destructive) measurements after time∆ t allows one to distinguish a glass from liquid or amor-phous solid (see Fig. 3). Note that the marker moleculesare quantum mechanically distinguishable which changesthe dynamics in the deep quantum regime when exchangestatistics is included, which is not relevant in the presentcase.In conclusion, we have shown that a bilayer system ofpolar molecules can, with properly chosen parameters,exhibit a glass phase in a regime where quantum fluctua-tions are of the order of thermal fluctuations. Thus ultra-cold ensembles of polar molecules could provide a widelytunable paradigmic model system for classical and quan-tum glass physics, where theory and experiment can meetin a yet unexplored parameter regime of glass physics,providing in particular a stimulus for the theoretical de-velopments. We thank K. Binder and M. Baranov for fruitful dis-cussions. Work supported by the Austrian Science Fund(FWF): P 25454-N27, SFB FOQUS, Marie Curie Ini-tial Training Network COHERENCE, and ERC SynergyGrant UQUAM. [1] References in special issue: D. S. Jin and J. Ye, Chem.Rev. , 4801 (2012).[2] L. D. Carr, D. DeMille, R. V. Krems, and J. Ye, New. J.Phys. , 231 (2008).[4] S. Ospelkaus, K.-K. Ni, G. Quemener, B. Neyenhuis, D.Wang, M. H. G. de Miranda, J. L. Bohn, J. Ye, and D. S.Jin Phys. Rev. Lett. , 030402 (2010).[5] B. Yan, S. A. Moses, B. Gadway, J. P. Covey, K. R. A.Hazzard, A. M. Rey, D. S. Jin, J. Ye, arXiv:1305.5598(2013).[6] J. Deiglmayr, A. Grochola, M. Repp, K. M¨ortlbauer, C.Gl¨uck, J. Lange, O. Dulieu, R. Wester, and M. Wei-dem¨uller, Phys. Rev. Lett. , 133004 (2008).[7] J. G. Danzl, E. Haller, M. Gustavsson, M. J. Mark, R.Hart, N. Bouloufa, O. Dulieu, H. Ritsch, and H.-C. N¨agerl,Science , 1062 (2008).[8] M. A. Baranov, M. Dalmonte, G. Pupillo and P. Zoller,Chem. Rev. , 5012 (2012).[9] T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, and T.Pfau, Rep. Prog. Phys. , 126401 (2009).[10] H. P. B¨uchler, E. Demler, M. Lukin, A. Micheli, N.Prokof’ev, G. Pupillo, and P. Zoller, Phys. Rev. Lett. ,060404 (2007).[11] G. E. Astrakharchik, J. Boronat, I. L. Kurbakov, and Y.E. Lozovik, Phys. Rev. Lett. , 060405 (2007).[12] N. R. Cooper and G. V. Shlyapnikov, Phys. Rev. Lett.103, 155302 (2009).[13] M. Lewenstein, A. Sanpera, V. Ahufinger, “UltracoldAtoms in Optical Lattices: Simulating quantum many-body systems”, Oxford University Press, Oxford (2012).[14] G. Pupillo, A. Griessner, A. Micheli, M. Ortner, D.-W.Wang, and P. Zoller, Phys. Rev. Lett. , 050402 (2008).[15] D.-W. Wang, M. D. Lukin, and E. Demler, Phys. Rev.Lett. , 180413 (2006).[16] A. Pikovski, M. Klawunn, G. V. Shlyapnikov, and L.Santos, Phys. Rev. Lett. bx105, 215302 (2010).[17] A. C. Potter, E. Berg, D.-W. Wang, B. I. Halperin andE. Demler, Phys. Rev. Lett. , 220406 (2010).[18] M. Knap, E. Berg, M. Ganahl, and E. Demler, Phys.Rev. B , 064501 (2012).[19] K. Binder and W. Kob, “Glassy Materials and Disor-dered Solids”, World Scientific Publishing, London (2011).[20] P. G. Debenedetti and F. H. Stillinger, Nature , 259(2001).[21] H. Shintani and H. Tanaka, Nature Phys. , 200 (2006).[22] For glasses in AMO setups as a result of disorder see:T. Roscilde and J. I. Cirac, Phys. Rev. Lett. , 190402(2007), L. Sanchez-Palencia and M. Lewenstein, NaturePhys. , 87 (2010), S. Gopalakrishnan, B. L. Lev, and P.M. Goldbart, Phys. Rev. Lett. , 277201 (2011). [23] D. Chandler and J. Garrahan, Annu. Rev. Phys. Chem. , 191 (2010).[24] H. Tanaka, T. Kawasaki, H. Shintani, and K. Watanabe,Nature Mat. , 324 (2010).[25] G. H. Fredrickson and H. C. Andersen, Phys. Rev. Lett. , 1244 (1984).[26] S. Whitelam, L. Berthier, and J. P. Garrahan, Phys. Rev.E , 026128 (2005).[27] B. Coluzzi, G. Parisi and P. Verrocchio, Phys. Rev. Lett. , 306 (2000).[28] V. Lubchenko, and P. G. Wolynes, Adv. Chem. Phys. , 95 (2007).[29] G. L. Hunter and E. R. Weeks, Rep. Prog. Phys. , 627 (2000).[31] F. Ebert, P. Keim, and G. Maret, Eur. Phys. J. E ,161 (2008).[32] V. N. Novikov and A. P. Sokolov, Phys. Rev. Lett. ,065701 (2013).[33] E. Rabani and D. R. Reichman, Annu. Rev. Phys. Chem. , 157 (2005).[34] T. E. Markland, J. A. Morrone, B. J. Berne, K. Miyazaki,E. Rabani, and D. R. Reichman, Nature Phys. , 134(2011).[35] T. E. Markland, J. A. Morrone, K. Miyazaki, B. J. Berne,D. R. Reichman, and E. Rabani, J. Chem. Phys. ,074511 (2012).[36] M. Boninsegni, L. Pollet, N. Prokof’ev, and B. Svistunov,Phys. Rev. Lett. , 025302 (2012).[37] M. Repp, R. Pires, J. Ulmanis, R. Heck, E. D. Kuhnle,and M. Weidem¨uller, and E. Tiemann, Phys. Rev. A ,010701 (2013);S.-K. Tung, C. Parker, J. Johansen, C.Chin, Y. Wang, and P. S. Julienne, Phys. Rev. A ,010702 (2013).[38] see Supplementary Material.[39] N. Matveeva and S. Giorgini, arXiv :1306.5588 (2013).[40] M. Ceriotti, M. Parrinello, T. E. Markland, and D. E.Manolopoulos, J. Chem. Phys. , 124104 (2010).[41] T. Gleim, W. Kob, and K. Binder, Phys. Rev. Lett. ,4404 (1998).[42] J. Ye, private communication. For the atomic quantumgas microscope see: W. S. Bakr, J. I. Gillen, A. Peng, S.F¨olling and M. Greiner, Nature , 74 (2009); C. Weit-enberg, M. Endres, J. F. Sherson, M. Cheneau, P. Schauß,T. Fukuhara, I. Bloch, S. Kuhr, Nature , 319 (2011) Supplemental Material to
From Classical To Quantum Glasses with Ultracold Polar Molecules
I. NUMERICAL METHODSA. Classical Dynamics
We use a second-order Langevin integrator [1] to study the dynamics of the bilayer system in the classical regime.The time evolution is given by the Langevin equation m ∂ x ∂t = −∇ V ( x , ..., x ) − γm ∂ x ∂t + (cid:112) k B T γmη ( t ) . (5)Here, W corresponds to a Wiener process and η ( t ) = ˙ W . This equation can be solved numerically in second orderaccuracy with the following updates v ( t + ∆ t/
2) = v ( t ) + 12 dtF ( x ( t )) − dtγv ( t ) + (6)+ 12 √ dtσζ n − dt γ ( F ( x ( t )) − γv ( t )) −− h / γσ ( 12 ζ n + 1 √ η n ) x ( t + ∆ t ) = x ( t ) + dtv ( t + ∆ t/
2) + dt / σ √ η n (7) v ( t + ∆ t ) = v ( t + ∆ t/
2) + 12 dtF ( x ( t + ∆ t )) − (8) − dtγv ( t + ∆ t/
2) + 12 √ dtσζ n −− dt γ ( F ( x ( t + ∆ t )) − γv ( t + ∆ t/ −− h / γσ ( 12 ζ n + 1 √ η n ) (9)Here, η and ζ are two independent noise variables from a Wiener process and F ( x ( t ) = (cid:80) ij ∇ V ij . Note that in thelimit of zero-friction γ →
0, the equations of motion reduce to a energy conserving (micro-canonical) propagator.
B. Semi-classical Dynamics
For completeness we repeat here the derivation of the numerical integrator used in the semi-classical regime givenin Ref. [2]. The partition function is rewritten in terms of Feynman path integrals (cid:104) A (cid:105) = 1 Z tr [ e − βH A ] , (10)with Z = 1(2 π ¯ h ) f (cid:90) d f p (cid:90) d f q e − β n H n ( q , p ) . (11)Here, f = N n and β n = β/n . The full Hamiltonian consists of system Hamiltonian and interaction Hamiltonian H n ( q , p ) = H n ( q , p ) + V n ( q ) . (12)The system part is H n ( q , p ) = N (cid:88) i =1 n (cid:88) j =1 (cid:32) ( p ( j ) i ) m i + 12 m i β n ¯ h [ q ( j ) i − q j − i ] (cid:33) . (13)While sophisticated semi-classical methods have been developed to sample this Hamiltonian without dissipation, thecoupling to a thermal bath is challenging as the frequencies in the path space separate from the frequencies from thethermal bath. In the framework of path integral Langevin dynamics, this is solved by splitting the dynamics into twoparts: First, the full Hamiltonian evolves according to the path integral formalism in absence of a thermal bath. Ina second step, the thermal bath acts on the slowest and higher normal modes of the paths separately, on the level ofthe equations of motion.The energy conserving Liouvillian part of the dynamics (without dissipation) is written as a sum of system andinteraction L = L + L V . (14)Here, L = − [ H n ( q , p ) , . ] and L V = − [ V n ( q , p ) , . ]. Using Totter’s theorem, the propagator for a small ∆ t can bewritten as e ∆ tL = e (∆ t/ L V e ∆ tL e (∆ t/ L V . (15)The semi-classical treatment of the heat bath requires the normal mode represent of the Hamiltonian which is derivedfrom the transformation ˜ p k = n (cid:88) j =1 p ( j ) i C jk a nd ˜ q k = n (cid:88) j =1 q ( j ) i C jk , (16)where the matrix elements of c jk are C jk = (cid:112) /n, if k = 0 (cid:112) /n cos(2 πjk/n ) , if 1 ≤ k ≤ n/ − (cid:112) /n ( − j , if k = n/ (cid:112) /n sin(2 πjk/n ) , if n/ ≤ k ≤ n − . (17)In the normal mode representation the Hamiltonian reads as H n ( q , p ) = N (cid:88) i =1 n − (cid:88) k =0 (cid:32) (˜ p ( j ) i ) m i + 12 m i sin( kπ/n ) β n ¯ h (˜ q ki ) (cid:33) . (18)Dissipation is then added on the level of the equations of motion in the space of normal modes. This leads to thefollowing equations of motion e ∆ tL = e (∆ t/ L γ e (∆ t/ L V e ∆ tL e (∆ t/ L V e (∆ t/ L γ (19)While L and L V are associated with momentum and position dynamics in the path integral, L γ is defined in thespace of normal mode variables. This requires four transformation between normal mode and path integral space ineach time step. The Trotter step associated with L V is p ( j ) i ( t + ∆ t ) = p ( j ) i ( t ) − ∆ t ∂V∂q ( j ) i (20)The evolution of the path integral with L defined in normal mode space corresponds to˜ p ( k ) i ( t + ∆ t ) = ˜ p ( j ) i ( t ) cos( ω k ∆ t ) − ˜ q ( j ) i ( t ) m i ω k sin( ω k ∆ t ) . (21)and ˜ q ( k ) i ( t + ∆ t ) = ˜ p ( j ) i ( t ) 1 m i ω k sin( ω k ∆ t ) − ˜ q ( j ) i ( t ) cos( ω k ∆ t ) . (22)In the last step, L γ describes thermalization in the normal mode space and leads to the updates˜ p ( k ) i ( t + ∆ t ) = ˜ p ( k ) i ( t ) e − ∆ t/ γ ( k ) + (cid:114) m i β n (1 − e − ∆ tγ ( k ) ) ξ ( k ) i . (23)Fig. 4a and b show typical snapshots from a path integral simulations in the large r d limit (classical) compared to a (a) (b) FIG. 4: Snapshots from a path integral simulation with d B = − d A , K = 0 .
02 deep in the crystalline region r d = 40 (a) andclose to the superfluid region r d = 25 (b) (see Fig. 2 in the main text) . The projection of the paths onto real space of layer A (red) and layer B (blue) indicate the extension of the wave function. phase where quantum fluctuations become important. II. STRUCTURE FACTOR AND INTERMEDIATE SCATTER FUNCTION
The dynamical slowing down of the glass phase corresponds to the development of a plateau in the time dependentintermediate scatter function F ( k, t ) (classical) and Φ( k, t ) (path integral). This time dependent order parameterallows one to measure the relaxation time of a characteristic k -vector. This characteristic k -vector ( k ∗ ) is associatedwith the first peak in the static structure factor S ( k ) and taken as an input parameter for the time dependentintermediate scatter function of the Kubo transformed density Φ( k ∗ , t ) (or F ( k ∗ , t )). The static structure factor S ( k )is defined as S ( k ) = 1 N (cid:104) (cid:88) l,m e − i k ( R l − R m ) (cid:105) . (24)Here, R i is the position of the classical particle or in the quantum case the centroid of the path integral R i = 1 /n (cid:80) n r ni ,where the average runs over all imaginary time slices. Finally, the time-dependent density-density correlation functionis Φ( k ∗ , ∆ t ) = 1 N ¯ hβ (cid:90) ¯ hβ dλ (cid:104) ρ † ( t + ∆ t + iλ ) ρ ( t ) (cid:105) (25)Here, ρ ( k, t ) = (cid:80) Nj e i kr j and the integral runs over the imaginary path integral time. The ensemble averages (cid:104)·(cid:105) aretaken over a large number of realizations. l d B analyticcrystal FIG. 5: The extension of the particle wave function of a free particles with unit temperature and mass is in units of ¯ h : λ dB = (cid:112) / π = 0 . a = 1 (density ρ = 2 / √ III. COMPARISON OF PATH INTEGRAL AND WIGNER FUNCTION FORMALISM
The extension of the wavefunction of a free particle with mass m at temperature T is λ dB = (cid:115) h πmk B T . (26)In the path integral formalism, the average extension of the particles is calculated from λ = 2 (cid:104) N ¯ hβ (cid:90) ¯ hβ dλ | r ( t + iλ ) − ˆ r j ( t ) | (cid:105) . (27)To check for consistency we compare λ dB with the numerical results. Fig. 5 depicts the particle extension as afunction of the interaction strength D . In the free particle limit D → λ dB . The extension of the wavefunction can also be calculated from the width of the Wigner function witha thermal characteristic function after integrating out the momentum. The Wigner function of a thermal state is [3] W ( α, α ∗ ) = (28)= 1 π (cid:90) d λ exp[ − λα ∗ + λ ∗ α ] × exp[ − | λ | ] exp[ −| λ | e ¯ hω/k B T − π tanh[ ¯ hω k B T ] exp {− | α | tanh[ ¯ hω k B T ] } Using α = ip ¯ h √ η + x (cid:114) η η = √ km/ ¯ h (30)we rewrite Eq. 28 in terms of position and momentum W ( x, p ) = h W ( α, α ∗ ). Integrating out the momentum leadsto W ( x ) = (cid:90) ∞−∞ W ( x, p ) dp = (31)= (cid:114) η π (cid:114) tanh[ ¯ hωk B T ] exp {− tanh[ ¯ hωk B T ] ηx } W ( x ) to the Gaussian distribution of a thermalwavefunction P ( x ) = 1 (cid:112) πλ dB exp[ − x λ dB ] (32)Expanding the tanh in Eq. (31) and comparing the coefficients of Eq. 31 and 33 we find¯ h πmk B T = ¯ hωk B T η (33)which is consistent with Eq. (30).In a crystal the potentials of neighboring sites act as a trapping potential which localizes the wavefunctions (seeFig. 5). (a) (b) (c)
FIG. 6: (Color online) Three different representations of a typical configuration in the glass phase with partial crystalline order.(a) Projection of the path integral density to real space with particles in the crystal layer (red) and particles from the defectlayer (blue). The regions of large delocalization are correlated with defect particles while particles in crystalline regions aremore localized. (b) Maximum of the probability positions of particles as spheres colored by the value of the time independentmedium range order parameter ˆΨ ranging from unordered regions (red) to maximal ordered regions (blue). (c) Particlescolored by the time dependent mean squared displacement ˆ D i = D i − D min D max − D min ranging from frozen regions (red) to regions oflarge displacement (blue) . IV. HEATING DUE TO THE QUENCH
The system is described by the Hamiltonian Eq. (1) in the main text. Here, we explicitly include the time dependentparameter λ which allows one to vary the interlayer separation H ( λ ( t )) = (cid:88) i p m + (cid:88) A D (cid:112) x + y + (34) (cid:88) B D (cid:112) x + y + (cid:88) i,j D s ( λ ( t )) x + y + s ( λ ( t )) − (cid:112) x + y + s ( λ ( t ) . Here, the sum over A and B indicates interactions between particles of the same layer and the last term the interactionbetween particles in different layers. The layer separation is then linearly switched from s i nit to s f inal with s ( λ ( t )) = λ ( t ) s i nit + (1 − λ ( t )) s f inal (35)1 * /T P t=850 m st=170 m st = 95 m s FIG. 7: Probability distribution of the effective temperature T ∗ /T after a quench from s (0) = a to s ( t ) = 0 . a for variousquench velocities. Densities and temperature a chosen as in case anti-parallel dipoles depicted in Fig. 1 in the main text. Fortotal quench times larger than t = 850 µs the temperature is basically unchanged while a fast quench with t < µs leads toconsiderable non-adiabatic dissipative work. The distributions is sampled from 50000 runs. by the protocol λ ( t ) = v l t. (36)The non-equilibrium work performed during the quench is the sum of reversible work and irreversible work W = W r ev + W i rr or explicitly W = (cid:90) dλ ∂λ∂t ∂H ( λ ( t )) λ (37)The irreversible work from the quench increases with the speed of the protocol and heats up the system non-adiabatically. The protocol is a follows: First, the molecules are prepared at temperature T as described in themain text. By a rapid quench in the interlayer separation, the system is driven out of equilibrium. During thequench, the system evolves micro-canonically and work is performed from the change in the layer separation s ( λ ). Wedefine an effective temperature in the micro-canconical ensemble as the average kinetic energy T ∗ = 2 KN dk B . (38)Here, K = Nm (cid:104) v (cid:105) is the kinetic energy and d = 2 the dimensionality of the system. Fig. 7 depicts the probabilitydensity of the effective temperature in the final state. For a slow quench with switching time larger than t > µs the temperature increase is negligible. V. DYNAMICAL HETEROGENEITY VS. STRUCTURAL ORDER
The local diffusion D i = (cid:104) [ r i ( t + ∆ t ) − r i ( t )] (cid:105) reveals a remarkable connection between structure and dynamicsin the glass phase as was recently shown in colloidal glasses [6, 7]. In particular, for some parameters the mobilityof individual particles D i is correlated with the local structure of the neighbors. This local order is measured in twodimensions by the averaged hexatic order parameterˆΨ = 1 N b (cid:88) N b | Ψ i | (39)with Ψ i = 1 N b (cid:88) N b e i φ ij . (40)2Here, the sum runs over all neighbors N b including all particles within a distance r n of particle i .The dynamical D i and the structural order parameter ˆΨ i are homogeneous in a liquid or a solid but are locallyheterogeneous and anti-correlated in the glass phase. Fig. 6(b) and (c) show the mean position of the particlewavefunctions where colors represent the value of ˆΨ i and ˆ D i , respectively, with the normalized diffusion ˆ D i = D i − D min D max − D min .In the path integral picture, structure and dynamics is also correlated with the extension of the particles inimaginary time. Fig. 6(a) shows the projection of the paths from the path integral simulation. The extensionof the quantum fluctuations are calculated from the second moment of the normal modes of the path integral λ = 2 (cid:104) N ¯ hβ (cid:82) ¯ hβ dλ | r ( t + iλ ) − ˆ r j ( t ) | (cid:105) which is consistent with the analytic result of the extension of the Wignerfunction from a thermal quantum characteristic function (see section III). The extension of a free particle at finitetemperature is given by the thermal deBroglie wavelength λ = [ h / (2 πk B T m )] / while in a crystal the extension issmaller due to the interaction of neighboring particles acting as a trapping potential which is to first order quadratic.In the glass, the structure and therefore the effective potential at each site is heterogeneous which is also apparent inthe wavefunction extension. [1] E. Vanden-Eijnden and C. Ciccotti, Chem. Phys. Lett. , 310 (2006).[2] M. Ceriotti, M. Parrinello, T. E. Markland, and D. E. Manolopoulos,
J. Chem. Phys. , 124104 (2010).[3] C. W. Gardiner and P. Zoller,
Quantum Noise , Springer (2010).[4] T. E. Markland, J. A. Morrone, B. J. Berne, K. Miyazaki, E. Rabani, and D. R. Reichman,
Nat. Phys. , 134 (2011).[5] T. E. Markland, J. A. Morrone, K. Miyazaki, B. J. Berne, D. R. Reichman, and E. Rabani, J. Chem. Phys. , 074511(2012).[6] H. Tanaka, T. Kawasaki, H. Shintani, and K. Watanabe,
Nature Materials , 324 (2010).[7] K. Watanabe, T. Kawasaki,and H. Tanaka, Nature Materials10