Tuning towards dynamic freezing using a two-rate protocol

# Tuning towards dynamic freezing using a two-rate protocol

## Abstract

We study periodically driven closed quantum systems where two parameters of the system Hamiltonian are driven with frequencies and . We show that such drives may be used to tune towards dynamics induced freezing where the wavefunction of the state of the system after a drive cycle at time has almost perfect overlap with the initial state. We locate regions in the plane where the freezing is near exact for a class of integrable and a specific non-integrable model. The integrable models that we study encompass Ising and XY models in , Kitaev model in , and Dirac fermions in graphene and atop a topological insulator surface whereas the non-integrable model studied involves the experimentally realized one-dimensional (1D) tilted Bose-Hubbard model in an optical lattice. In addition, we compute the relevant correlation functions of such driven systems and describe their characteristics in the region of plane where the freezing is near-exact. We supplement our numerical analysis with semi-analytic results for integrable driven systems within adiabatic-impulse approximation and discuss experiments which may test our theory.

###### pacs:
73.43.Nq, 05.70.Jk, 64.60.Ht, 75.10.Jm

## I Introduction

The behavior of closed quantum systems driven out of equilibrium has attracted a lot of attention in recent times pol1 (); dziar1 (); dutta1 (). Such out-of-equilibrium systems explore a major region (if not all) of it’s Hilbert space leading to a gamut of novel phenomena which have no equilibrium counterparts. For example, they exhibit universal scaling of excitation density when driven slowly through a quantum critical point kz1 (); pol2 (); ks1 (); ks2 (); pol3 (); sondhi1 (). Another such class of phenomenon involves dynamic transition where the free energy density of a quantum system after a quench develops non-analyticities leading to sharp cusp-like features in it’s Lochsmidt echo pattern at specific times mukh1 (); pol4 (); dytr1 (); dytr2 (); dutta2 (); these transitions have no equilibrium analogs and can not be described by any local order parameter. Furthermore, the universality and the possibility of application of renormalization group (RG) methods to such driven systems has been a subject of recent research ehud1 (); sang1 (); aditi1 (); such analysis also leads to a several interesting features involving RG flow and universality of driven systems which are qualitatively different from their equilibrium counterparts.

A class of such dynamical phenomenon involve periodically driven quantum systems leading to multiple passages through their critical points during a drive cycle stu1 (); dyn1 (); dyn2 (); adas1 (); rev1 (). Such systems exhibit a wide class of dynamical phenomena which do not occur in their aperiodically driven counterparts. For example, they display signatures of Stuckelberg interference phenomenon in their excitations densities, dynamical correlation functions, and thermodynamic quantities such as work distribution dyn2 (); ks3 (). Furthermore, periodically driven integrable systems show a separate class of dynamical transitions which originates from the change in topology of their Floquet spectrum and leaves its imprint on temporal behavior of local correlation functions sen1 (). Such driven systems also exhibit the phenomenon of dynamics induced freezing adas1 (); pekker1 (); uma1 (); rev1 (); this phenomenon manifests itself in a near unity overlap of the system wavefunction after single or multiple drive period(s) with the initial wavefunction at . Such near-unity overlap signifying freezing occurs at specific discrete frequencies and is therefore qualitatively different from the trivial freezing that occurs for ultrafast drive frequencies.

In all of the phenomena mentioned above, the drive protocols used involve quench, ramp, or periodic change of a single parameter of the system Hamiltonian as a function of time. More recently, a separate class of protocols involving linear or power law ramp of two Hamiltonian parameters with distinct rate has been studied in details jay1 (); rev1 (); uma1 (). The novel features of dynamics induced by such two-rate ramp protocols include a generalization of the Kibble-Zurek scaling laws and possibility of reduction of excitation production as the system is ramped through a quantum critical point. The latter feature leads to better fidelity for quantum state preparation and provide a controlled way of reduction of heating during dynamics in critical systems rev1 (). To the best of our knowledge, the study of dynamics of closed quantum critical systems for such two rate protocols has been carried out for ramp protocols only; no such studies have been undertaken for periodic drives.

