Generation of atypical hopping and interactions by kinetic driving

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.

At the top of page 7, equation (12) should read:

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 [1]. 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 [2][3][4][5] and to treat mathematical problems [6]. 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 [7,8] and time crystals [9,10]. 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 [11,12], 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 [13,14] and to drive the Mott transition [15][16][17], 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 [18][19][20][21]. An alternative method of driving, considered more recently is, to oscillate the interaction term with time [22][23][24][25]. 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' [26]). 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 [27] by engineering the appropriate occupationdependent 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 [28][29][30][31] 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 6. 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 [32] 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.

Model
We consider a one-dimensional Bose-Hubbard (BH) model with N sites a a x x 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 Tperiodic function with zero time-average, and we use the specific form . We will take ÿ=1, and for convenience we characterize the amplitude of the driving by the dimensionless parameter k w = J ac . 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 k  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 . To obtain the quasienergy spectrum for a given value of κ, we evaluate the time evolution operator for one period, U(T,0), and diagonalise it to obtain its eigenvalues λ j . These allow us to then calculate the quasienergies from the relation The quasienergies are only defined up to integer multiples of the driving frequency, since ò j and ò j +nω are not distinguishable, which produces a Brillouin zone structure in the quasienergy coordinate [33]. We choose to set n=0, fixing the quasienergies to lie in the 'first Brillouin zone', [0, ω). As we consider the case of high-frequency driving, ω 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 3 that the spectrum of the effective Hamiltonian that we derive, which is naturally well-ordered and possesses a welldefined ground state, indeed reproduces the quasienergy spectrum, confirming that this interpretation is justified.
In figure 1 we show the quasienergy spectrum obtained for a 6-site, 6-particle model, driven at a highfrequency of ω=250U. For κ=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 κ is increased from zero, the degeneracies between the excited states are broken, and the For small values of κ the quasienergies group into well-defined flat bands separated by the interaction energy U. As κ increases the bands fan out and begin to mix. The vertical dashed line marks the phase transition between the Mott insulator and the superfluid, estimated from the vanishing of the Mott gap, D ¥ , as extrapolated to the thermodynamic limit (see figure 6).
bands begin to fan out and overlap. Eventually for κ>0.42 the ground state approaches the excited states, but little other structure, such as crossings or avoided crossings, is visible.

