Generation of atypical hopping and interactions by kinetic driving

TL;DR

This paper demonstrates how periodic modulation of hopping in a Bose-Hubbard model can engineer unconventional Hamiltonians, leading to exotic many-body states like a superfluid with opposite momentum condensates.

Contribution

It introduces a method of Floquet engineering that suppresses single-particle hopping while enabling higher-order processes, creating novel Hamiltonians and phases.

Findings

Suppression of nearest-neighbor hopping via periodic driving.

Emergence of a superfluid with opposite nonzero momenta.

Transition from Mott insulator to exotic superfluid phase.

Abstract

We study the effect of time-periodically varying the hopping amplitude in a one-dimensional Bose-Hubbard model, such that its time-averaged value is zero. Employing Floquet theory, we derive a static effective Hamiltonian in which nearest-neighbor single-particle hopping processes are suppressed, but all even higher-order processes are allowed. Unusual many-body features arise from the combined effect of nonlocal interactions and correlated tunneling. At a critical value of the driving, the system passes from a Mott insulator to a superfluid formed by two quasi-condensates with opposite nonzero momenta. This work shows how driving of the hopping energy provides a novel form of Floquet engineering, which enables atypical Hamiltonians and exotic states of matter to be produced and controlled.

Equations

$$H(t)=\sum _{x=0}^{N-1}\left[\frac{U}{2}{n}_x(n_x-1)-J(t)\left(a_{x+1}^{†}{a}_x+H.c.\right)\right],$$

$$H_{\mathrm{eff}}=\sum _{w,x,y,z=0}^{N-1}{Q}_{wxyz}{a}_w^{†}{a}_x^{†}{a}_y{a}_z,$$

$$Q_{wxyz}=$$

$$M(\kappa ;k_l,k_m,k_n)$$

$$\rho (k)=\sum _{j,l}{e}^{\mathrm{i}k\left(r_j-r_l\right)}{P}_{jl}^{(1)},$$

$$|{\Psi }_{\mathrm{dip}}⟩=\frac{1}{\sqrt{2N}}\sum _{x=0}^{N-1}\left(D_x+D_x^{†}\right)|\mathrm{Mott}⟩,$$

$$|\xi ⟩=\alpha |\mathrm{Mott}⟩+\beta |{\Psi }_{\mathrm{dip}}⟩,$$

$$⟨H_{\mathrm{eff}}⟩=$$

$$\displaystyle\begin{aligned} &\langle{\rm Mott}|H_{\rm eff}|\Psi_{\rm dip}\rangle=2N\sqrt{2}\Big{(}\sum_{x\neq 0}+\sum_{x\neq 0,\,2}\Big{)}Q_{2xx0}\ ,\end{aligned}$$

$$⟨\xi |\rho (k)|\xi ⟩=N\left(1+\frac{4\alpha \beta }{\sqrt{N}}\cos (2k)+\frac{3{\beta }^2}{N}\cos (4k)\right)$$

$$⟨\rho (p)\rho (q)⟩=\sum _{j,k,l,m}{e}^{\mathrm{i}p\left(r_j-r_k\right)}{e}^{\mathrm{i}q\left(r_l-r_m\right)}{P}_{jklm}^{(2)},$$

$$⟨\rho (p)\rho (q)⟩=N\left(N-2\right)+N^2{\delta }_{pq},$$

$$|\mathrm{Fock}⟩=\frac{1}{(N/2)!}{\left(a_{\pi /2}^{†}\right)}^{N/2}{\left(a_{-\pi /2}^{†}\right)}^{N/2}|0⟩,$$

$$|\mathrm{cat}⟩=\frac{1}{\sqrt{2N!}}\left[(a_{\pi /2}^{†}{)}^N+(a_{-\pi /2}^{†}{)}^N\right]|0⟩.$$

$$\Delta (N)=E(N,N+1)-2E(N,N)+E(N,N-1),$$

$$\Delta (N)={\Delta }_{\infty }+d/N+\mathcal{O}(1/N^2),$$

$${\Delta }^{\mathrm{ic}}(N)=E(N,N-1)-2E(N,N-2)+E(N,N-3).$$

$$N(r)=\frac{1}{N}\sum _{x=1}^N⟨n_x{n}_{x+r}⟩-⟨n_x⟩⟨n_{x+r}⟩,$$

$$N(r\to \infty )=\frac{-1}{2K_b{\left(\pi r\right)}^2}+A\frac{\cos \left(2\pi r\right)}{r^{2/K_b}}+\dots .$$

$$\frac{1}{2\pi K_b}={\lim }_{k\to 0}\frac{N(k)}{k}.$$

$$A(\omega ,k)=\sum _n{\left|⟨{\Psi }_n|a_k^{†}|{\Psi }_g⟩\right|}^2\delta (\omega -{\omega }_n).$$

$$F=\cos k_l+\cos k_m-\cos k_n-\cos k_p.$$

$$p,q⟷p-k,q+k,$$

$$F=(\cos p+\cos q)(1-\cos k)+(\sin p-\sin q)\sin k.$$

$$\frac{\pi }{2},\frac{\pi }{2}⟷\frac{\pi }{2}+k,\frac{\pi }{2}-k,$$

$$a_x^{†}=\frac{1}{\sqrt{N}}\sum _{l=0}^{N-1}{e}^{-\mathrm{i}{k}_{l}x}{a}_{k_l}^{†},$$

$$a_{k_l}^{†}=\frac{1}{\sqrt{N}}\sum _{x=0}^{N-1}{e}^{\mathrm{i}{k}_{l}x}{a}_x^{†}.$$

$$\begin{array}{rl}H(t)& =-2J_{\mathrm{ac}}\cos (\omega t)\sum _{l=0}^{N-1}\cos (k_l)a_{k_l}^{†}{a}_{k_l}\\ & +\frac{U}{2N}\sum _{l,m,n,p=0}^{N-1}{\delta }_{k_l+k_m,k_n+k_p}{a}_{k_p}^{†}{a}_{k_n}^{†}{a}_{k_m}{a}_{k_l}.\end{array}$$