In this work, we aim to fill up this gap in the literature by studying the effect of two-rate periodic drive protocols on both integrable and non-integrable closed quantum systems. The drive protocol chosen involves periodic variation of two parameters of the system Hamiltonian with frequencies and . For integrable models, we choose a class of models represented by free fermionic Hamiltonian , where is the two component fermion field, denotes fermionic annihilation operator, and is a matrix Hamiltonian density given by

 H→k = [λ1(ω1t)−b→k]τ3+λ2(rω1t)Δ→kτ1 (1)

where and are Pauli matrices in particle-hole space. Such free fermionic Hamiltonians represents Ising and XY models in subir0 () and Kitaev model in kit1 (); sen2 (); feng1 (). In addition, it also describes Dirac-like quasiparticles in graphene grarev1 () and on the surface of a topological insulators toprev1 (). In what follows we shall carry out numerical analysis of this model in the context of XY model in ; we note, however, that our results shall be valid for any other representations of .

The XY model in a transverse field constitute a model for of half-integer spins on a one-dimensional (1D) chain whose Hamiltonian is given by

 HXY=∑⟨ij⟩,α=x,yJαSαiSαj−h∑iSzi, (2)

where are nearest coupling between and the components of the spins, indicates that is one of the neighboring sites of , and is the transverse field. A standard Jordan-Wigner transformation dutta1 () maps Eq. 2 to Eq. 1 in with the identification

 bk = (Jx+Jy)cos(k),λ1=−h Δk = −isin(k),λ2=(Jx−Jy) (3)

For non-integrable system, we choose the recently experimentally realized ultracold boson in their Mott insulating phase in the presence of an artificial electric field created either by shifting the center of the confining harmonic trap or by applying the spatially varying Zeeman magnetic field to the spinor bosons subir1 (); subir2 (); ksen1 (); ksen2 (); uma1 (); bakr1 (). The Hamiltonian of such bosons can be written as

 H1 = −J∑⟨ij⟩(b†ibj+h.c)+∑i[U2ni(ni−1)−Enii]

where denotes boson annihilation operator on site , is the on-site interaction, is the magnitude of the effective electric field in units of energy subir1 (), is hopping amplitude of bosons between the nearest site, and is the boson number operator. It has been shown in Ref. subir1, that the low-energy properties of is captured by a dipole model whose Hamiltonian is given by

 Hd=−J√n0(n0+1)∑ℓ(f†ℓ+fℓ)+(U−E)∑ℓnℓ, (5)

where is the dipole annihilation operator living on the link between sites and , and denotes the boson occupation number of the parent Mott phase. These dipoles constitute bound states of a boson and a hole on adjacent sites and of the lattice and is the dipole number operator residing on the link between these sites. The dipole operators satisfy two additional constraints; first, at most one dipole can occur on any link () and second, the maximal number of dipoles on two adjacent links is also unity (). The model has two ground states; for , the state correspond to the dipole vacuum which is same as the parent Mott state of the bosons while for , the state corresponds to maximal number of dipoles which is a density wave state for the bosons. It was shown in Ref. subir1, that these two ground states are separated by a transition at which belongs to the Ising universality class. These phases of the model and the signature of the Ising transition between them have been verified in a recent experiment bakr1 (). The non-equilibrium dynamics of the model for single parameter quench, ramp and periodic protocols have also been studied ksen1 (); ksen2 (); uma1 (). In what follows, we shall vary the electric field and hopping amplitude of these models according to periodic protocols characterized by frequencies and respectively.

