Floquet-Induced Superfluidity with Periodically Modulated Interactions of Two-Species Hardcore Bosons in a One-dimensional Optical Lattice

We consider two species of hard-core bosons with density dependent hopping in a one-dimensional optical lattice, for which we propose experimental realizations using time-periodic driving. The quantum phase diagram for half-integer filling is determined by combining different advanced numerical simulations with analytic calculations. We find that a reduction of the density-dependent hopping induces a Mott-insulator to superfluid transition. For negative hopping a previously unknown state is found, where one species induces a gauge phase of the other species, which leads to a superfluid phase of gauge-paired particles. The corresponding experimental signatures are discussed.

Recently, time-dependent and driven optical lattices have opened an era of exploring exotic dynamical quantum states . For instance, assisted Raman tunneling and shaking were proposed to induce a densitydependent complex phase in the hopping elements, which may allow the experimental realization of anyonic physics [26][27][28][29]. On the other hand, a fast time-periodic modulation of the interaction [30,31] will lead to an effective hopping matrix element depending on the density difference [34][35][36][37], which gives rise to pair superfluidity in one dimension (1D) [34], while superfluidity is suppressed in higher dimensions [35]. Experimental realizations of time-periodic driving [37][38][39][40][41][42][43][44][45][46] demonstrate that signatures of interesting effective models can be observed before heating or decoherence destroys the so-called Floquet states.
In this Letter we propose a realization of a densitydependent hopping model of two interacting boson species in 1D via time-periodic driving. The corresponding quantum phase diagram is determined using a combination of advanced numerical methods. We find that a reduction of the density-dependent hopping by driving, counter-intuitively, causes a MI to SF quantum phase transition. For larger driving we obtain negative effec-tive hopping, which gives rise to an exotic SF phase of gauge-dressed composite particles.
The model for hard-core bosons with two hyperfine states (marked by "a" and "b") in a 1D optical lattice is given in terms of corresponding creation and annihilation operatorsâ † l ,â l ,b † l ,b l at each sitê where we have included a possible Rabi coupling J Ω .
Here U is the repulsive interaction between densities of opposite speciesn a l =â † lâl andn b l =b † lbl , which can be achieved by a magnetic field just below the interspecies Feshbach resonance in a deep lattice potential (see Supplemental Materials for details [47]). We consider a time-periodic modulation of the Feshbach resonance U =Ū + δU cos(ωt) without Rabi transitions J Ω = 0, oralternatively -an oscillating Rabi amplitude J Ω = J 0 Ω cos(ωt) for constant U . The concrete details for the experimental realization of the two alternative setups are described in the Supplemental Materials [47], which require a lattice depth of s ∼ 20 to obey the hardcore constraint, corresponding to J/h ∼ 20Hz for Rubidium-87 [47]. Choosing driving frequency and amplitude around 1 kHz, we have ω ≫ J,Ū , while transitions to higher bands are still suppressed. Using Floquet theory in this limit [32][33][34][35][36][37], in both cases a time-independent effective Hamiltonian with density-dependent hopping can be reached by adiabatically increasing the driving amplitudê the local densities of the opposite specieŝ with matrix elements J for n a l − n a l+1 = 0 JJ 0 [K] for |n a l − n a l+1 | = 1 (4) and analogously forĴ a l . Here J 0 [K] denotes the zerothorder Bessel function of the first kind and the dimensionless driving amplitude K = δU/ ω oralternatively -K = 2J 0 Ω / ω gives the modulation strength in units of ω [47]. As illustrated in Fig. 1(a) the effect of driving is therefore to suppress hopping of hard-core type-a bosons by J 0 [K] if the occupation of type-b bosons is different and vice versa. The suppression decreases with increasing driving K from J 0 [0] = 1 to the negative minimum value of J 0 [3.8717] ≈ −0.4024. Further tuning parameters of the model are possible, e.g. by an asymmetry in the pulse sequence [39], which makes this setup an interesting general platform. In this Letter we will focus on the phase diagram of the model (2) at half-filling n a l = n b l = 1 /2. In this case, the undriven system J 0 [0] = 1 is known to be in the Mott state for anyŪ > 0 without a quantum phase transition [48]. However, as we will see below, the selective reduction of hopping elements by driving will destroy the MI state.
An interesting point is reached at the zeros of the Bessel function since for J 0 [K] = 0 the hopping between neighboring double occupied and empty sites is not possible in this case as shown in Fig. 1 (a). Because the Hamiltonian no longer distinguishes between double occupied and empty sites, we can denote both of them with pseudo-spin up |↑ (for n a l = n b l ). Likewise hopping between neighboring single occupied sites is forbidden regardless if they are type-a or b, so both can be denoted with pseudo-spin down | ↓ (for n a l = n b l ). The corresponding total occupation numbers for the four different possible local states (a, b, double, empty) are all conserved and the resulting Hamiltonian for half-filling is expressed exactly aŝ whereŜ +/− l andŜ z l represent the respective pseudo-spin-1/2 operators. There is a macroscopic degeneracy 2 L increasing with the number of sites L, since each pseudospin state represent two different but equivalent local states for each site. The xy−model in Eq. (5) is exactly solvable, whereŪ provides a Zeemann splitting between |↑ and |↓ . ForŪ > 4J, the system is saturated with only single occupied sites and a finite charge gap corresponding to the MI phase. WhenŪ ≤ 4J, the ground state is in a gapless xy phase without SF response indicated by a blue vertical line in Fig. 1(b). Details of the solution and correlations at the degenerate line are discussed in the Supplemental Materials [47].
To obtain the full quantum phase diagram at halffilling, we now use a combination of three independent advanced numerical simulation methods. The density matrix renormalization group (DMRG) method [49][50][51][52] is used to measure properties of finite-size chains, such as the charge gap ∆ c , the superfluid density ρ s , and correlation functions using up to M = 4096 states. With the further development of the DMRG to infinite systems (iDMRG) [53][54][55], we can moreover determine the fidelity susceptibility χ F and the entanglement entropy S directly in the thermodynamic limit. Last but not least the stochastic series expansion algorithm of the quantum Monte Carlo (QMC) method with parallel tempering [56][57][58] is used to calculate the compressibility κ close to the zero temperature limit.
As shown in Fig. 2(a) for J = 0.4Ū we now observe signatures of a quantum phase transition at half-filling as a function of the effective hopping J 0 [K], which is reduced by the driving amplitude K. Because the phase transition is of Berezinskii-Kosterlitz-Thouless (BKT) type [59], finite-size effects are only logarithmically small. Therefore, measuring the transition point numerically by physical observables is very tricky and inaccurate, so we employ a combination of methods. Only in the full thermodynamic limit, the charge gap increases from zero, the global compressibility goes to zero, the entanglement entropy drops from infinity to a finite value and the fidelity susceptibility becomes extremely sharp at the transition point. The superfluid density ρ s can be obtained using DMRG from the second-order response [E 0 (θ) − E 0 (0)]/θ 2 of the ground-state energy E 0 to a twist-angle θ [60]. The response ρ s is finite and increasing for small J 0 [K], which shows that the system is indeed in a superfluid phase for this part of the phase diagram. The increase of ρ s with effective hopping J 0 [K] Fig. 2(a) is not surprising, since for smaller J 0 [K] the hopping of type-a bosons is blocked by a changing occupation of type-b and vice versa. However, for larger J 0 [K] a maximum and sudden drop to ρ s → 0 as J 0 [K] → 1 signals a quantum phase transition to the well-established Mott state in the undriven system [48,61]. To pinpoint the transition point, we consider the fidelity susceptibility χ F (x) = −2 ln F (x 1 , x 2 )/δ 2 , which is defined via the overlap of ground states F (x 1 , x 2 ) = ψ 0 (x 1 )|ψ 0 (x 2 ) with δ = |x 1 − x 2 | andx = (x 1 + x 2 )/2 for two close values x 1 and x 2 of the parameter J 0 [K]. A peak in χ F is a clear signal of a quantum phase transition [62,63], which occurs at J 0 [K] c = 0.624 (6). In addition, the entanglement entropy S = −Trρ r ln ρ r is obtained from the partial trace of the reduced density matrix for half the system [64][65][66], which shows a distinct drop in the vicinity of the transition point. Using QMC we find the compressibility κ = N 2 − N 2 for L = 100 sites at low temperatures which vanishes in the deep Mott phase. The charge gap ∆ c = E p + E h − 2E 0 is found by DMRG from the energies of systems with one additional particle E p and one additional hole E h relative to the ground state and becomes finite in the MI. After finite-size scaling analysis on J 0 [K] c by level-spectroscopic technique [47,67], we find it matches well with the maxima in χ F within errorbars, so we use the latter to obtain the full phase diagram in Fig. 1 At first sight it is strange that the reduction in hopping J 0 [K] can induce a SF state, since normally weaker hopping makes the MI more stable. However, in this case the density-dependent processes in Fig. 1 are responsible for a virtual exchange, which reduces the energy of an alternating density order ababab... to second order 4J 2 J 2 0 [K]/Ū [14]. Therefore, by selectively tuning away those processes via periodic driving, the alternating order and the corresponding MI are actually destabilized, which in turn enables a SF for finiteŪ . For J 0 [K] = 0 the system has no a-b-density correlations, which leads to the degeneracy discussed above.
It is instructive to analyze the characteristic correlation functions for the different phases as shown in Fig. 2(b) for J = 0.4Ū . The single-particle correlation G c (r) = â † 0â r shows a typical power-law decay in the r on the other hand shows a slow power-law decay in either phase.
We now turn to negative effective hopping J 0 [K] < 0. The corresponding phase diagram and the superfluid density are shown in Fig. 1(b) and Fig. 3, respectively. At first sight the results look perfectly symmetric around J 0 [K] = 0, which would suggest that negative hopping has the same effect as positive hopping. However, the underlying states for positive and negative values are quite different, which becomes clear by looking at the signature as a function of momentum k, where w(k) stands for the Fourier transformation of the Wannier function in a 1D optical lattice with lattice spacing equal to one [68]. As shown in the inset of Fig. 3 the MD shows an interference pattern with sharp peaks at k = 0 (modulo 2π) for positive values J 0 [K] = 0.35, which originates from the phase coherence of bosons in the normal SF. However in the region J 0 [K] < 0, no sharp interference pattern is observed. Both the symmetry in ρ s and the difference in the MD interference pattern can be explained by a gauge transformation which defines new quasi-particles of type-β β l =b l exp(iπn a l ) and analogous for type-α. We see that the hopping terms in Eqs.  is therefore equivalent to a transformationb l →β l and a l →α l in Eq. (7). Accordingly, the energies and phase transition lines are identical for positive and negative J 0 [K], but the superfluid density for negative sign corresponds to a response of gauge-paired particlesα,β and is therefore called a gauge-dressed SF with a different MD shown in the inset of Fig. 3. The transition to such an exotic condensed density can also be captured by a Gutzwiller mean-field argument which is discussed in the Supplemental Materials [47]. Note, that the symmetry transformation to new gauge-paired particles in Eq. (7) is independent of the dimensionality and geometry of the lattice.
Thus, the gauge dressed SF is characterized by a lattice gauge exp(iπn a l ) provided by one species (type-a) which couples to the hopping of the other species (typeb) and vice versa. As can be seen from Eq. (7) the gauge dressed hopping becomes dominant in the strongly driven region J 0 [K] < 0, resulting in a superfluid response from gauge dressed particles. The quantum phase transition to a MI is analogous to an ordinary SF and happens at exactly the same critical value of J/Ū in Fig. 1(b) as for corresponding positive J 0 [K] > 0 since the gauge does not change the energy response to a twist-angle θ. The gauge dressed SF is therefore different from pair superfluidity, where correlated hopping is observed due to a strong coupling of the hopping directly to the density [34,69]. The so-called counterflow SF is another type of correlated hopping [13][14][15], where hopping of particles of one species is facilitated by holes of the opposite species. In contrast, in the new gauge dressed SF the hopping is facilitated by gauges exp(iπn a l ), which can also be viewed as particles that are their own anti-particles, analogous to a Majorana description.
For the experimental realization of these phases several critical questions must be solved. First of all accessing the steady state by adiabatic ramping of the driving amplitude from the ground state is only possible, when no dense avoided level crossing of the Floquet quasi-energy take place. Our analysis of the quasi-energy spectrum in the Supplemental Materials [47] ensures that there are no critical avoided level crossings in the relevant parameter range. Secondly, a measurement can be affected by the unitary transformation into the effective Floquet basis, if the operators do not commute with the Kick operator [47]. For stroboscopic measurements at times of integer multiples of period T = 2π/ω we show in the Supplemental Materials [47] that this effect is reduced by J/ ω in the high frequency limit and calculate the corrections from higher order terms, in order to predict the experimental mapping of the phase diagram by time-of-flight and compressibility measurements. For a separate check of the predictions we also performed real-time simulations for a small lattice L = 6 [47] which clearly show the stability of the effective Hamiltonian and the feasibility of real-time dynamic measurements on finite time-and length-scales.
In conclusion, we proposed a setup of a 1D lattice with two species of hard-core bosons and time-periodically modulated fields, which can be described by densitydependent tunneling with an interesting quantum phase diagram. By controlling the driving amplitude, density dependent hopping processes are selectively tuned away, which are responsible for an alternating density a-b order. This in turn leads to a transition from the MI to a SF at half-filling in contrast to the undriven case. By tuning away these terms completely at J 0 [K] = 0, a highly degenerate state is obtained corresponding to an exactly solvable model without a-b correlations. For many-body systems the study of nearly degenerate points is a very active research area, e.g. in the context of frustrated models, spin ice, and spin liquids. Much theoretical activity is devoted to studying novel quantum states, which are dominated by the quantum fluctuations near degenerate points, but we are not aware of any such studies for driving-induced degeneracy. In this case, dynamical effects will likely dominate the quantum correlations, which opens an interesting research field beyond our current abilities. For even larger driving amplitudes, negative hopping parameters J 0 [K] < 0 lead to a new gauge dressed SF with a novel type of pairing mechanism, where an atom of one species and a gauge phase of the other are bound to contribute to a nonzero superfluidity. This gauge dressed SF has different correlations from an ordinary SF, as shown in Fig. 3 for the momentum distribution. Nonetheless, an exact hidden transformation to the positive hopping case can be found.
We thank Youjin Deng, Shaon Sahoo, Oliver Thomas, and Zhensheng Yuan for useful discussion. This research was supported by the Special Foundation from NSFC for theoretical physics Research Program of China