$$R(t)=\exp \left(-2\mathrm{i}\frac{J_{\mathrm{ac}}}{\omega }\sin (\omega t)\sum _{k_l=0}^{N-1}\cos \left(k_l\right)a_{k_l}^{†}{a}_{k_l}\right),$$

$$\tilde{H}(t)=\mathrm{i}\dot{R}(t)R^{†}(t)+R(t)H(t)R^{†}(t).$$

Full text

I Introduction

Periodically driving a quantum system provides a convenient tool to manipulate and control its properties by “Floquet engineering”. In this technique, Floquet theory is used to describe the effects of a high-frequency driving in terms of an effective static Hamiltonian Eckardt (2017). This method gives a high degree of control over the parameters of the effective Hamiltonian, which allows these systems to be used as both quantum simulators Feynman (1982); Jaksch and Zoller (2005); Lewenstein et al. (2007); Georgescu, Ashhab, and Nori (2014) and to treat mathematical problems Creffield and Sierra (2015). It is also possible for the effective Hamiltonian to have properties very different from the original model, which enables the realization of systems with exotic properties which cannot be produced in other ways, such as Floquet topological insulators Iomin and Fishman (2000); Lindner, Refael, and Galitski (2011) and time crystals Zhang et al. (2017); Choi et al. (2017).

In principle, any term of the Hamiltonian can be periodically varied to yield an effective model. The earliest forms of driving consisted of applying an external potential that oscillated periodically in time. This has the effect of modifying the tunneling dynamics of the system Dunlap and Kenkre (1986); Grifoni and Hänggi (1998), since the tunneling terms do not commute with the potential. Manipulating the effective tunneling in this way has been used, for example, to produce localization Lignier et al. (2007); Kierig et al. (2008) and to drive the Mott transition Eckardt, Weiss, and Holthaus (2005); Creffield and Monteiro (2006); Zenesini et al. (2009), by setting the effective hopping energy to zero. This form of control has also been used to produce artificial magnetic fluxes by inducing phases on the hopping terms, and tuning them to mimic the required Peierls phases Kolovsky (2011); Creffield and Sols (2013); Struck et al. (2012); Hauke et al. (2012). An alternative method of driving, considered more recently is, to oscillate the interaction term with time Rapp, Deng, and Santos (2012); Di Liberto et al. (2014); Gaul et al. (2011); Meinert et al. (2016). This again produces a modification of the tunneling terms in the system, with the unusual feature that the effective (nearest-neighbor) hopping depends on the occupation of the sites involved (so-called “correlated hopping” or “bond-charge interaction” Dutta et al. (2015)). Correlated hopping models of this type can also be produced by applying a resonant driving potential, which gives the additional possibility of simulating anyon physics Sträter, Srivastava, and Eckardt (2016) by engineering the appropriate occupation-dependent Peierls phases on the tunneling elements.

If we consider a Hamiltonian to be composed of potential, interaction, and kinetic terms, driving of the first two terms has thus already been intensively studied. In this work we consider the remaining possibility: to drive the kinetic term, which for a lattice system corresponds to oscillating the tunneling. While the effect of varying the tunneling with time has been studied previously in cold atom systems Stöferle et al. (2004); Kollath et al. (2006); Dirks et al. (2014); Cardarelli, Greschner, and Santos (2016) in the context of lattice modulation spectroscopy, the oscillation is usually taken to be of small amplitude around a larger constant value. Here, in contrast, we consider the nearest-neighbor hopping to oscillate between positive and negative values with a zero time-average. Although this form of driving may seem rather unusual we discuss possible means of achieving it in Section VI. This “kinetic driving” produces an unusual effective Hamiltonian in which there is no sharp distinction between hopping and interaction terms, and where long-range correlated and assisted hopping processes occur. We consider a one-dimensional system, and show that for small values of the driving parameter the dynamics is frozen: the system forms a Mott insulator. As the driving is increased, however, states containing bound doublon-hole pairs, which we term “dipoles”, begin to play an important role. Eventually the system undergoes a quantum phase transition Sachdev (2001) at which the Mott gap vanishes, and forms a many-body cat state, consisting of the superposition of two quasi-condensates with opposite, non-zero momenta.

II Model

We consider a one-dimensional Bose-Hubbard (BH) model with $N$ sites

[TABLE]

where $a_{x}(a_{x}^{\dagger})$ are the usual bosonic annihilation (creation) operators, $n_{x}$ is the number operator, and $U$ is the Hubbard interaction term. The tunneling amplitude between nearest-neighbor sites, $J(t)$, is taken to be a $T$-periodic function with zero time-average, and we use the specific form $J(t)=J_{\mathrm{ac}}\cos(\omega t)$ 111We have also considered squarewave and triangular driving functions, and find that they give very similar behaviour.. We will take $\hbar=1$, and for convenience we characterize the amplitude of the driving by the dimensionless parameter $\kappa=J_{\mathrm{ac}}/\omega$. The majority of our results are obtained for commensurate filling, with the number of bosons equal to $N$. This gives a well-defined Mott state in the limit $\kappa\rightarrow 0$.