The main results that we obtain from such a study are the following. First, we develop an analytical framework for studying dynamics for integrable models (Eq. 1) subjected to arbitrary two-rate periodic protocol characterized by frequencies and using an adiabatic-impulse approximation nori1 (); the formalism developed constitutes a generalization of the adiabatic-impulse approximation method to two-rate periodic protocols. We show that in the low-frequency regime, where this approximation is expected to be accurate, our results match those obtained using exact numerics. Second, we chart out, both numerically and using adiabatic impulse approximation, the overlap between the wavefunctions at the end of a drive period ) and the initial ground state wavefunction as a function of and . Our results demonstrate that the second drive frequency can be used as a tuning parameter to obtain dynamics induced freezing, i.e., a near-perfect overlap between these wavefunctions. In particular, for the integrable models studied here, we show that seems to be the best choice for obtaining dynamic freezing and provide a semi-analytic explanation of this result. Third, we compute equal-time correlation functions of the XY model as a function of and and discuss the signature of dynamics induced freezing in such correlation functions. Fourth, we study the dynamics of the tilted Bose-Hubbard model driven out of equilibrium via the two-rate protocol described above using exact diagonalization (ED) for finite size systems with . We show that such a drive leads to near-perfect dynamics induced freezing for several for , provide a qualitative explanation of this phenomenon using a numerically supported effective two state model, and compute the dipole-dipole correlation functions which bear the signature of such freezing phenomenon. Finally, we chart our experiments which can be used to test our theory.

The plan of the rest of the paper is as follows. In Sec. II, we analyze the integrable fermionic Hamiltonian (Eq. 1) using both adiabatic-impulse approximation and exact numerics. This is followed by Sec. III, where we study the dynamics of the dipole model (Eq. 5). Finally we summarize our results, discuss their experimental implications, and conclude in Sec. IV.

## Ii Integrable models

In this section we shall study the dynamics of the free fermion Hamiltonian given by Eq. 1 subjected to the two rate periodic protocol. We first analyze such dynamics using an adiabatic-impulse approximation in Sec. II.1. This is followed by discussion of results from exact numerics in Sec. II.2, carried out using appropriate parameters for 1D XY model (Eqs. 2 and 3), and their comparison with those obtained from adiabatic-impulse approximation developed in Sec. II.1.

The basic assumption behind the adiabatic-impulse approximation lies in dividing the dynamics of a system subjected to a drive into two distinct regimes nori1 (). The first is the adiabatic regime where the rate at which the system Hamiltonian changes is small compared to the instantaneous energy gap; here the dynamics merely lead to accumulation of a phase of the system wavefunction. The second constitute the impulse regime where the rate of change of the Hamiltonian parameter is comparable to or larger than the instantaneous energy gap; in this regime, excitations are produced since the system can no longer follow the instantaneous ground state. In the context of the Hamiltonian given by Eq. 1, the latter region occurs when the system reaches the critical point. This approximation therefore becomes accurate for low-frequency drives where the impulse regime occurs near the critical point of the Hamiltonian. In what follows, we shall proceed with the assumption that (Eq. 1) has an isolated critical point which the system crosses during the evolution; this assumption fails when the has a gapless region such as for the Kitaev model kit1 (); feng1 (); sen1 (). In this case, the adiabatic-impulse approximation can not be used in its present form. In this section, we shall only discuss cases where , given by Eq. 1, does not have such gapless regions.

To treat the dynamics of Eq. 1 subjected to a two-rate periodic protocol using the adiabatic-impulse approximation, we first sketch the adiabatic and the impulse regions in Fig. 1. We denote the instantaneous ground state at a given momentum to be the ground sate of the (Eq. 1): , where

 u0→k(t) = λ2(t)Δ→kD→k(t),v0→k(t)=−(λ1(t)−b→k)+E→k(t)D→k(t) E→k(t) = √(λ1(t)−b→k)2+λ22(t)Δ2→k D→k(t) = √(λ1(t)−b→k+E→k(t))2+λ22(t)Δ2→k (6)

The corresponding instantaneous excited state wavefunction for any given is given by . Using Eq. 6, one can write down the wavefunction of the driven system at any time and for a given as

 |ψ→k(t)⟩ = b′1→k(t)|ψ→k(t)⟩1+b′2→k(t)|ψ→k(t)⟩2 (7)

with . In what follows, we shall use the notation and .

To obtain expressions for , we first identify the times and (Fig. 1) at which the system comes to an avoided level crossing for a given momenta . To this end, we note that the instantaneous energy gap of the system is given by

 δE→k = 2√(λ1(ω1t)−b→k)2+λ22(rω1t)Δ2→k. (8)