APPENDIX
Here, we discuss two proposals of the experimental realization, higher orders in the effective Hamiltonian, the role of the Kick operator, avoided level crossing of the quasi-energy spectrum, real-time dynamics of the original time-dependent Hamiltonian, finite size scaling, and the Gutzwiller mean field method.

EXPERIMENTAL REALIZATION
Initially we assume to have ultra-cold atoms in one hyperfine state confined in an optical dipole trap. [23] The trap consists of a pair of counter-propagating laser beams with wave length λ L along the x-axis, so the potential function reads with the lattice depth V 0 and the wave vector k r = 2π/λ L . Using an initially off-resonant radio-frequency magnetic field, we suggest to adiabatically ramp the Zeeman fields to zero and decrease the radio-frequency coupling strength to a certain value for a while. Then we propose to suddenly turn off the coupling strength projecting the BEC into an equal superposition of two hyperfine states. In the single-band approximation, the movement of the atoms appears in form of hoppings between two neighboring sites, e.g. the minima of the lattice potential with spacing λ L /2, and the related Hamiltonian readsĤ whereâ l (b l ) andâ † l (b † l ) are the annihilation and creation operators, J a (J b ) the hopping coefficient of atoms with hyperfine level "a" ("b") and the site index l runs over the whole lattice. Obviously, J a = J b = J because the hopping processes are independent of the hyperfine internal states of atoms. [9] The depth of the optical lattice potential can affect the on-site repulsive interaction between the ultra-cold atoms, [9] which is independent of the hyperfine states unless we are close to a Feshbach resonance. Therefore, increasing the lattice depth V 0 will generate large intra-and inter-species repulsive interactions independent of the hyperfine stateŝ wheren a l =â † lâ l (n b l =b † lb l ) denote the particle number operator of the species "a" ("b") and U L ≫ J. Furthermore, we suggest to add a static magnetic field and tune its amplitude B to be very close to the Feshbach resonance point B 0 , where two-species atoms form s-wave bound states, while it is far away from intra-species Feshbach resonance points for both species. [7,70] On the side of the negative scattering length, an attractive inter-species interaction emerges, namelyĤ with U F ≫ J. It can compensate the inter-species repulsion U L and results in a total finite inter-species repulsion U = U L − U F , which is assumed to be of the order of J. Meanwhile a large intra-species repulsion U L still leads to a hard-core constraint, which means more than one atom from the same species is forbidden at the same lattice . Similarly an atom jumps from "m" to "b" by emitting a photon from a Stokes laser (blue line, frequency ωS, linearly polarized, coupling strength ΩS). [72] (c) Schematic picture of a two-level system directly coupled via a pump laser. An atom jumps from the hyperfine state "a" to "b" by absorbing a photon from a pump laser (red line, frequency ωP , fast rotating linearly polarized, coupling strength ΩP ). One realization of the pump laser is the output of a circularly polarized laser passing a 1/4 wave plate which is connected to a fast rotating mechanical motor (frequency ω sets several kHz).
site. Taking 87 Rb atoms for example, we can choose the magnetic field a little smaller than the inter-species Feshbach resonance point 1259.96 G and far away from the intra-species ones, which amounts to be 685.43 G for |F = 1, m F = 1 ("a") and 661.43 G for |F = 1, m F = 0 ("b"). [70] To fulfill the hardcore constraint, the lattice depth V 0 must be choosen significantly larger than the recoil energy For Rubidium-87 we therefore have E r /h ∼ 3.5kHz in a 400nm lattice, which gives J/h ∼ 17Hz. The rotating frequency of the time-periodic driving must be much larger than J and U , but not too large to avoid "photon-assisted hopping" between different energy bands of the optical lattice. In the following discussion of possible realizations we therefore assume a rotating frequency ω of the order of 1 kHz, which will also be the order of magnitude of the driving amplitude.