Effective model
As shown in detail in appendix A and discussed, for example, in [34,35], 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 [36,37], and consists of a sum of 4-operator terms where the transition amplitudes Q wxyz are given by  Here  0 is the zeroth Bessel function, and k u =2π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 eff , the Hamiltonian (2) involves a variety of matrix elements, their number and value depending on the size of the lattice.
From equations (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 behavior. We have verified this by explicit simulation of the full time-dependent Hamiltonian (1). We have also compared the energy spectrum of H eff with the quasienergies obtained from the full time-dependent model for several lattice sizes and values of κ, and find excellent agreement between the two 3 .
We can also see that although the time-averaging removes all nearest-neighbor single-particle hopping, higher-order processes are permitted in H 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 [38]. However, for small κ only a few hopping processes are significant, as shown in figure 2, which makes analyzing the system tractable in this range.

Momentum density function
In figure 3 we show the momentum density function for an eight-site system, is the one-particle reduced density matrix. For κ=0 the system is in the Mott state and ρ(k) is completely flat, reflecting the absence of phase coherence. As κ is increased, however, peaks develop at k=±π/2, suggesting the formation of condensates at those momenta 4 . 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 figure 4 we can see that the natural orbitals are equally occupied for κ=0, just as expected for the Mott state. Comparing figure 3 with figure 4, the formation of the 3 Comparing the quasienergies ò j and the eigenenergies of the effective Hamiltonian E j eff by evaluating  c = å -| | E j j 2 eff 2 , we find that for a fixed value of κ it is bound by χ 2 <5×10 −7 . The ordering and the values of the quasienergies are thus essentially indistinguishable from the energy levels of the effective Hamiltonian. 4 For simplicity of language, here and in the following we employ the term 'condensate' although, strictly speaking, in one-dimension only a 'quasi-condensate' exists.  . Momentum density ρ(k) for different values of κ, for a system of 8 particles on 8 sites. When κ=0, the system is in the Mott state and ρ(k) is a constant. As κ is increased, the momentum density function develops peaks at k=±π/2, indicating the formation of condensate states with non-zero momentum. In this plot we smoothly interpolate between the 8 allowed values of lattice momentum. 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 ±π/2. This contrasts with the behavior 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 eff over the displayed range of κ 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 Inspection of the numerical results suggests approximating the ground state by , and a b > | | | |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 , the energy is minimized for ab = -( ) sgn 1. As the Mott and dipole states are only connected by matrix elements of the type Q xx 2 0 , we obtain and we have verified that this is indeed positive (see middle figure 2) over the range of κ that we investigate, for chains of up to 90 sites. The negative sign of ab ( ) also causes the momentum density to develop peaks at k=±π/2, in agreement with the numerical results shown in figure 3. The presence of density peaks at non-zero 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 equation (1), to second-order hopping with a negative effective mass. This favors the occupation of singleparticle states with a phase difference of π between next-nearest neighbors, which are states with momentum ±π/2. For higher driving strengths, however, we find that this simple picture of equation (7) breaks down due to the increasing weight of other configurations in the ground state. Figure B1 in appendix B shows the ground state energy as a function of the driving parameter κ. We observed that, comparing with the exact calculation, the dipole approximation gives good results up to κ;0.3. For that value of κ, figure 3 shows that the momentum distribution is beginning to peak at momenta k=±π/2. Later in this section we will see that the insulator-superfluid transition occurs at κ c ;0.42. Thus, although the dipole approximation accounts for a partial condensation of bosons at momenta ±π/2 (see equation (10)), its range of validity remains confined to the insulator region. . Eigenvalues of the one-particle reduced density matrix. At κ=0 all the eigenvalues are equal, implying the even spread of the Mott state over all orbitals. As κ increases, two (degenerate) eigenvalues increase at the expense of the others, signifying the macroscopic occupation of these two natural orbitals, and thus the formation of a fragmented condensate. The degeneracy of each set of eigenvalues is labeled in the figure.

Two-particle momentum density
The behavior 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 [39], in particular the two-particle reduced density matrix, = á ñ , for two values of the driving parameter κ for an 8-site system. Above: κ=0. The system is a perfect Mott insulator, and ρ (2) (p, q) consists of a uniform background, with a line of peaks along p=q. Below: κ=0.5. ρ (2) (p, q) now shows two peaks, centered on (−π/2, −π/2) and (π/2, π/2). This corresponds to the system being in a cat state; a superposition of N particles peaked at momentum −π/2, and N particles peaked at π/2.
which for momenta discretised over the first Brillouin zone gives perfect agreement with the numerical results.
Below, in figure 5(b), we show the result for κ=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)=(π/2, π/2) and (p, q)=(−π/2,−π/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 r r á ñ ( ) ( ) p q , centered on the origin. In addition to confirming the fragmentation of the condensate, figure 5(b) provides information about the nature of the fragmentation. A Fock state in momentum space, having the approximate form would present four peaks in the two-particle momentum density, located at the points (±π/2, ±π/2). However, as figure 5(b) 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 This finding has considerable implications for the observation of this state in experiment [40]. While a single measurement of the momentum distribution of a Fock state would yield an equal weight for the two momenta k=±π/2, the cat state would be collapsed to one momentum or the other with equal probability. This unusual behavior would thus provide a clear experimental signal of the exotic superfluid state that this form of driving will produce.

Vanishing of the Mott gap
The Mott phase is characterized by a non-zero energy gap, Δ, for adding or subtracting a particle from the system. Conversely, a superfluid state will be gapless. As a result, evaluating Δ and locating the point at which Δ=0, provides a direct method of identifying the point at which the phase transition occurs. To calculate Δ we follow the procedure described in [41], where it was applied to the standard BH 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 where E(m, n) denotes the ground state energy of an m-site system holding n bosons. In order to evaluate Δ in the thermodynamic limit, we diagonalise the many-body effective Hamiltonian (2) for several different lattice sizes (N=5, K, 8), and make a least-square fit to the function , 16 2 to extract the value of Δ in the 1/N=0 limit. Although these systems may appear rather small, [41] 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 figure 6(a) we plot Δ(N) as a function of 1/N for several values of κ, and indeed see an approximately linear behavior. We show in figure 6(b) the behavior of D ¥ , obtained in this way, as κ is varied. Although the deviations from linearity in the fitting procedure introduce rather large error bars, the tendency of D ¥ is clear. In the absence of driving (κ=0), D = ¥ U as expected. As κ is increased, D ¥ reduces, and reaches zero at approximately κ=0.42. Thereafter D ¥ 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 Analyzing this system in the same way (computing the cases N=5, K, 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 [42], which states that for non-commensurate fillings a onedimensional system falls into a different universality class than for the commensurate case, and will be superfluid for all non-zero values of κ. This property is also evident in the behavior of the Luttinger parameter K b (shown in figure 7(b)), as we will discuss in the next section.

Luttinger liquid parameters
The unifying concept for one-dimensional interacting systems, of both bosonic and fermionic type, is the Tomonaga-Luttinger liquid [43,44]. A system of this type exhibits two important features, namely, that correlation functions follow a power-law decay [45], 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 behavior of all the correlation functions can be expressed, and the group velocity, v, of the collective modes.

Correlation function
For the system we study, K b can be extracted from the density-density correlation function [46] å = á ñ-á ñá ñ = + + ( ) ( ) N r N n n n n 1 , 1 8 x N x x r x x r 1 which asymptotically decays as [42] p p Rather than directly fitting N(r) to equation (19), however, more stable results are obtained by evaluating its Fourier transform, N(k), and estimating K b from its derivative [47] In figure 7(a) we show N(k) for a range of values of κ. For large κ, for which the system is gapless, the results show a clearly linear behavior for small k, indicating that this procedure is valid. As κ is reduced, however, N(k) is suppressed and has a quadratic behavior for small k, corresponding to the system being gapped and thus not being describable in terms of a Luttinger liquid.
For small values of κ the system is in the Mott insulator regime, and consequently K b should diverge. In figure 7(b) 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 κ 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 [42,48].
We can also note that since Luttinger liquids with K b <2 are dominated by superfluid correlations [42], the transition will be to a superfluid state, exactly as we would expect from the vanishing of the gap seen previously. From figure 7(b) we can see that K b crosses the critical value at approximately κ=0.48. This is slightly larger than the value κ=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 κ=0, and will immediately make a phase transition as soon as κ>0. At this point K b will again take a universal value, in this case being exactly twice that of the commensurate transition [42]. In figure 7(b) 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 1 b in the limit of small κ. 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.

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 ω Here Y ñ | n is the nth excited state in the N+1 particle system, Y ñ | g is the ground state of N particles and w = + - , where E n (N, N+1) is the energy of the nth excited state in the N+1 particle system. To smooth out this spectral function we convolve equation (21) with a narrow Lorentzian. The presence of sharp peaks in A(ω, k) indicate the presence of well-defined quasiparticle excitations, and plotting the maxima of A(ω, k) allows us to determine their dispersion relation ω(k).
In the upper panels of figure 8 we show the results obtained for a system well in the superfluid regime (κ=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(ω, k) reveals that this momentum dependence corresponds to  (20)). For commensurate filling (black triangles), κ, K b diverges as κ reduces, indicating a Mott insulator. When K b passes through the universal value K b =1/2, the system makes a transition to a superfluid state. By contrast, an incommensurate system of 6 bosons in an 8-site lattice (red circles) is superfluid for all non-zero values of κ. In this case, as k  0, K b approaches unity, as expected for a system in this universality class. two regions of linear dispersion, centered on k=±π/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 5 . 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 figure 8 show the corresponding spectral functions when the system is in the Mott insulator regime (κ=0.1). In contrast to the previous result, the peaks in A(ω, 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.

Scattering processes
As shown in figures 3 and 4, and discussed in section 4, it is clear that the superfluid is composed of a pair of macroscopically occupied degenerate orbitals, with momenta ±π/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±π/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 equation (4)  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 «k k 0, 0 , , with | | k small but non-zero. In the standard case, there are energetic reasons for favoring the Figure 8. Spectral density (left panels) and line shape (right panels) A(ω, k) of a single quasiparticle excitation. The red lines in the left panels indicate the dispersion relation ω(q); the velocity of excitations (the other Luttinger parameter) being given by the slope of these curves at k=±π/2. Upper plots are for the superfluid regime κ=0.53, the lower panels for the Mott insulating regime κ=0.1 for N=8. 5 By 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 ±π/2) such that the total density fluctuation operator can be written as . 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 κ (such that  kF x 4 1 , with x 1 =2.404 the first zero of  ( ) x 0 ) this translates into small values of F.
A generic pair collision of the type - has an associated matrix element which is proportional to 6 cos 1 cos sin sin sin . 24 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 k sin will tend to cancel out, while the k cos term will average to a non-zero value. According to this criterion, the favored momenta will be those which minimize cos , which requires =q p cos cos . 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=±π/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 against the background of a macroscopically occupied π/2-momentum state, and similarly for the condensate at momentum −π/2. In both cases, the main processes are of the type (25) with small but non-zero 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 ±π/2 with depletion clouds packed near those two momenta. To the extent that the two depletion clouds can be viewed as peaked around momenta ±π/ 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 figure 4). This is also consistent with our results for the two-particle reduced density matrix, shown in figure 5.
The thermodynamic limit of Hamiltonian (A12) 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 rhs of (25) with | | k comparable to π/ 2, the result being a peculiar one-dimensional system.

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 Ω, then the tunneling term is renormalized to an effective value. For sinusoidal driving, for example, , where K is the amplitude of the shaking and  0 is the zeroth Bessel function. We show this behavior in figure 9, and we can see that when the Bessel function is zero (at K/ Ω;2.404) the effective tunneling vanishes.
We can now consider altering the parameters of the shaking [49] at a timescale much slower than the period of the shaking τ=2π/Ω. The effective hopping will now vary according to K(t)/Ω. In particular we can slowly oscillate K about the root of the Bessel function at a frequency ω, where ω=Ω, which will give a timedependent effective tunneling ( ) J t eff . For a suitable choice of the variation K(t), this will thus produce the timedependent Hamiltonian (1). In order for equation (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=ω=Ω. 6 Umklapp processes with final momenta These conditions can be met in systems of periodically-driven cold atoms such as rubidium-87 [13,20] and potassium-39 [50]. 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=Ω, and gives a wide window of available frequencies for the frequency ω 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 [51,52] 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 BH model is known to exhibit this behavior [53], 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 figure 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 dc , between lattice sites, so long as J dc /U is lower than the critical value for the standard Mott transition. Consequently we believe that its experimental realization is completely feasible.

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 ±π/2. While the creation of similar condensates has been been theoretically predicted [54,55] and experimentally observed [56] 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. Figure 9. The most frequently considered case of potential driving is to 'shake' the lattice potential sinusoidally with time. This produces an effective tunneling where Ω is the frequency of the shaking and K is its amplitude. If K/Ω is periodically varied about the first zero of the Bessel function, J eff will oscillate with same periodic time dependence in the shaded region, yielding the Hamiltonian (1).
Time averaging then provides us with the effective Hamiltonian  The dipole approximation works well for κ<0.25, as can be seen in figure B1. Even though the ansatz is not too far from the exact energy even for κ=0.4, it is unable to describe the phase transition since its density-density correlation 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.
and the average á ñ a · is taken with respect to a state that only contains species α. 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 [42] 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, where one assumes that the joint system of A and B is in a state such that equation (C9) is satisfied. Using equation (C5) one can readily conclude that K b =Q.  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 equations (C10) and (C11). Numerical evidence discussed in section 4.2 strongly suggests that our ground state is of the cat type (C10), where A, B are operators that create many bosons peaked around±π/2, respectively. We now need to justify whether this constitutes two independent boson species in the sense of equations (C1) and (C3). We focus on the homogeneous case for which the average density is position independent. In general, if we define . This permits us to infer equation (C1), which is crucial for the derivation of equation (C8). Thus, the cat-type state in equation (C10) with † A and † B creating bosons near ±π/ 2 can rightly be viewed as formed by two independent condensates in the sense of equations (C1) and (C3). 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.