This gap reaches a minima at at an avoided level crossing for which . This leads to

 ˙λ1(ω1t1→k)[λ1(ω1t1→k)−b→k]+˙λ2(ω2t1→k)λ2(ω2t1→k)Δ2→k = 0,

where . For simple enough protocols, it is possible to solve this equation to obtain an analytic expression for . For example if , one obtains and . In other, more complicated cases, one needs to solve Eq. II.1 numerically to obtain .

Let us now consider the evolution of the system in region I (Fig.  1) for . In this region, the system is in the adiabatic regime and the wavefunction merely gathers a kinetic phase during evolution. The evolution operator in this regime is thus given by nori1 ()

 U1(t,0) = e−iΩ1→k(t)τ3,Ω1→k(t)=12∫t0δE→k(t′)dt′ (10)

Similarly, in region II for and III, for , the evolution operators are given by

 U2(t,t1→k) = e−iΩ2→k(t)τ3,Ω2→k(t)=12∫tt1→kδE→k(t′)dt′ U3(t,t2→k) = e−iΩ3→k(t)τ3,Ω3→k(t)=12∫tt2→kδE→k(t′)dt′

Around , the system enters the impulse regions and excitations are produced. The key aspect of the adiabatic-impulse approximation is to approximate the impulse region to be exactly at . Near these points, the effective Hamiltonian has the form

 Ha→k(t) = [Λ1→k+˙λ1(ω1ta→k)δta→k]τ3 +[Λ2→k+˙λ2(ω2ta→k)δta→kΔ→k]τ1 Λ1→k = λ1(ω1ta→k)−b→k Λ2→k = Δ→kλ2(ω2ta→k),δta→k=t−ta→k, (12)

where for the two impulse regions. We note that in contrast to standard applications nori1 (); ks3 (), the impulse Hamiltonian for the two-rate protocol is not of the Landau-Zener form; in contrast, it has time-dependence both in the diagonal and off-diagonal terms. To apply the Landau-Zener result, we therefore cast it in the standard form following the method suggested in Ref. jay1, . To this end, we first write in terms of a new set of Pauli matrices as

 Ha→k(t) = μa1→k(δta→k−t′a→k)~τ3+μa2→k~τ1 (13)

A comparison of Eqs. 12 and 13 shows that for both equations to be simultaneously correct, one needs

 ~τ3μa1→k = ˙λ1(ω1ta→k)τ3+˙λ2(ω2ta→k)Δ→kτ1 ~τ1μa2→k = (Λ1→k+˙λ1(ω1ta→k)t′a→k)τ3 (14) +(Λ2→k+˙λ2(ω2ta→k)t′a→kΔ→k)τ1

The identities and , then leads to

 μa1→k = √˙λ21(ω1ta→k)+˙λ22(ω2ta→k)Δ2→k μa2→k = [(Λ1→k+˙λ1(ω1ta→k)t′a→k)2 +(Λ2→k+˙λ2(ω2ta→k)t′a→kΔ→k)2]1/2 t′a→k = −(˙λ1(ω1ta→k)Λ1→k+˙λ2(ω2ta→k)Λ2→k)/μ21→k. (15)

For the Hamiltonian given by Eqs. 13 and 15, the probability of excitation production can be read off as pol1 ()

 pa→k = exp[−2πδa→k],δa→k=|μa2→k|2/(2|μa1→k|) (16)

Using Eq. 16, one can write down the unitary evolution matrices which propagates the wavefunction across nori1 (). The elements of can be expressed in terms of and the Stokes phase as

 Σa = ⎛⎜ ⎜⎝√1−pa→ke−i~ϕas→k−√pa→k√pa→k√1−pa→kei~ϕas→k⎞⎟ ⎟⎠, (17) ~ϕas→k = −π4+δa→k(lnδa→k−1)+ArgΓ(1−iδa→k)