Periodically modulated inter-species interaction
A straight-forward, but technologically challenging time-periodic driving can be realized by an oscillating magnetic field B(t) =B + δB cos ωt near B 0 , whereB denotes the time-average strength of the magnetic field, δB represents the oscillation amplitude of the magnetic field and ω stands for the oscillating frequency as shown in Fig. 4 (a). Thus the relevant s-wave scattering length can be written as where ∆ is the width of the Feshbach resonance and a bg represents the background scattering length, which is determined by the lattice depth. If we choose δB ≪ |B − B 0 |, we can further perform a Taylor series expansion of a s (t) with respect to a small value of δB/(B − B 0 ) and get where the coefficients of the leading orders are a Note that a pure cosine oscillation of a s is in principle also possible for larger amplitude δB if the waveform of the magnetic field is adjusted correspondingly. The inter-species interaction energy U (t) is proportional to the related scattering length, which means U (t) =Ū + δU cos(ωt), where the time-average energyŪ is proportional to a (0) s and the oscillation amplitude δU is proportional to a (1) s if we neglect higher-order terms. In this case, the system can be described by the following Hamiltonian:Ĥ Periodically modulated Rabi oscillation Fast oscillating fields are possible but challenging, so an alternative experimental realization in a static magnetic field is useful. To this end we propose to gradually switch on a pair of Raman laser beams, which are coupled to the atomic cloud with the frequency difference δω R . An atomic transition from the hyperfine state "a" to "b" by a two-photon emission-absorption has the standard Λ-form [72] shown in Fig. 4 (b). One pump laser initiates that atoms jump from the hyperfine state "a" to the intermediate state "m" by absorbing a photon with frequency ω P and coupling strength Ω P , while the other Stokes laser triggers atoms to jump from "m" to "b" by emitting a photon with frequency ω S and coupling strength Ω S . [72] Both coupling strengths depend on the projection of a transition dipole moment onto the polarization of the Raman laser beams, [72] namely where d P/S denote the transition dipole momenta and E P/S represent the electric field of the respective Raman laser beams. The transition dipole momenta only depend on the initial and final states and are unchanged during the two-photon transition processes. At the stimulated Raman transition, the frequency difference δω R exactly coincides with the hyperfine splitting between two levels ∆ ab . In order to avoid a resonant excitation of the intermediate state, the detuning ∆ < 0 of the Raman beam from the one-photon transition has to be much larger than the linewidth of the excited level and it should also be much larger than the Raman detuning δ from the two-photon resonance. In this way, the system can be considered as an effective two-level systems "a" and "b". Its Hamiltonian readŝ in the basis of the two hyperfine levels with the energies E P/S = Ω 2 P/S /∆ and the so-called "Rabi frequency" J Ω = Ω P Ω S /∆. After the second quantization, we obtain the Hamiltonian for the Rabi oscillation where we neglect the small shift of the chemical potential δ between two hyperfine levels.
In order to produce a time-periodic oscillating Rabi coupling strength, we let the polarization direction of the pump laser circulate in time, e.g. E x = A cos(ωt) and E z = A sin(ωt) with the amplitude of the polarization A. For realization, a circularly-polarized pump laser may pass through 1/4 wave plate to get a linearly-polarized laser beam as output, the polarization direction of which is 45 degree shifted to the optical axis of the wave plate. Sequentially, the wave plate is connected to a mechanical motor with rotating frequency ω in the kHz range, which is much lower than the laser frequency of several hundreds of THz (10 12 Hz). Therefore the polarization of the pump laser is also rotating and effectively provides a time-modulated Rabi coupling. In order to avoid coupling to other magnetic sublevels when the linear polarization is rotated, we assume a sufficiently strong Zeeman splitting. Therefore assuming d z = 0, we get with Ω 0 P = Ad x . As a result, the effective Rabi-frequency turns out to be time-periodic J Ω (t) = J 0 Ω cos(ωt) (19) with the amplitude J 0 Ω = Ω 0 P Ω S /∆. We can use, for instance, an acousto-optic modulator (AOM) to change both the amplitude and the polarization of the pump laser. [73,74] In an extra scheme, we suggest to prepare a cloud of ultra-cold atoms with equally-weighted species but with different total angular momenta and long life time. They directly couple to a fast rotating linearly polarized pump laser beam with rotating frequency in the kHz range. In this way, we can generate the same time-periodic Rabi oscillation in the Eq. (9), see Fig. 4 (c).
In general, the full Hamiltonian reads then