As $H(t)$ is time dependent it does not have a set of static eigenvectors with associated energies. However, its time-periodicity means that it can instead be described in terms of Floquet functions and quasienergies, which solve the Floquet equation $(i\partial_{t}-H(t))|\psi_{j}(t)\rangle=\epsilon_{j}|\psi_{j}(t)\rangle$. To obtain the quasienergy spectrum for a given value of $\kappa$, we evaluate the time evolution operator for one period, $U(T,0)$, and diagonalise it to obtain its eigenvalues $\lambda_{j}$. These allow us to then calculate the quasienergies from the relation $\lambda_{j}=\exp(-iT\epsilon_{j})$. The quasienergies are only defined up to integer multiples of the driving frequency, since $\epsilon_{j}$ and $\epsilon_{j}+n\omega$ are not distinguishable, which produces a Brillouin zone structure in the quasienergy coordinate Holthaus (2016). We choose to set $n=0$, fixing the quasienergies to lie in the “first Brillouin zone”, $[0,\omega)$. As we consider the case of high-frequency driving, $\omega$ is so large (and the Brillouin zone is thus so wide) that the quasienergies do not wrap around the zone edges, and so there is no ambiguity about the quasienergy ordering. Consequently we can simply treat the quasienergy spectrum within this Brillouin zone like a standard energy spectrum, and in particular we can consider the lowest quasienergy to be that of the system’s ground state. We shall see later in Section III that the spectrum of the effective Hamiltonian that we derive, which is naturally well-ordered and possesses a well-defined ground state, indeed reproduces the quasienergy spectrum, confirming that this interpretation is justified.

In Fig. 1 we show the quasienergy spectrum obtained for a 6-site, 6-particle model, driven at a high frequency of $\omega=250U$. For $\kappa=0$ the system is equivalent to a standard BH model with $J=0$. The ground state is thus a perfect Mott state with each site having unit occupation, while the excited states fall into separate flat bands that can be classified by the number and type of multiple site occupations. For example, the set of next highest states have one hole and one doublon (and so have an energy of $U$), the next set has two holes and two doublons, and so on. As $\kappa$ is increased from zero, the degeneracies between the excited states are broken, and the bands begin to fan out and overlap. Eventually for $\kappa>0.42$ the ground state approaches the excited states, but little other structure, such as crossings or avoided crossings, is visible.

III Effective model

As shown in detail in Appendix A and discussed, for example, in Refs. Goldman and Dalibard (2014); Creffield et al. (2016), the procedure for obtaining an effective Hamiltonian consists of Fourier transforming the creation and annihilation operators, making a unitary transformation to the interaction picture, and then averaging the transformed Hamiltonian over one period of the driving. The result, transformed back to real space, constitutes the lowest-order term in the Magnus expansion Bukov, D’Alessio, and Polkovnikov (2015); Eckardt and Anisimovas (2015), and consists of a sum of 4-operator terms

[TABLE]

where the transition amplitudes $Q_{wxyz}$ are given by

[TABLE]

and

[TABLE]

Here ${\cal J}_{0}$ is the zeroth Bessel function, and $k_{u}=2\pi u/N$ denote the lattice momenta. Unlike the case of standard potential driving, in which the renormalized Hamiltonian in the high-frequency limit is described by just a single parameter $J_{\mathrm{eff}}$, the Hamiltonian (2) involves a variety of matrix elements, their number and value depending on the size of the lattice.

From Eqs. (2)-(4) it is clear that the only remaining energy scale in the problem is $U$. Thus changing $U$ in the original Hamiltonian (1) simply results in a rescaling of time, without introducing different behaviour. We have verified this by explicit simulation of the full time-dependent Hamiltonian (1). We have also compared the energy spectrum of $H_{\mathrm{eff}}$ with the quasienergies obtained from the full time-dependent model for several lattice sizes and values of $\kappa$, and find excellent agreement between the two 222Comparing the quasienergies $\epsilon_{j}$ and the eigenenergies of the effective Hamiltonian $E^{\mathrm{eff}}_{j}$ by evaluating $\chi^{2}=\sum|\epsilon_{j}-E^{\mathrm{eff}}_{j}|^{2}$, we find that for a fixed value of $\kappa$ it is bound by $\chi^{2}<5\times 10^{-7}$. The ordering and the values of the quasienergies are thus essentially indistinguishable from the energy levels of the effective Hamiltonian..

We can also see that although the time-averaging removes all nearest-neighbor single-particle hopping, higher-order processes are permitted in $H_{\mathrm{eff}}$, and can be of unusual type, with a variety of elementary processes that blur the distinction between interaction, correlated hopping, and assisted hopping. For example, while $Q_{0000}$ and $Q_{2121}$ are standard density-density interaction terms, $Q_{1100}$ represents hopping of a doublon, $Q_{1120}$ provides another instance of correlated hopping, while $Q_{wxxy}$ represents assisted hopping in which the jump between sites $w$ and $y$ is influenced by the occupation of site $x$. It is important to note that hopping can in principle be arbitrarily long-ranged, as has also been seen in studies of the driven-tunneling Kitaev model Benito et al. (2014). However, for small $\kappa$ only a few hopping processes are significant, as shown in Fig. 2, which makes analyzing the system tractable in this range.

IV Results

IV.1 Momentum density function

In Fig. 3 we show the momentum density function for an eight-site system,

[TABLE]

where $P^{(1)}_{jl}=\langle a_{j}^{\dagger}a_{l}\rangle$ is the one-particle reduced density matrix. For $\kappa=0$ the system is in the Mott state and $\rho(k)$ is completely flat, reflecting the absence of phase coherence. As $\kappa$ is increased, however, peaks develop at $k=\pm\pi/2$, suggesting the formation of condensates at those momenta 333For simplicity of language, here and in the following we employ the term “condensate” although, strictly speaking, in one-dimension only a “quasi-condensate” exists.. To further explore the nature of these peaks, we calculate the eigensystem of the one-particle reduced density matrix. Its eigenstates are known as natural orbitals, and the occupation of each natural orbital is given by its corresponding eigenvalue, a condensate being formed if an orbital is macroscopically occupied. In Fig. 4 we can see that the natural orbitals are equally occupied for $\kappa=0$, just as expected for the Mott state. Comparing Fig. 3 with Fig. 4, the formation of the peaks in the momentum density is clearly correlated with the enhanced occupation of two (degenerate) natural orbitals, indicating the formation of a fragmented condensate. Numerical inspection shows that this pair of degenerate orbitals very much resemble plane waves of momentum $\pm\pi/2$. This contrasts with the behaviour of the standard Mott-superfluid transition (described by the undriven BH model), for which a single condensate develops at $k=0$.