where denotes the gamma function. Using Eqs. 10, II.1, and 17, we finally obtain

 (b′1→kb′2→k) = U3(T,t2→k)Σ2→kU2(t2→k,t1→k)Σ1→kU1(t1→k,0)(10) = (q+→k−q∗−→kq−→kq∗+→k)(10)=(q+→kq−→k) q+→k = e−iΩ3→k[√(1−p1→k)(1−p2→k)e−i∑2a=1(Ωa→k+~ϕas→k) −√p1→kp2→kei(Ω2→k−Ω1→k)]
 q−→k = eiΩ3→k[√(1−p2→k)p1→kei(~ϕ2s→k+Ω2→k−Ω1→k) (18) +√(1−p1→k)p2→ke−i(~ϕ1s→k+Ω1→k+Ω2→k)] |ψ→k(T)⟩ = (q+→ku0→k−q−→kv0→kq+→kv0→k+q−→ku0→k)=(u→k(T)v→k(T))

Eq. 18 is the central result of this section. It allows us to express the final wavefunction at the end of a drive period for an arbitrary two-rate drive protocol. We note that the expressions of obtained in Eq. 18 constitutes a generalization of its single-rate protocol counterpart nori1 (). Furthermore, it allows us to obtain semi-analytic expressions for excitation density and all fermionic correlators at the end of the drive. In the context of XY model in a transverse field where the transverse spin correlators can be expressed in terms of fermionic correlators subir1 (), Eq. 18 allows us to obtain analytic expressions for the transverse magnetization and other relevant spin correlation functions. In what follows, we shall be interested in equal-time fermionic correlation functions

 C→i−→j = L−d∑→k⟨ψ→k(T)|(c†→kc→k+h.c.)|ψ→k(T)⟩ (19) ×cos[→k⋅(→i−→j)] = L−d∑→k|v→k(T)|2cos[→k⋅(→i−→j)] F→i−→j = L−d∑→k⟨ψ→k(T)|(ckc−k+h.c.)|ψ→k(T)⟩ ×sin[→k⋅(→i−→j)] = L−d∑→k(u∗→k(T)v→k(T)+h.c.)sin[→k⋅(→i−→j)]

and the excitation density

 nd=1−F,F=L−d∑→k|⟨ψ0→k|ψ→k(T)⟩|2 (20)

where is the wavefunction overlap at the end of the drive which approaches close to unity for dynamics induced freezing, the sum over momenta extends over half of the Brillouin zone, and is the linear dimension of the system. These quantities can be expressed in terms of and the instantaneous eignstates at , and , using Eqs. 18 and 7, as

 nd = 1−F=L−d∑→k|q−→k|2 (21) C→i−→j = L−d∑→k|q−→ku0→k+q+→kv0→k|2cos[→k⋅(→i−→j)] F→i−→j = L−d∑→k[(q∗+→ku0→k−q∗−→kv0→k)(q−→ku0→k+q+→kv0→k) (22) +h.c.]sin[→k⋅(→i−→j)]

In the next section, we shall use Eq. 22 to discuss the behavior of magnetization, excitation density and correlation functions of the transverse field XY models as computed using adiabatic-impulse technique and compare with the exact numerics.

### ii.2 Results for integrable models

In this section, we discuss numerical results obtained by exact numerical solution of the driven Hamiltonian Eq. 1 and compare those obtained Sec. II.1 using generalized adiabatic-impulse approximation. To do this, we take advantage of the fact that represents a two-level Hamiltonian for each pair. Thus for every given , it is possible to write a wavefunction which satisfies

 i∂tu→k(t) = (λ1(t)−b→k)u→k(t)+λ2(t)Δ→kv→k(t) i∂tv→k(t) = −(λ1(t)−b→k)v→k(t)+λ2(t)Δ→ku→k(t), (23)

where we first choose the parameters , , and to be those for XY model in a transverse field (Eq. 3) with , , and . Note that our choice of parameters is equivalent to scaling all energy parameters of the Hamiltonian with and all time scales with . A numerical solution of Eq. 23 enables us to obtain which leads to the defect density and the correlation functions and using Eqs. 20 and 19 at the end of a single drive period .