EFFECTIVE HAMILTONIAN AND KICK OPERATORS
For the wave function of a system, which is described by a time-periodically driven HamiltonianĤ(t) =Ĥ(t + T ) with the period T = 2π/ω being determined by the driving frequency ω, the time-dependent Schrödinger equation holds where the steady-state after long times obeys the Floquet theory. [75] As a result, the wave function is of the form |ψ(t) = e −iǫt/ |φ(t) with time-periodic Floquet modes |φ(t) = |φ(t + T ) , where ǫ represents the corresponding quasi-energy. With this the time-dependent Schrödinger equation (21) becomeŝ with the Floquet HamiltonianĤ(t) =Ĥ(t) − i ∂/∂t. In the extended Hilbert space including the time dimension, Eq. (22) is the eigenvalue equation of the Floquet HamiltonianĤ(t). In the high frequency regime, where ω is large, it is natural to choose the eigenfunctions of the operator −i ∂/∂t as the basis and to treat the whole Hamiltonian H(t) perturbatively. Then, we can use a nearly-degenerate perturbative method to solve Eq. (22). [76] Although the full solution of Eq. (22) is still difficult because of the intricate quasi-energy spectrum, we can calculate the timeindependent effective Hamiltonian and the kick operator order by order. Thus, the whole dynamics of the system can be described by the effective Hamiltonian and the corresponding kick operator. However in our cases we can not perform directly perturbative calculations because of the extra energy scales δU and J 0 Ω ∼ ω. Therefore, we need to apply a rotationV at the preliminary step in order to eliminate these extra terms by moving them to the phases of the respective hopping terms. [77] In the rotating frame after the unitary transformation, the wave function turns out to obey the transformed time-dependent Schrödinger equation the quickly rotating phases have been absorbed in the transformation.