Examining the ground state of $H_{\mathrm{eff}}$ over the displayed range of $\kappa$ reveals that it is strongly dominated by the Mott state, with the next most important contributions coming from “dipole states”. These are states containing a doublon-hole pair separated by a single site

[TABLE]

where $D_{x}=a^{\dagger}_{x}a_{x+2}$. Inspection of the numerical results suggests approximating the ground state by

[TABLE]

where $\alpha,\beta\in\mathbb{R}~{}$ and $|\alpha|>|\beta|$ is expected. Although inspired by the numerical simulation of a finite-size system, we emphasize that the state (7), as well as the following discussion, has a well-defined thermodynamic limit.

The expectation value of (2) in the state (7) is

[TABLE]

which suggests that if ${\rm Re\,}\langle{\rm Mott}|H_{\rm eff}|\Psi_{\rm dip}\rangle>0$, the energy is minimized for $\mathrm{sgn}(\alpha\beta)=-1$. As the Mott and dipole states are only connected by matrix elements of the type $Q_{2xx0}$, we obtain

[TABLE]

and we have verified that this is indeed positive (see middle Fig. 2) over the range of $\kappa$ that we investigate, for chains of up to 90 sites. The negative sign of $\left(\alpha\beta\right)$ also causes the momentum density

[TABLE]

to develop peaks at $k=\pm\pi/2$, in agreement with the numerical results shown in Fig. 3.

The presence of density peaks at nonzero momenta can be further understood physically by noting that the predominantly positive character of the matrix elements. $Q_{yxx,y+2}$ corresponds, in the usual convention of Eq. (1), to second-order hopping with a negative effective mass. This favors the occupation of single-particle states with a phase difference of $\pi$ between next-nearest neighbors, which are states with momentum $\pm\pi/2$. For higher driving strengths, however, we find that this simple picture of Eq. (7) breaks down due to the increasing weight of other configurations in the ground state.

Figure 10 in Appendix B shows the ground-state energy as a function of the driving parameter $\kappa$. We observed that, comparing with the exact calculation, the dipole-approximation gives good results up to $\kappa\simeq 0.3$. For that value of $\kappa$, Fig. 3 shows that the momentum distribution is beginning to peak at momenta $k=\pm\pi/2$. Later in this section we will see that the insulator-superfluid transition occurs at $\kappa_{c}\simeq 0.42$. Thus, although the dipole approximation accounts for a partial condensation of bosons at momenta $\pm\pi/2$ [see Eq. (10)], its range of validity remains confined to the insulator region.

IV.2 Two-particle momentum density

The behaviour of the one-particle reduced density matrix clearly indicates that the system makes a transition to a fragmented condensate. However, to characterize this state more precisely, higher-order correlation functions are required Mueller et al. (2006), in particular the two-particle reduced density-matrix, $P^{(2)}_{jklm}=\langle a_{j}^{\dagger}a_{k}a_{l}^{\dagger}a_{m}\rangle$. In Fig. 5 we show the Fourier transform of this quantity

[TABLE]

for two different values of the driving strength $\kappa$. For $\kappa=0$ the system is a perfect Mott insulator, and we can see from Fig. 5a that $\langle\rho(p)\rho(q)\rangle$ manifests peaks along the line $p=q$, and is flat elsewhere. For the case of a perfect Mott state, $\langle\rho(p)\rho(q)\rangle$ can be calculated analytically, yielding the result

[TABLE]

which for momenta discretised over the first Brillouin zone gives perfect agreement with the numerical results.

Below, in Fig. 5b, we show the result for $\kappa=0.5$, for which the system is well within the regime in which the condensate develops. Unlike the Mott state, only two peaks are visible, located at $(p,q)=(\pi/2,\pi/2)$ and $(p,q)=(-\pi/2,-\pi/2)$. This again highlights the fragmented nature of this condensate; in the standard case of the Mott transition, the superfluid would present just a single peak in $\langle\rho(p)\rho(q)\rangle$, centered on the origin.

In addition to confirming the fragmentation of the condensate, Fig. 5b provides information about the nature of the fragmentation. A Fock state in momentum space, having the approximate form

[TABLE]

would present four peaks in the two-particle momentum density, located at the points $(\pm\pi/2,\pm\pi/2)$. However, as Fig. 5b shows two peaks the ground state is clearly not of this type. The absence of the anti-diagonal peaks in the two-particle momentum density instead suggests that the ground state has a Schrödinger-cat structure, which in its idealized form would be

[TABLE]

This finding has considerable implications for the observation of this state in experiment Di Liberto et al. (2011). While a single measurement of the momentum distribution of a Fock state would yield an equal weight for the two momenta $k=\pm\pi/2$, the cat state would be collapsed to one momentum or the other with equal probability. This unusual behaviour would thus provide a clear experimental signal of the exotic superfluid state that this form of driving will produce.

IV.3 Vanishing of the Mott gap

The Mott phase is characterized by a non-zero energy gap, $\Delta$, for adding or subtracting a particle from the system. Conversely, a superfluid state will be gapless. As a result, evaluating $\Delta$ and locating the point at which $\Delta=0$, provides a direct method of identifying the point at which the phase transition occurs. To calculate $\Delta$ we follow the procedure described in Ref. Sowiński (2012), where it was applied to the standard Bose-Hubbard model. We first obtain the ground-state energy for the case of commensurate filling, when the number of lattice sites, $N$, is equal to the number of bosons. We then define the energy gap to be the difference

[TABLE]

where $E(m,n)$ denotes the ground-state energy of an $m$-site system holding $n$ bosons. In order to evaluate $\Delta$ in the thermodynamic limit, we diagonalise the many-body effective Hamiltonian (2) for several different lattice sizes ($N=5,\dots,8$), and make a least-square fit to the function

[TABLE]