The results of the numerical procedure described above is shown in Figs. 2, 3, and 4. In Fig. 2(a), we plot as a function of for while in Fig. 2(b), we plot as a function of for . These plots clearly demonstrate that approaches zero close to and ; at this point, has almost perfect overlap with for all which signifies dynamic freezing. We note that such a freezing does not occur for a single parameter drive protocol (which corresponds to ) for any ; thus our results demonstrate that the second drive frequency, or equivalently , for a two-parameter drive protocol can act as a tuning parameter to obtain such freezing. A plot of as a function of the drive amplitude (Fig 3(a)) and the initial phase (Fig 3(b)) for several representative values of indicates that indeed reaches a minimum of for all ; decreases rapidly with increasing for while it displays a slow oscillatory decay for other values of . Further reaches a minimal value for as shown in Fig. 3(b). A comparison of these exact numerical results with those obtained using adiabatic-impulse approximation (Eq. 21) for is shown in Fig. 4. We find that the adiabatic-impulse approximation turns out to match the exact numerical results for for all and for any .

From Figs. 2, 3, and 4, it becomes apparent that the dynamics of is qualitatively different from that at since it shows a rapid decay of with both and . We now present a reason for this behavior. We first note that at , the dynamics is controlled by given by

 Hr=1→k=(Af(ω1t)−b→k)τ3+Δ→kf(ω1t)τ1, (24)

where we have chosen for Figs. 2..4. It turns out that can be cast onto a form jay1 ()

 Hr=1→k=α→k(f(ω1t)−t1→k)~τ3+α2→k~τ1 (25)

where are new Pauli matrices which are related to by

 α1→k~τ3 = Aτ3+Δ→kτ1 α2→k~τ1 = (At1→k−b→k)τ3+Δ→kt1→kτ1. (26)

Using the properties that and (where denotes anticommutator) one obtains

 α1→k = −√A2+Δ2→k t1→k = Ab→k/α21→k,α2→k=−Δ→kb→k/α1→k (27)

We note that such a transformation allows us to write in terms of new Pauli matrices and with all time dependence being shifted to the diagonal term; such transformation exist for non-linear drive protocols only for , i.e., when the off-diagonal and the diagonal terms of are driven by identical drive protocol. The avoided level crossings occur at and where . Linearizing around these point, we find that within adiabatic-impulse approximation, the Landau-Zener transition probability can be read-off from Eq. 25 to be

 p→k = e−π|α2→k|2/|α1→k˙f(ω1t0→k)| (28)

We note that the maximal excitation production occurs from modes around . Assuming that around and , we find that for

 p→k≃e−b2→k|→k−→k0|2z/(ω1A3|√1−b2→k/A2|) (29)

From Eq. 29 we find that as or increases; thus the excitation probability after one cycle of the drive rapidly decays to zero with increasing or . It is instructive to compare this decay to the case where the off-diagonal term in time-independent. In this case one can read off using the standard results provided in Ref. nori1, : . Thus we find approaches unity with increasing drive amplitude much more rapidly for which explains the almost monotonic decay of and hence for .

Finally we plot the change in transverse magnetization at , for the XY model as a function of and in Fig. 5. We note that this magnetization can be easily obtained from the expressions of the correlation function given in Eq. 19

 mz(T) = 2L−d∑→k|v→k(T)|2−1=2C0(T)−1. (30)

We note vanishes if dynamic freezing occurs; such a vanishing of is clearly seen along the line in Fig. 5.

## Iii Titled Bose-Hubbard model

In this section, we shall study the two-rate periodic drive protocol for the tilted Bose-Hubbard model whose Hamiltonian is given by Eq. LABEL:tiltham1. We shall analyze the dynamics of this model for ; in this regime, the dipole Hamiltonian given by Eq. 5 may be used to describe the dynamics of the bosons. In what follows, we shall make the electric field and the hopping time dependent according to the protocol

 E(t) = U−E0cos(ω1t) J(t) = J0(1+cos(rω1t)) (31)

with . The electric field may be made time-dependent by varying the Zeeman field used to generate it as a function of time bakr1 () while the detailed method creating a time dependent in realistic experimental situation has been discussed in Ref. rs1, . Following the convention used in previous studies on this model, we shall choose , and so that the dynamics takes the system through the critical point located subir1 ().