Periodically modulated inter-species interaction
For the periodically modulated inter-species Hamiltonian (14) the rotation operator is given byV (t) = exp[−iK U ln a ln b l ] with dimensionless modulation strengths K U = δU/ ω andK U = K U sin(ωt). With this the Hamiltonian after the rotation results in HereĤ r is again periodic and can thus be expanded into a Fourier seriesĤ r = where J n denotes the nth order Bessel function of first kind. Now we can calculate the effective Hamiltonian order by order. [76,77] The zeroth-order effective Hamiltonian iŝ Furthermore, the first-order effective Hamiltonian vanisheŝ where we use the propertyĤ (n) r = (−1) nĤ (−n) r . All second-order corrections consist of many terms, which are accompanied by the prefactor (J/ ω) 2 and are not listed here. Because of J/ ω ≪ 1 we conclude that all higherorder corrections to the zeroth-order effective Hamiltonian are small.
In order to describe the whole dynamics of the system, we also calculated the kick operator up to first order: The effective Hamiltonian and the kick operator calculated above are defined in the rotating frame, but we are interested in observables in the lab frame. The link between both frames is provided by the fact that the time evolving operatorsÛ (t 2 , t 1 ) in the laboratory frame andÛ r (t 2 , t 1 ) in the rotating frame are connected by a rotation transformation, namelyÛ (t 2 , t 1 ) =V (t 2 )Û r (t 2 , t 1 )V † (t 1 ). In this paper, we are only interested in the stroboscopic dynamics at time t = nT , so we concludeÛ (t 2 , t 1 ) =Û r (t 2 , t 1 ) and any observable turns out to be the same in both frames. Furthermore, if one prepares the system in the ground state of the non-driven Hamiltonian, then adiabatically turning on the driving has the consequence that the ground state of the system will follow the instantaneous stroboscopic Floquet Hamiltonian. [77] Thus the time-evolving wave-function consists of the ground state of the effective Hamiltonian and a phase factor from the kick operator Up to the first order, we getK(nT ) =K (1) (nT ) = 2Ĥ (1) r /i ω. The expectation value of an observableÔ results then in where Â e ≡ ψ e |Â|ψ e . Thus, the expectation value of an observable in the lab coincides with that of the dressed observable e