to extract the value of $\Delta$ in the $1/N=0$ limit. Although these systems may appear rather small, Ref. Sowiński (2012) showed that this technique indeed provides a robust and well-controlled means of extrapolating to the thermodynamic limit, even for systems of this size. In Fig. 6a we plot $\Delta(N)$ as a function of $1/N$ for several values of $\kappa$, and indeed see an approximately linear behaviour. We show in Fig. 6b the behaviour of $\Delta_{\infty}$, obtained in this way, as $\kappa$ is varied. Although the deviations from linearity in the fitting procedure introduce rather large error bars, the tendency of $\Delta_{\infty}$ is clear. In the absence of driving ($\kappa=0$), $\Delta_{\infty}=U$ as expected. As $\kappa$ is increased, $\Delta_{\infty}$ reduces, and reaches zero at approximately $\kappa=0.42$. Thereafter $\Delta_{\infty}$ remains zero within the error bars, indicating that the Mott gap has closed, and that the system is in the superfluid regime.

As a further test of the method, we consider the case of an incommensurate filling, focusing on the case of a system of $N$ sites with $N-2$ particles, for which the gap is given by

[TABLE]

Analyzing this system in the same way (computing the cases $N=5,\dots,8$) yields a value for the gap that remains zero (within the error bars) for all driving strengths, indicating that the system is always superfluid. This finding anticipates a result of Luttinger liquid theory Cazalilla et al. (2011), which states that for non-commensurate fillings a one-dimensional system falls into a different universality class than for the commensurate case, and will be superfluid for all non-zero values of $\kappa$. This property is also evident in the behaviour of the Luttinger parameter $K_{b}$ (shown in Fig. 7b), as we will discuss in the next section.

IV.4 Luttinger liquid parameters

The unifying concept for one-dimensional interacting systems, of both bosonic and fermionic type, is the Tomonaga-Luttinger liquid Tomonaga (1950); Luttinger (1963). A system of this type exhibits two important features, namely, that correlation functions follow a power-law decay Haldane (1981), and that the low-energy excitations are collective modes with a linear dispersion. A Tomonaga-Luttinger liquid can thus be exactly specified by just two parameters; the interaction parameter $K_{b}$, in which the asymptotic behaviour of all the correlation functions can be expressed, and the group velocity, $v$, of the collective modes.

IV.4.1 Correlation function

For the system we study, $K_{b}$ can be extracted from the the density-density correlation function Kühner, White, and Monien (2000)

[TABLE]

which asymptotically decays as Cazalilla et al. (2011)

[TABLE]

Rather than directly fitting $N(r)$ to Eq. (19), however, more stable results are obtained by evaluating its Fourier transform, $N(k)$, and estimating $K_{b}$ from its derivative Ejima, Fehske, and Gebhard (2011)

[TABLE]

In Fig. 7a we show $N(k)$ for a range of values of $\kappa$. For large $\kappa$, for which the system is gapless, the results show a clearly linear behaviour for small $k$, indicating that this procedure is valid. As $\kappa$ is reduced, however, $N(k)$ is suppressed and has a quadratic behaviour for small $k$, corresponding to the system being gapped and thus not being describable in terms of a Luttinger liquid.

For small values of $\kappa$ the system is in the Mott insulator regime, and consequently $K_{b}$ should diverge. In Fig. 7b we can see that the calculated values of $K_{b}$ are indeed large, but remain finite since the system is far from the thermodynamic limit. As $\kappa$ is increased the value of $K_{b}$ drops. For any commensurate filling, the transition between the Mott state and the Luttinger liquid will be of Kosterlitz-Thouless type, and at the transition $K_{b}$ will take a universal value. For unit filling this is given by $K_{b}=1/2$, and accordingly we can use this as a criterion for identifying the point of the phase transition Giamarchi (2003); Cazalilla et al. (2011).

We can also note that since Luttinger liquids with $K_{b}<2$ are dominated by superfluid correlations Cazalilla et al. (2011), the transition will be to a superfluid state, exactly as we would expect from the vanishing of the gap seen previously. From Fig. 7b we can see that $K_{b}$ crosses the critical value at approximately $\kappa=0.48$. This is slightly larger than the value $\kappa=0.42$ obtained in the previous section from the analysis of the gap, but still represents good agreement given that the previous result was extrapolated to the thermodynamic limit, while this was obtained on a finite system of 8 sites.

Doping the Mott insulator away from commensurate filling has the effect of immediately producing superfluidity. Accordingly, for an incommensurate filling the system will only be a Mott insulator at $\kappa=0$, and will immediately make a phase transition as soon as $\kappa>0$. At this point $K_{b}$ will again take a universal value, in this case being exactly twice that of the commensurate transition Cazalilla et al. (2011). In Fig. 7b we show the values of $K_{b}$, calculated in the same way, for an incommensurate filling of 6 bosons on an 8-site lattice. As expected, $K_{b}\rightarrow 1$ in the limit of small $\kappa$. The excellent accuracy of this limiting value of $K_{b}$, obtained with no fitting parameters, provides a strong indication that this method of extracting $K_{b}$ is indeed reliable.

IV.4.2 Spectral function

A final check of the Luttinger liquid picture is to confirm that the low-lying excitations of the system have a linear dispersion, yielding a well-defined value for the other Luttinger parameter, the group velocity $v$. To do this we calculate the zero temperature spectral functions for the emission of a single particle with momentum $k$ and energy $\omega$

[TABLE]

Here $|\Psi_{n}\rangle$ is the $n$-th excited state in the $N+1$ particle system, $|\Psi_{g}\rangle$ is the ground state of $N$ particles and $\omega_{n}=E_{n}(N,N+1)-E_{0}(N,N)$, where $E_{n}(N,N+1)$ is the energy of the $n$-th excited state in the $N+1$ particle system. To smooth out this spectral function we convolve Eq. (21) with a narrow Lorentzian. The presence of sharp peaks in $A(\omega,k)$ indicate the presence of well-defined quasiparticle excitations, and plotting the maxima of $A(\omega,k)$ allows us to determine their dispersion relation $\omega(k)$.

