Dynamic Stimulation of Quantum Coherence in Lattice Bosons
DDynamic Stimulation of Quantum Coherence in Systems of Lattice Bosons
Andrew Robertson, Victor M. Galitski, and Gil Refael Joint Quantum Institute and Condensed Matter Theory Center, Department of Physics,University of Maryland, College Park, Maryland 20742-4111, USA Department of Physics, California Institute of Technology,1200 E. California Boulevard, Pasadena, California 91125, USA (Dated: October 29, 2018)Thermal fluctuations tend to destroy long-range phase correlations. Consequently, bosons in alattice will undergo a transition from a phase-coherent superfluid as the temperature rises. Contraryto common intuition, however, we show that non-equilibrium driving can be used to reverse thisthermal decoherence. This is possible because the energy distribution at equilibrium is rarely optimalfor the manifestation of a given quantum property. We demonstrate this in the Bose-Hubbard modelby calculating the non-equilibrium spatial correlation function with periodic driving. We showthat the non-equilibrium phase boundary between coherent and incoherent states at finite bathtemperatures can be made qualitatively identical to the familiar zero-temperature phase diagram,and we discuss the experimental manifestation of this phenomenon in cold atoms.
PACS numbers: 64.60.Ht, 03.75.Lm, 03.75.Kk
A system of bosons confined to a lattice has long rep-resented an alluring opportunity to study the interplaybetween two phenomena at the heart of the academicand industrial interest in many-body quantum mechan-ics: particle tunneling and phase coherence. The Bose-Hubbard model (BHM) describing such systems in thetight-binding approximation is much richer than its sim-ple mathematical form betrays. It admits such noveltiesas dynamic localization [1, 2], photon-assisted tunneling[3, 4], as well as an archetypal example of a quantumphase transition between superfluid and insulator-likestates [5, 6]. Thermal fluctuations destroy these quan-tum effects, but it is comparatively unknown that delib-erately driving the system out of equilibrium can mod-erate or reverse entirely the destructive effect of raisingthe temperature.Despite its obscurity, it has been known since the1960’s that pushing a system out of equilibrium can en-hance its quantum properties. In 1966, Wyatt et. al. [7]showed that illuminating a microbridge could stimulateits superconductivity. Eliashberg explained this result in1970 by calculating the nonequilibrium quasiparticle dis-tribution induced by the radiation [8]. Blamire et. al.later demonstrated that superconducting transition tem-peratures could be enhanced in this manner by severaltimes their equilibrium values [9]. More recently, the ideaof non-equilibrium phase transitions has garnered inter-est in studies of optically trapped atoms [10] and inducedtopological structure [11].The reason for enhancement goes as follows. The quan-tum properties of a system depend on the energy distri-bution of excitations. It is easily verified that the equilib-rium distribution is rarely optimal for the enhancementof a chosen property. A brief survey of the model at equi-librium reveals this to be the case for the BHM. Indeed,an analogue of photon-assisted tunneling (PAT) [2, 4] has already been theorized for the BHM at T = 0. However,a fully non-equilibrium treatment including the effectsof temperature (and how they may be mitigated) is notknown to us.In this Letter, we shall show that harmonically drivinga system of lattice bosons connected to a thermal reser-voir can increase the region in parameter space wherethe quantum coherent phase exists. Even for finite tem-peratures of the bath, the phase diagram of the BHMcan be made qualitatively identical to the T = 0 di-agram. We shall demonstrate this by defining non-equilibrium correlation functions (cid:104) a † i ( t ) a j ( t (cid:48) ) (cid:105) within theKeldysh and Floquet formalisms [12–16]. Divergence ofthe real part of this quantity corresponds to infinite long-range correlations. This will define the phase boundarybetween superfluid and incoherent states for our non-equilibrium system. We shall find these functions per-turbatively [17, 18] in the small quantity J/U and arriveat a Dyson equation that has both ordinary and entry-wise (Hadamard) matrix products. This novel structurewill then be solved by column vectorization for the sta-tionary non-equilibrium correlation function.To see how superfluidity may be enhanced by meansof a non-equilibrium pulse, let us briefly review the BHMat equilibrium. The hamiltonian of the BHM is H + H J = 12 U (cid:88) i a † i a i (cid:16) a † i a i − (cid:17) − µ (cid:88) i a † i a i − (cid:88) ij J ij a † i a j , (1)where nearest neighbor tunneling of strength J ij = J is assumed on a 2D square lattice with T (cid:28) U . Thisbath temperature models the coupling to an environ-ment [19] that dissipates energy. Let us first considerthe T = 0 case where µ/U is close to some integer M sothat µ/U = M ± δ for some 0 < δ < /
2. Letting thetunneling J be infinitesimally small, the ground state isa Mott insulator. The energy gap for adding a parti- a r X i v : . [ c ond - m a t . qu a n t - g a s ] A p r cle or hole is proportional to δ and (1 − δ ) respectively.As δ is tuned to zero (unity), the state with an extraparticle (hole) becomes degenerate with the Mott insu-lator. Thus, excess particles(holes) will be free to hopamong sites with no energy barrier. At low tempera-tures, they will condense producing superfluidity [6]. As J is increased, the low lying excitations become long-range collective particle (hole) tunneling events betweenthe system and reservoir. These events promote particlefluctuations, but they tend to stabilize the phase acrosssites through the number-phase conjugate relationship.The manifold where the energies needed for these long-range excitations vanishes defines the phase boundary [5]at T = 0. When T is finite, sites will have a thermal prob-ability p n = e − βε n of being occupied by n bosons. Siteswill have different energies, and their phases will rotateat different rates. This aggravates the phase fluctuationsthat destroy superfluidity. The region in J vs. µ spacethat permits superfluidity is thus reduced (especially nearinteger values of µ/U ) as the temperature is increased.We will see that the long-range phase-coherence demon-strated by our non-equilibrium correlation functions canbe interpreted physically as the artificial closing of theMott gap thereby allowing for the excitation of long-range tunneling modes.To effect an enhancement of superfluidity, we shallneed a perturbation that pushes the system out ofequilibrium. A thermal bath is also necessary bothto counterbalance the heating induced by the pertur-bation as well as to ensure that the concept of tem-perature remains well-defined. We shall add threeterms to the Hamiltonian so that H ( t ) = H + H J + H bath + H coup + H V ( t ). The bath and cou-pling Hamiltonians H bath = (cid:80) iα ε α b † iα b iα and H coup = (cid:80) iα g α (cid:16) b † iα + b iα (cid:17) a † i a i model the coupling of each siteto an infinite bath of oscillator degrees of freedom such asthe collective modes of a larger condensate in the whichthe lattice is immersed. Since the only crucial functionof the bath is to balance the energy input from the driv-ing, many types of dissipation are possible and the re-sults of this paper are likely to extend to these alter-nate models. While future work is needed to verify thisassertion, low frequency irradiation of a Josephson ar-ray should easily corroborate or deny this intuition ex-perimentally. We shall model the strength of our bathcoupling [19] by a purely local and Ohmic parameter g α = ηε α exp {− ε α / Λ } . The driving term that will forcea departure from equilibrium is given by H V ( t ) = V (cid:88) i a † i a i cos ( k · x i − Ω t ) , (2)In practice, H V ( t ) describes what is called a Bragg pulse[10, 20, 21]. In the limit where energy differences betweeninternal atomic levels are much larger than U , J , and T ,the potential in Eq. (2) is created by the superposition of two lasers offset to each other in both frequency andwave vector. The spatial dependence of our perturbationis necessary for nontrivial results because a constant per-turbation simply multiplies the correlation function by aspatially constant phase factor [15]. This factor is irrele-vant to the tunneling of particles, and it can be gaugedaway by a time-dependent transformation of our opera-tors. The precise form of the spatial dependence in oursystem, however, seems not to be very important andother schemes are certainly possible [2]. We have chosenthis one because of its experimental simplicity.It is difficult to say much about the nature of the inco-herent and coherent nonequilibrium phases or how theyrelate to the Kosterlitz-Thouless transition. Our focuswill be to determine the boundary between the two non-equilibrium phases. To do this we shall approximate thecorrelation function (cid:104) a † i ( t ) a j ( t (cid:48) ) (cid:105) and look for divergencesof this quantity over long distances. We shall do thisperturbatively in the small quantity J/U following thegeneral method in [17, 18]. Anticipating a nonequilib-rium formalism and considering only the first order inthe self-energy, the correlation function can be writtenas an infinite sum of simple chain diagrams defined onthe foward-backward Keldysh contour, C . The evolu-tion on C is important in nonequilibrium problems be-cause it dispenses with the need to know the state of thesystem at t = ∞ for the calculation of expectation val-ues. The contour-time-ordering is accounted for by allow-ing green’s functions to have a matrix structure [12, 14].That is, if we define ˆ G = (cid:18) G A G K G R (cid:19) where A, K, R referto advanced, Keldysh, and retarded Green’s functions,then the non-equilibrium Dyson series can be written asa sum of matrix products of correlation functions.ˆ G ij ( t, t (cid:48) ) = ˆ g ij ( t, t (cid:48) ) − (cid:88) i i (cid:48) J i i (cid:48) (cid:90) ∞−∞ dt ˆ g ii (cid:48) ( t, t ) ˆ G i j ( t , t (cid:48) ) , (3)where ˆ g ij refers to the correlation functions with re-spect to the hamiltonians H + H bath + H coup + H V ( t )at bath temperature T . Because these hamiltoni-ans are just sums of single-site terms proportional toproducts of density operators n i = a † i a i , they areeasy to diagonalize in the occupation basis. Addi-tionally, the bath can be decoupled from the system[14, 19] via a canonical transformation a i → e S a i e − S = a i e (cid:80) α gαεα ( b † iα − b iα ) = a i X i that uses the time-derivativesof the transformed fields to cancel the coupling to thebath. Averages of the form (cid:104) T c a † i ( t ) a j ( t (cid:48) ) (cid:105) simply trans-form to (cid:104) T c a † i ( t ) a j ( t (cid:48) ) (cid:105)(cid:104) T c X † i ( t ) X j ( t (cid:48) ) (cid:105) = ˆ g ij ( t − t (cid:48) ) ◦ ˆ f ij ( t − t (cid:48) ) where ◦ denotes a Hadamard or entrywiseproduct given by (cid:16) ˆ A ◦ ˆ B (cid:17) αβ = A αβ B αβ where α , β des-ignate components in the 2 × H V ( t ) produces a simplephase factor [15] equal to e i V Ω [ sin( k · x i − Ω t ) − sin ( k · x i − Ω t (cid:48) )].We conclude that the non-equilibrium function ˆ g ij isa product of the time-dependent factors produced by H bath + H coup + H V ( t ) and the bare function with re-spect to H . Expanding the non-equilibrium prefactor interms of Bessel functions, the exact expression for ˆ g ij isˆ g ij ( t, t (cid:48) ) = ˆ g bare ij ( t − t (cid:48) ) ◦ ˆ f ( t − t (cid:48) ) × ∞ (cid:88) mn = −∞ ( − m + n J n (cid:18) V Ω (cid:19) × J m (cid:18) − V Ω (cid:19) e i [ n Ω t + m Ω t (cid:48) − ( n + m ) k · x i ] . (4)The functions ˆ f and ˆ g bare ij are 2 × H respectively. As they are equilibriumfunctions, they depend only on the difference of theirtime-arguments. All of the non-equilibrium informationis stored in the expansion of the phase prefactor whichdepends on t and t (cid:48) separately. This is the source ofthe J (cid:0) V Ω (cid:1) dependence of the tunneling renormalizationfamiliar from studies of dynamic localization [2].The non-equilibrium phase factor makes ˆ G ij a functionof τ = ( t + t (cid:48) ) / t − t (cid:48) . How-ever, due to the time-periodicity of Eq. (2), ˆ G ij ( τ, ∆)is a function of τ only up to period 2 π/ Ω. This dis-crete time symmetry of the Hamiltonian allows us to de-compose every matrix function in Eq. (3) as ˆ G ij ( t, t (cid:48) ) = π (cid:80) N e iN Ω τ (cid:82) ∞−∞ e − iω ∆ ˆ G ij ( ω, N ). Following the tech-nique illustrated in [15], equation (3) can now be writ-ten in terms of the functions ˆ g ij ( ω, N ), but it will beburdensome to work with it because the equation forˆ G ij ( ω, N ) will include contributions from ˆ G ij ( ω, N (cid:48) ) forall N (cid:48) . Fortunately, we can mathematically representthis coupling as simple matrix multiplication if we trans-form to the so-called Floquet representation [15, 16] de-fined as ˆ G ij ( ω ) mn = ˆ G ij (cid:0) ω + m + n Ω , m − n (cid:1) . We maynow think of ˆ G ij ( ω ) mn as an infinite square matrix inthe two-dimensional space of Floquet indices m and n .Each element of this matrix is itself a 2 × G i ( ω, q ) mn = (cid:80) j e i q i · ( x j − x i ) ˆ G ij ( ω ) mn . Finally, we explicitly accountfor the plane-wave contribution to the full Green’s func-tion by defining e i ( n − m ) k i · x i ˆ G ( ω, q ) mn = ˆ G i ( ω, q ) mn .Having made these transformations, we arrive at an ex-tremely simple form for the nonequilibrium Dyson equa-tion. ˆ G ( ω, q ) = ˆ g ( ω ) (cid:104) − J ( q , k ) ◦ ˆ G ( ω, q ) (cid:105) , (5)where for given matrices A and B in the Floquet space(indexed by m, n ), AB indicates an ordinary matrixproduct while A ◦ B denotes a Hadamard in the Floquetspace rather than Keldysh space. The matrix J ( q , k ) is a generalized lattice dispersion in Floquet space givenby J mn ( q ) = J (cid:80) ν cos [ q ν + ( m − n ) k ν ] where ν = 1 , G because there are now two types of in-verses corresponding to the two types of products. Thisdouble-product structure in a Dyson equation seems sofar unknown in any other context, but the situation canbe salvaged by column vectorization (CV): the mappingof matrix A to a vector (cid:126)A consisting of the first columnof A stacked on the next column and so on. We can thenmake use of the convenient identities relating Hadamardand ordinary matrix products through CV to solve Eq.(5) and rewrite it as (cid:126) G ( ω, q ) = (cid:16) n p × n p − (cid:2) n p × n p ⊗ ˆ g ( ω ) (cid:3) D (cid:110) (cid:126)J ( q ) (cid:111)(cid:17) − (cid:126)g ( ω ) , (6)where ⊗ indicates a Kronecker product, and D (cid:110) (cid:126)A (cid:111) de-notes the diagonal matrix with entries given by those of (cid:126)A . The identity matrix of size k × k is given by k × k ,while n p is the size of the matrix ˆ g in the Floqet ( m, n )space. It signifies how many higher harmonics we wishto include, or equivalently, the maximum time-resolutionof our treatment. If we needed infinite time-resolution,we would of course let n p → ∞ . However, each off-diagonal element ( m, n ) will be weighted by ( J/U ) m − n while ˆ g mn → m + n , so we expect that n p need not be large to capture the relevant stationarybehavior.To determine the phase boundary, we are only inter-ested in the stationary behavior of the system given byˆ G . Inverting the block-diagonal n p × n p matrix in Eq.(6), the real part of the Keldysh component of our corre-lation function, Re G K , can be displayed. Fig. 1(a) showsthe system at equilibrium ( V = 0) including the effectsof Ohmic dissipation. The phase boundary is given bythe points where Re G K diverges, and our results matchthose of Ref. [19]. Fig. 1(b) is an example of dynamicenhancement of superfluidity. Note the similarity of Fig.1(b) to the T = 0 phase diagram at equilibrium. In Fig.1(c), we see the effect of making the pulse energy ¯ h Ω tentimes smaller. Again there is superfluid enhancement,but notice that the valleys come to much finer pointscloser to the µ/U axis. This advocates an interpreta-tion wherein our perturbation excites phase-stabilizingcollective modes and artificially closes the energy gap. Asmaller perturbation energy ¯ h Ω is resonant with a smallergap. Thus, the phase boundary is much closer to the µ/U axis at integer values of µ/U where the gap goes to zero.It becomes qualitatively identical to the T = 0 equilib-rium phase diagram.The system’s behavior in other regions of V, Ω spaceare catalogued in Fig. 2. For instance, at high frequencies(¯ h Ω (cid:29) U ), we see the familiar phenomenon of dynamic
0. 1. 2.0.0.020.04 Μ (cid:144) U J (cid:144) U (a)
0. 1. 2.0.0.020.04 Μ (cid:144) U J (cid:144) U (b)
0. 1. 2.0.0.020.04 Μ (cid:144) U J (cid:144) U (c) FIG. 1:
Numerical density plots of Re G K for bath temperature βU = 25, coupling width Λ /U = 3, and strength η = 0 .
01. (a)Equilibrium: V = 0. (b) Superfluid enhancement: V/U = 0 . h Ω /U = 0 . k = πa x , n p = 5. (c) Superfluid enhancement: V/U = 0 .
1, ¯ h Ω /U = 0 . k = πa x , n p = 5. Note the similarity tothe T = 0 equilibrium diagram. localization [1, 2]. All Wigner components other than N = 0 can be ignored, and the tunneling is renormalized J → J J ( V /
Ω) by the zeroth Bessel function of the firstkind coming from the expansion of the non-equilibriumphase factor in Eq. (4). If
V /
Ω is tuned to a zero of J ( x ), we have dynamic suppression of tunneling. Wefind a featureless phase diagram (not shown in Fig. 2)because there is no region of J vs. µ space that admitslong-range phase-coherence.This phenomenon, familiar from driven Josephson ar-rays, can be understood as a cancellation of the dynamicphase acquired due to one period of driving with thatof the hopping between sites. The bosons become local-ized with no long-range correlations, and Re G K divergesnowhere. If we now tune V toward zero, the phase bound-ary returns to its equilibrium form. If instead Ω is tunedlower, it will eventually be small enough to be resonantwith the low energy modes available when µ/U is closeto an integer. We will again have dynamic enhancementof superfluidity. However, further from the valley, therewill be no available low-energy modes, and we will seesuppression of tunneling. The result is larger Mott lobesthat almost touch the µ/U axis.We have demonstrated the enhancement of the super-fluid region in parameter space by driving. The exper-imental signature of this effect is similar to what hasbeen found in time-of-flight experiments [22]. When µ and J are tuned to a point within the enhanced super-fluid region close to integer values of µ/U , there will bewell-defined peaks in momentum space when the pertur-bation is on ( V (cid:54) = 0) and a featureless interference pat-tern corresponding to destroyed phase coherence whenthe perturbation is off ( V = 0).The authors are grateful to Ehud Altman, WilliamPhillips, and Hong Ling for illuminating discussions ofthis topic. This research was supported by JQI-PFC(A.R.), DARPA and US-ARO (V.G.), and the PackardFoundation, Sloan Foundation, and NSF grants PHY-0456720 and PHY-0803371 (G.R.) U UV FIG. 2:
The T = 0 . U equilibrium phase boundary (dotted line)plotted with examples from different regions of driving parameterspace. [1] A. Zenesini et al., Laser Physics , 1182 (2010).[2] C. E. Creffield and T. S. Monteiro, Phys. Rev. Lett. ,210403 (2006).[3] C. Weiss and H.-P. Breuer, Phys. Rev. A , 023608(2009).[4] A. Eckardt et al., Phys. Rev. Lett. , 200401 (2005).[5] J. K. Freericks and H. Monien, Euro. Phys. Lett. , 545(1994).[6] M. P. A. Fisher et al., Phys. Rev. B , 546 (1989).[7] A. F. G. Wyatt et al., Phys. Rev. Lett. , 1166 (1966).[8] G. M. Eliashberg, JETP Lett. , 114 (1970); [Pis’maZh. Eksp. Teor. Fiz. , 186 (1970)].[9] M. G. Blamire et al., Phys. Rev. Lett. , 220 (1991).[10] A. Robertson and V. M. Galitski, Phys. Rev. A ,063609 (2009).[11] N. Lindner, G. Refael, and V. Galitski,arXiv:1008.1792v2 (to appear in Nature Physics).[12] A. Kamenev and A. Levchenko, Advances in Physics ,197 (2009).[13] L. V. Keldysh, JETP Lett. , 1018 (1965); [Zh. Eksp.Teor. Fiz. , 1515 (1964)].[14] G. D. Mahan, Many-particle physics (Plenum Press, NewYork, 1981).[15] T. Brandes, Phys. Rev. B , 1213 (1997).[16] N. Tsuji, T. Oka, and H. Aoki, Phys. Rev. B , 235124(2008).[17] W. Metzner, Phys. Rev. B , 8549 (1991).[18] M. Ohliger and A. Pelster, ArXiv e-prints (2008),0810.4399.[19] D. Dalidovich and M. P. Kennett, Phys. Rev. A ,053611 (2009).[20] P. B. Blakie, R. J. Ballagh, and C. W. Gardiner, Phys.Rev. A , 033602 (2002).[21] A. M. Rey et al., Phys. Rev. A , 023407 (2005).[22] M. Greiner et al., Nature415