2Ĥ
(1) r ωÔ e − 2Ĥ (1) r ω , which is determined by the effective Hamiltonian. As J/ ω is small, we only need to keep the two lowest orders As a concrete example we take the expectation value ofâ kâq , which has been used for calculating the density distribution in momentum space â kâq = â kâq e + 2J 1 [K] Â e J/ ω witĥ We read off that the correction operatorÂ includes finite local terms. Note that one can always reduce the effect of the correction by tuning the value of J/ ω.

Periodically modulated Rabi oscillation
We now deal with the Hamiltonian (20), which describes a periodically modulated Rabi oscillation. In case of δU = 0 and U =Ū , the rotation transformation is given byV = exp[−iK Ω (â † lb l +h.c.)] with dimensionless modulation strengths K Ω = J 0 Ω / ω andK Ω = K Ω sin(ωt), so we get From this we concludê Thus, the rotated Hamiltonian results in where we haveĤ and the interaction term remains unchanged, e.g.Ĥ U r (t) =Ū ln a ln b l . Also hereĤ r (t) is time-periodic and can be expanded into a Fourier series, namelyĤ for even and odd orders, respectively. Now we can use the formalism of the high-frequency expansion to calculate order by order the effective Hamiltonian [76,77] in the rotating frame, namelyĤ e = +∞ n=0Ĥ (n) e . The zeroth order effective Hamiltonian turns out to beĤ (39) and the first order vanishesĤ We also calculate the kick operator [76,77]K e = +∞ n=0K (n) e up to first order: With the same reasoning as described in the last subsection we only check the dynamics of the system stroboscopically at t = nT , so the expectation value of an observableÔ is given by The expectation value of an observable in the lab frame coincides with exp(2Ĥ (1) r / ω), which is the one of the dressed observable being determined by effective Hamiltonian. As J/ ω is small, we only need to keep the two lowest orders As an example we also consider the expectation value ofâ kâq , which is used for calculating the density distribution in momentum space â kâq = â kâq e + 2J 1 [K] Â e J/ ω witĥ We read off that the correction operatorÂ includes finite local terms. Again one can always reduce the effect of the correction by tuning the value of J/ ω to lower values. In general we conclude for both experimental realizations that the zeroth-order effective Hamiltonian readŝ with K = K U = 2K Ω . As the first-order effective Hamiltonian vanishes, only this zeroth-order effective Hamiltonian is relevant. Although the value of any observable in the lab differs from the ground state expectation value of the effective Hamiltonian, the difference between them is always accompanied by a prefactor J/ ω, which can be decreased by increasing the driven frequency ω.

REAL-TIME DYNAMICS
In the following we check the real-time dynamics for a small system and answer the question if it is feasible to complete the preparation of the sample and the measurement before the thermalization sets in. To this end we consider two steps for switching on K, i.e. K U or K Ω , andŪ , respectively. At the first step, we initialize the system staying at the ground state of the non-driven model with a small on-site repulsionŪ i . Then for t > 0 the amplitude of the time-periodic circularly-polarized Raman laser beams is gradually switched on following a linear function of time t where the modulation lasts for the duration of t 1 = n 01 T until K reaches a desired value K f , so the effective speed amounts to v K = K f /t 1 . It persists as a working Rabi oscillation before the measurement. In the second step, we turn on gradually the on-site interaction as a linear function of time t where the modulation lasts for the duration t 2 − t 1 = n 12 T until the on-site interaction reaches a desired valueŪ f , so the effective speed reads vŪ = (Ū f −Ū i )/(t 2 − t 1 ). After the modulation duration we expect that the low-energy behavior of the system can be described by the effective Hamiltonian with the desired parameter valuesŪ f and K f . And then the state should persist for a while in order to complete the measurement. Here we exploit the fourth-order Runge-Kutta method to solve the time-dependent Schrödinger equation (21) and to determine the time-evolving wave function |ψ(t) . In order to understand how close it is to the desired one |ψ e , we measure the time-dependent overlap P = | ψ(t)|ψ e |. Furthermore, we measure the time-dependent energy e 0 (t) = ψ(t)|Ĥ T +Ĥ U |ψ(t) and the structure factors  (21). Here we choose J = 1 as an energy unit, a relatively high frequency ω = 20, KΩ = 0, and balanced filling N a = N b = 3. Initially, we let the system stay at the ground state of the Hamiltonian with a small on-site repulsion strengthŪ i = 1. Then KU linearly grows until t1 and persists as a working Rabi oscillation. At the moment t1, we start to linearly increaseŪ to a desired valuē U f at t2, which persists until the measurement. We consider two cases: t1 = 100 T and t2 = 200 T for the fast switching-on (upper four clusters of panels) while t1 = 300 T and t2 = 600 T for the slow switching-on (lower four clusters of panels). For each case we show the time-evolving behavior related to four different parameter sets: (1) and (5) for K f U = 1 andŪ f = 2; (2) and (6) for K f U = 4 andŪ f = 2; (3) and (7) for K f U = 1 andŪ f = 6; (4) and (8) for K f U = 4 andŪ f = 6. For each scheme we plot the modulated parameters KU (black lines) andŪ (red lines) as a function of time t/T in panel (a). And we plot the time-evolving overlap P = | ψ(t)|ψe | and the energy e0(t) = ψ(t)|ĤT +ĤU |ψ(t) in panel (b). In panel (c) and (d), we exhibit the structure factors ofñ b k andñ β k at the moments t = 0 (blue circles), t = t1 (magenta squares), and t = t2 (green diamonds), respectively.
whereβ l =b l exp(iπn a l ) andβ † l =b † l exp(−iπn a l ) are the annihilation and creation operators for the gauge-dressed particles. In the Figs. 5 and 6, we systematically study the real-time dynamics for 6 sites for a fast and a slow switching-on. Although the middle process is complicated, the time-evolving overlap is close to unity and the energy e 0 is almost constant at the end of the modulations. This means that we can obtain the ground state of the effective Hamiltonian with the desired physical parameters following our scheme of the sample preparation. Besides we note that a slow switching-on always works better than the fast one, so thus we suggest that experimentalists need to tune the parameters as slow as possible before the thermalization happens.