In the upper panels of Fig. 8 we show the results obtained for a system well in the superfluid regime ($\kappa=0.53$). For each value of momentum the spectral function shows a single sharp peak, whose location shifts as the momentum changes. The contour plot of $A(\omega,k)$ reveals that this momentum dependence corresponds to two regions of linear dispersion, centered on $k=\pm\pi/2$. We can thus picture the ground state as consisting of two independent Luttinger liquids of bosons, condensed at these two different momenta but otherwise identical, with the collective excitations being density fluctuations about these two condensates 444By independent Luttinger liquids, we mean here two sets of bosons $A$ and $B$, occupying different regions of the one-particle Hilbert space (here around momenta $\pm\pi/2$) such that the total density fluctuation operator can be written as $\delta n_{A}(x)+\delta n_{B}(x)$ while ${\langle\Psi|\delta n_{A}(x)\delta n_{B}(x^{\prime})|\Psi\rangle}=0$.. In Appendix C we argue that $K_{b}=1/2$ remains the (commensurate case) critical value of this double Luttinger liquid system.

The lower panels of Fig. 8 show the corresponding spectral functions when the system is in the Mott insulator regime ($\kappa=0.1$). In contrast to the previous result, the peaks in $A(\omega,k)$ barely show any dependence on momentum, and the flat dispersion relation obtained clearly indicates that the system has a gap of approximately $U$, exactly as would be expected for a Mott insulator.

V Scattering processes

As shown in Figs. 3 and 4, and discussed in Section IV, it is clear that the superfluid is composed of a pair of macroscopically-occupied degenerate orbitals, with momenta $\pm\pi/2$. The existence of this fragmented condensate has been explained on physical grounds by noting that the dominant matrix elements linking the Mott state and the dipole states correspond to single-particle hopping processes with “negative effective mass” between next-nearest neighbor sites. Here we examine the various scattering amplitudes in more detail to try to understand the dominant processes coexisting with the macroscopic occupation of $\pm\pi/2$ momentum eigenstates, and thus to obtain a microscopic picture of the physics occurring.

An important piece of information that we can take from the simulations of this system is that the ground state only involves Fock states of zero total crystal momentum, i.e. the many-body ground state is an eigenstate of the total momentum with eigenvalue zero. Quite generally, particles of initial momenta $k_{l},k_{m}$ collide to final momenta $k_{n},k_{p}$, with $k_{n}+k_{p}=k_{l}+k_{m}$, the connecting matrix elements being given by Eq. (4) [see also Eq. (37)] The Bessel function appearing there can be written as ${\cal J}_{0}\left(2\kappa F\right)$, with

[TABLE]

Our goal here is to elucidate some basic features of the many-body ground state.

First we note that the sheer macroscopic occupation of a given momentum state does not favor a particular value of the privileged momentum, as $F$ is independent of $p$ when $k_{l}=k_{m}=k_{n}=k_{p}=p$. So the preference for a macroscopically occupied momentum will be determined by the ability of the condensate to connect with other configurations. More specifically, the success of a particular condensate choice will depend on how it relates to its depletion cloud. So the situation here is quite different from that of a standard condensate (such as the one yielded by the undriven BH model, for example). There the ground state is known to be dominated at zero temperature by the macroscopic occupation of the momentum $p=0$ combined with pair collisions of the type $0,0\leftrightarrow k,-k$, with $|k|$ small but nonzero. In the standard case, there are energetic reasons for favoring the macroscopic occupation of $p=0$, the depletion cloud merely contributing to the robustness of the resulting superfluid.

By contrast, here the connection between the condensate and its depletion cloud is essential to decide the macroscopically occupied momenta. A reasonable criterion is that the occupation of a given one-particle momentum state is favored in the ground state if it can intervene in processes yielding large matrix elements, as a suitable choice of the relative sign between the intervening many-body configurations permits a strong energy decrease. For moderate $\kappa$ [such that $4\kappa F\lesssim x_{1}$, with $x_{1}=2.404$ the first zero of ${\cal J}_{0}(x)$] this translates into small values of $F$.

A generic pair collision of the type

[TABLE]

has an associated matrix element which is proportional to 555Umklapp processes with final momenta $(p-k,q+k+2\pi)$ yield the same $F$ value.

[TABLE]

For a given momentum pair $(p,q)$, the transferred momentum $k$ spans symmetrically around $k=0$. Thus, on average, when determining the typical value of $F$ for given $(p,q)$, the term in (24) proportional to $\sin k$ will tend to cancel out, while the $\cos k$ term will average to a non-zero value. According to this criterion, the favored momenta will be those which minimize $|\cos p+\cos q|$, which requires $\cos q=-\cos p$. If now we note that the magnitudes of the matrix elements of boson collisions are enhanced by the macroscopic occupation of a single momentum state (bosonic enhancement), we conclude that the case $p=q$ must also be advantageous. These two favorable requirements (zero average $F$ and $p=q$) can be satisfied simultaneously only for $p=q=\pm\pi/2$, in accordance with the numerical findings.

Thus, both from numerical evidence and semiquantitative arguments, the ground state appears to be dominated by configurations connected to each other by processes involving pair collisions of the type

[TABLE]

against the background of a macroscopically occupied $\pi/2$-momentum state, and similarly for the condensate at momentum $-\pi/2$. In both cases, the main processes are of the type (25) with small but nonzero values of $|k|$, all yielding small values of $F$ in (22) and thus large matrix elements.

Therefore the picture is reinforced of two condensates at momenta $\pm\pi/2$ with depletion clouds packed near those two momenta. To the extent that the two depletion clouds can be viewed as peaked around momenta $\pm\pi/2$ and thus non-overlapping, it seems fair to speak of a fragmented condensate, as suggested by our numerical results on the reduced one-particle density matrix (see Fig. 4). This is also consistent with our results for the two-particle reduced density-matrix, shown in Fig. 5.