AVOIDED LEVEL-CROSSINGS
After a cloud of ultracold atoms is confined in the optical lattice and has reached equilibrium in the ground-state, the relevant parameters K U , K Ω , andŪ are adiabatically switched on in order to obtain the ground state of the effective Hamiltonian in the specified parameter regime. In previous studies it was found that the request to the adiabatic modulation is not achievable if avoided level-crossings between different Floquet bands occur. [78] In this section, we investigate avoided level-crossings in the quasi-energy spectrum for a small system as a function of K and U , respectively. Besides the perturbative treatment discussed above, for a small system size we can exactly diagonalize the general Floquet HamiltonianĤ(t) =Ĥ(t) − i ∂/∂t which obeys the eigenvalue equation where the Floquet mode |φ α (t) is a many-body state instead of the local Fock basis. The Floquet modes live in the Hilbert space of real dimensions D P . Because we have |φ α (t) = |φ α (t + T ) , each Floquet mode can be expanded by  positions of the avoided level-crossings, and thereby the valid parameter regime, which we can reach through the adiabatic switching, by investigating the positions of the avoided level-crossings. In Fig. 7, we depict the quasi-energy spectrum of 6-sites as a function of K U , K Ω , and U , respectively. We choose the integer-1 filling N a + N b = 6 and a relatively high frequency ω = 20. Specially we have chosen N a = N b = 3 when K Ω = 0. WhenŪ /J = 1 is fixed in panel (a) and (c), no problem of generating an adiabatic modulation of K U or K Ω occurs. Furthermore, no extremely dense avoided level-crossings are found. When K U = 1 is fixed in panel (b) and K Ω = 1 in panel (d), we find a bunch of dense avoided level-crossings occurring in the vicinity ofŪ /J ≈ 18. In comparison with the ground-state phase diagram, where the main interesting phases happen whenŪ /J < 4, the avoided level-crossings do not appear in this region. Thus we conclude that all the phases can be achieved by an adiabatic modulation ofŪ .

FINITE-SIZE SCALING
Because of logarithmic corrections, it is challenging to derive the accurate position of the BKT-type transition point from the Mott-insulator to the superfluid phase by the finite-size scaling of the charge gap at zero temperature or the compressibility χ at low temperatures. In Fig. 8, all the curves L∆ c collapse for small J 0 [K], which means that the charge gap ∆ c scales like 1/L deep in the superfluid region. In the deep Mott-insulator region, when J 0 [K] is large, ∆ c remains finite in the thermodynamic limit. At the anticipated critical point J 0 [K] c = 0.624(6), we find that the curves L∆ c with different system sizes get slowly close to each other, but reveal no level-crossings. Similarly at the low temperature β/Ū = 2L, the compressibility converges to a finite value and to zero in the deep superfluid and the deep Mott-insulator region, respectively. However the turning points for the finite system are slowly approaching Note that under the Jordan-Wigner transformation our model can be mapped to a density-dependent hopping Fermi-Hubbard model. Thus, with the help of the operator analysis involving the level-spectroscopic technique, [67] we can choose the level-crossing of two representative excited states to be the quasi-critical point for the finite system. In the superfluid region, the representative excitation is a particle or a hole if one adds or removes an atom. Whereas in the Mott-insulator region, the lowest-excitation is a pair of a particle "a" together with a hole "b" or vice versa.
The former has a gap ∆ + c = E p − E 0 +Ū /2 measured from the energy of system, where we put one more particle E p relative to the ground-state energy E 0 . The latter has a pseudo-spin gap ∆ 0 s = E 1 − E 0 with the first-excitation energy E 1 in the same Hilbert space N a = N b = L/2 for the ground state.
In Fig. 9 (1b) and (2b), the curves of two excitation gaps reveal level-crossings for various finite system sizes. They scale very well as a linear function of 1/L in the insets and give us the position of BKT-type critical points in the thermodynamic limit. The extrapolated values are also consistent with peaks of the fidelity susceptibility, which are obtained from iDMRG calculations.

GUTZWILLER MEAN-FIELD
Here we exhibit the details of applying the Gutzwiller mean-field (GWMF) method to the problem at hand. Because of half filling N a = N b = L/2 and the hard-core constraint, the probabilities to occupy a site by a particle pair or a hole, as well as that by a single atom a or b are the same. Therefore, taking into account the normalization condition, we perform for the wave function an ansatz in terms of the uniform product matrix state ,0 sin ϕ|0, 0 l + e iφ0,1 cos ϕ|0, 1 l + e iφ1,0 cos ϕ|1, 0 l + e iφ1,1 sin ϕ|1, 1 l , where |n a , n b l denotes the local basis at site l with n a and n b standing for the numbers of species a and b, respectively. Furthermore, ϕ and φ n a ,n b represent variational parameters, which are determined below. With this the average energy per-site yields where we have introduced the abbreviation δφ = φ 0,0 − φ 1,0 − φ 0,1 + φ 1,1 . Minimizing the energy determines the wave function of the ground state. As 1 − cos 2 2ϕ is always larger than zero, the choice of the value of δφ in the ground-state wave function depends on the sign of

INTEGRABLE POINT
For the effective Hamiltonian in Eq. (2) of the main text the hopping term of the hardcore species "a" consists of two partsĤ a =Ĥ (1) a +Ĥ (2) a , wherê Here we use the hardcore constraint conditionâ lâ † l +â † lâ l = 1. The second part depends on the J 0 [K], i.e. the normalized driven amplitude K, while the first one does not. Similarly, the hopping term of the species "b" readŝ Furthermore, the onsite interacting termĤŪ =Ū  . In (1b) and (2b), we adopt the level-spectroscopic technique to achieve finite-size scaling. In the first step, we calculate the excitation gaps ∆ + c = Ep − E0 +Ū /2 (red hexagon) and ∆ 0 s = E1 − E0 (black pentagon), where E0 and E1 are the ground state and the first excited state in the Hilbert space N a = N b = L/2, respectively, and Ep is the lowest energy of system, where we add one more particle relative to the ground state. Obviously we find a level-crossing between them called quasi-critical points such as (1b) J0[K]qc = 0.415 for L = 16 and J/Ū = 0.28, (2b) J0[K]qc = 0.924 for L = 24 and J/Ū = 0.4. In the second step, we plot these quasi-critical points as a function of 1/L in the inset and find that they can be linearly extrapolated to the thermodynamic limit quite well. We get the best extrapolation values On the site l we can define number operators of a single-"a", a single-"b", a hole, and an "ab" pair, respectively, namelyN . Because of the vanishing commutator Ĥ a(b,Ū ) ,N a(b,h,p) t = 0, the total numbers of single-"a", single-"b", hole, and "ab" pair are all conserved in any eigenstate of the Hamiltonian. Furthermore, hopping terms of the species "a" and "b" can be divided into four individual and equivalent exchange processes: either "a" from "b" or vacuum from "ab" pair. Furthermore, we find that the structure of subspaces is invariant, if we replace either "a" ("b") by "b" ("a") or replace the vacuum ("ab" pair) by "ab" pair (vacuum), which leaves N a t + N b t and N p t + N h t unchanged. This means that forŪ = 0 the Hamiltonian has a larger hidden symmetry D = Z L 2 . Let us therefore play a trick of preserving the hidden symmetry D as an inner one and regrouping the four states: both "a" and "b" belong to the group "spin-down ↓", while both vacuum and "ab" pair belong to the group "spin-up ↑". With this the Hamiltonian becomesĤ Usually, the onsite interacting term with finiteŪ breaks the exchange symmetry D such that the vacuum is inequivalent to the "ab" pair. However the symmetry D can be recovered in the case of the integer-1 filling Both vacuum and "ab" pair contribute toŪ /2, while neither "a" nor "b" have any contribution. And thus theŪ -term can be considered as an effective external magnetic field, which is applied to a redefined spin-1/2 and the effective Hamiltonian in the reduced Hilbert space readŝ where have introducedŜ z l = 2 (n a l − 1/2) n b l − 1/2 . By using the common Jordan-Wigner transformation, this becomes an integrable model in the language of spinless fermions. We know that the ground state energy is −2J Supposing that the twisted angle is θ, the hopping terms in the Hamiltonian becomeŝ where the exchange process between the single "a" ("b") and the "ab" pair carries a positive phase, while the one between the single "a" ("b") and the vacuum carries a negative phase. The hidden inner symmetry D disappears and, thus, the model can not be mapped to the effective spin-1/2 XY model. Again the onsite interacting term is invariant. The ground-state energy with infinitesimal twisted angle θ ≪ 1 can be expanded in the vicinity of θ = 0, namely where the first-order term disappears because for θ = 0 the system holds the time-reversal symmetry and has no residual "current". In fact we can prove that the energy response vanishes in case of integer-1 filling. From the effective model, the ground-state occurs when N a t + N b t = N p t + N h t = L/2. Therefore we conclude N p t = N h t = L/4. We have L/2 relevant Hilbert subspaces: L/4 subspaces are connected by the exchange processes between single occupations on the site-1 and holes on the site-L carrying a twisted phase exp(iθ), while the other L/4 subspaces are connected by ones between single occupations on the site-1 and pairs on the site-L carrying a twisted phase exp(−iθ). As a result, the residual twist phase in this group vanishes under the gauge transformation. Thus, the ground state has no energy response to the twisted phase on the boundary, which means that their superfluid density vanishes.