The thermodynamic limit of Hamiltonian (37) remains to be studied in greater detail. A plausible outcome is that, as the one-dimensional system becomes larger, the relative atom density of both depletion clouds grows logarithmically. Thus the two quasi-condensates shrink in relative size while their respective depletion clouds grow and thus overlap more strongly through pairs of the type shown in the r.h.s. of (25) with $|k|$ comparable to $\pi/2$, the result being a peculiar one-dimensional system.

VI Experimental implementation

As we noted earlier, the form of driving that we choose is rather unusual. In an optical lattice system it is straightforward to modulate the depth of the system to alter the amplitude of the hopping term, but we also require changing its sign so that its time-average vanishes. One method of achieving this would be to make use of the well-known lattice shaking technique, which indeed provides a means of precisely controlling the tunneling in this way. If a system of cold atoms, described by a lattice tunneling model, is periodically accelerated in space at a high shaking frequency $\Omega$, then the tunneling term is renormalized to an effective value. For sinusoidal driving, for example, $J_{\mathrm{eff}}=J{\cal J}_{0}\left(K/\Omega\right)$, where $K$ is the amplitude of the shaking and ${\cal J}_{0}$ is the zeroth Bessel function. We show this behaviour in Fig. 9, and we can see that when the Bessel function is zero (at $K/\Omega\simeq 2.404$) the effective tunneling vanishes.

We can now consider altering the parameters of the shaking Novičenko, Anisimovas, and Juzeliūnas (2017) at a timescale much slower than the period of the shaking $\tau=2\pi/\Omega$. The effective hopping will now vary according to $K(t)/\Omega$. In particular we can slowly oscillate $K$ about the root of the Bessel function at a frequency $\omega$, where $\omega\ll\Omega$, which will give a time-dependent effective tunneling $J_{\mathrm{eff}}(t)$. For a suitable choice of the variation $K(t)$, this will thus produce the time-dependent Hamiltonian (1). In order for Eq. (1) to be in the high-frequency regime itself, which we require for our analysis, we therefore need the driving frequencies to follow the hierarchy: $J\ll\omega\ll\Omega$.

These conditions can be met in systems of periodically-driven cold atoms such as rubidium-87 Lignier et al. (2007); Struck et al. (2012) and potassium-39 Reitter et al. (2017). In order to obtain a single-band BH model as a starting point, it is necessary to use optical lattice depths of $V_{0}>6E_{r}$, where $E_{r}$ is the recoil energy, in order to eliminate all except the nearest neighbor hopping terms. This permits shaking frequencies of several tens of kilohertz to be used without driving atoms into the next highest band, while the bare tunneling, $J$, is of the order 100 Hz. This amply satisfies the first condition, $J\ll\Omega$, and gives a wide window of available frequencies for the frequency $\omega$ at which the tunneling is modulated.

A further consideration is that one might expect a periodically-driven system to absorb energy from the driving, and thus eventually heat up to an “infinite temperature”, a phenomenon known as the eigenstate thermalization hypothesis. If such heating occurs, it would be doubtful that the coherent quantum effects we have studied would persist. However, in the high-frequency limit, when the driving frequency is much larger than the energy scales of the undriven Hamiltonian, it has been shown Abanin, De Roeck, and Huveneers (2015); Abanin et al. (2017) that a generic Floquet system first passes through a long-lived prethermal regime, in which the heating rate is an exponentially small function of the driving frequency. Thermalization will eventually occur, but on a much longer timescale. In the prethermal regime, the Floquet states of the system can be expected to have long lifetimes, and the dynamics of the system can be well-approximated by time evolution under the effective (Floquet) Hamiltonian. As the Bose-Hubbard model is known to exhibit this behaviour Bukov et al. (2015), and the driving we consider is well inside the high-frequency regime, we can be confident that the prethermal regime will be sufficiently large to allow the unusual effects we predict to be experimentally observed.

Standard time-of-flight techniques can be used to image the momentum density function, which should reveal the signatures of the fragmented condensate we predict, shown in Fig. 3. While fragmented condensates can be expected to be unstable to perturbations coupling the two fragments, we have verified that the state we observe is robust to the presence of weak diagonal impurities. Finally we have also verified that the effects we observe are not destroyed by the presence of a residual static hopping, $J_{\mathrm{dc}}$, between lattice sites, so long as $J_{\mathrm{dc}}/U$ is lower than the critical value for the standard Mott transition. Consequently we believe that its experimental realization is completely feasible.

VII Conclusions

We have shown that by driving the tunneling of the BH model with a zero average value, it is possible to produce an effective time-independent Hamiltonian with a number of surprising properties, a remarkable one being the macroscopic occupation, in the ground state, of momentum eigenstates $\pm\pi/2$. While the creation of similar condensates has been been theoretically predicted Rigol and Muramatsu (2004); Rodriguez et al. (2006) and experimentally observed Vidmar et al. (2015) in the context of the expansion of an initially confined Mott state, the effect we report is rather different, being an equilibrium property of the many-body ground state. The sensitivity of the physics discussed here to the presence of impurities and the choice of boundary conditions will be the subject of a future study. The form of Floquet engineering investigated in this paper opens the prospect of a new way of controlling the coherent dynamics of quantum systems, and producing unusual states of matter. Extending this work to higher dimensions and to fermionic systems remain fascinating subjects for future research.

Acknowledgements.

We would like to thank Wolfgang Ketterle and Germán Sierra for valuable discussions. This work has been supported by Spain’s MINECO through Grant Nos. FIS2013-41716-P and FIS2017-84368-P. One of us (FS) would like to acknowledge the support of the Real Colegio Complutense at Harvard and the Harvard-MIT Center for Ultracold Atoms, where part of this work was done.

Appendix A Effective Hamiltonian

To obtain the effective description of the system [Eq. (1)] it is useful to work in the momentum representation. We employ the discrete Fourier transform of the creation and annihilation operators

[TABLE]

where $k_{l}=2\pi l/N$, with the corresponding inverse transform

[TABLE]

The Hamiltonian then becomes

[TABLE]

Transforming into the interaction picture helps to find the non-vanishing time-averaged contributions to the effective Hamiltonian. It is straightforward to see that in Eq. (1) the hopping immediately vanishes when averaged over one period. Using the unitary transformation

[TABLE]

takes Eq. (28) into

[TABLE]

In performing the transformation it is essential to rearrange terms of the form

[TABLE]

Differentiation with respect to $\lambda$ results in

[TABLE]

which implies that

[TABLE]

Using the last equation in combination with Eq. (30) and Eq. (29) gives

[TABLE]

Time averaging then provides us with the effective Hamiltonian

[TABLE]

where $\mathcal{J}_{0}$ is the Bessel function of zeroth order and $\kappa=J_{\rm ac}/\omega$. Using the inverse transformation in Eq. (27) we find in position space

[TABLE]

where

[TABLE]

which is the form of the effective Hamiltonian given in Eqs. (2 -4).

Appendix B Dipole approximation

Using the effective Hamiltonian Eq. (2), we can calculate the expectation value of the energy of the dipole state (7) and then vary its parameters ($\alpha$ and $\beta$) to minimize this energy. This then provides a variational approximation for the ground state of the system. As the normalization of the state requires $\alpha^{2}+\beta^{2}=1$, we thus wish to minimize $\langle H_{\rm eff}\rangle_{\alpha}$, where $\beta=\sqrt{1-\alpha^{2}}$. The variational energy can then be written as

[TABLE]

where

[TABLE]

The energy is minimized for

[TABLE]

The dipole approximation works well for $\kappa<0.25$, as can be seen in Fig. 10.

Even though the ansatz is not too far from the exact energy even for $\kappa=0.4$, it is unable to describe the phase transition since its density-density correlation

[TABLE]

is always short ranged. This means that across the transition, higher-order excitations of the undriven Hamiltonian are essential to produce the long range density correlations, indicating a transition to a superfluid.

Appendix C Two-component Luttinger liquids

Let us assume we have two different Luttinger liquids (call them A and B) that coexist in the same region of space but which are independent in terms of their density fluctuations

[TABLE]

where

[TABLE]

and $\alpha=A,B$. Let us also assume that we can only measure the total density

[TABLE]

Since A and B are both Luttinger liquids, they separately satisfy Eq. (19)

[TABLE]

where

[TABLE]

and the average $\langle\cdot\rangle_{\alpha}$ is taken with respect to a state that only contains species $\alpha$. Here, $L$ is the number of sites and we assume $K_{b}<1$ is identical for the two liquids. We can apply the theory of Cazalilla et al. (2011) and note that, for the commensurate case, $K_{b}=1/2$ signals the onset of superfluidity for both A and B. Importantly, $K_{b}$ only depends on the commensurability of particles and sites and not on the average density. We expect the total density fluctuations to follow a similar law,

[TABLE]

and address the question of how $Q$ relates to $K_{b}$. Eqs. (46) - (48) imply

[TABLE]

where one assumes that the joint system of A and B is in a state such that Eq. (54) is satisfied. Using Eq. (50) one can readily conclude that $K_{b}=Q$.

Now let us see whether the above assumptions apply to a many-body state of cat type

[TABLE]

with $|\Phi_{A}\rangle\equiv A^{\dagger}|{\rm vac}\rangle$ and $|\Phi_{B}\rangle\equiv B^{\dagger}|{\rm vac}\rangle$ being orthonormal $\langle\Phi_{\alpha}|\Phi_{\alpha^{\prime}}\rangle=\delta_{\alpha,\alpha^{\prime}}$. Since $A^{\dagger}$ and $B^{\dagger}$ are, by construction, operators that only create A- and B-type bosons, respectively, we can write

[TABLE]

Clearly the cat state (55) satisfies (54), with $\langle O\rangle_{\alpha}=\langle\Phi_{\alpha}|O|\Phi_{\alpha}\rangle$. Thus we infer that the result $Q=K_{b}$ also applies to the above cat state.

It only remains to analyze to what extent our fragmented-condensate boson system conforms to the two Luttinger-liquid picture defined above. For that purpose we focus on whether our fragmented condensate satisfies Eqs. (55) and (56). Numerical evidence discussed in Section IV.2 strongly suggests that our ground state is of the cat type (55), where $A$, $B$ are operators that create many bosons peaked around $\pm\pi/2$, respectively. We now need to justify whether this constitutes two independent boson species in the sense of Eqs. (46) and (48). We focus on the homogeneous case for which the average density is position independent. In general, if we define

[TABLE]

we can write

[TABLE]

where $a_{k}^{\dagger}$ creates a boson of momentum $k$. Given that we are only interested in the long-distance (large $r$) behavior of the correlation functions, we can assume that only small $q$ values will be relevant in (58). Then we can divide the operators (57) and (58) into two contributions, labeled A [with $k\in(0,\pi)$] and B [with $k\in(-\pi,0)$], both restricted to small $q$. Thus, for instance,

[TABLE]

if we correspondingly define $\delta\hat{n}_{A}=\hat{n}_{A}-\frac{1}{L}\sum_{x}\hat{n}_{A}(x)$. This allows us to approximately decompose the total density fluctuation operator as

[TABLE]

Clearly $\hat{n}_{A}(x)$ and $\hat{n}_{B}(x)$ satisfy Eq. (56). This permits us to infer Eq. (46), which is crucial for the derivation of Eq. (53). Thus, the cat-type state in Eq. (55) with $A^{\dagger}$ and $B^{\dagger}$ creating bosons near $\pm\pi/2$ can rightly be viewed as formed by two independent condensates in the sense of Eqs. (46) and (48). This provides a plausible explanation of why our peculiar one-dimensional system, involving two condensates that cannot be distinguished in terms of density, has the same coefficient $Q=K_{b}$ as the standard, one-component Luttinger liquid.