Engineered Nearest-Neighbour Interactions with Doubly Modulated Optical Lattices

TL;DR

This paper proposes a method to engineer specific quantum interactions in optical lattices using doubly modulated fields, enabling control over quantum phases like Mott insulators and density waves.

Contribution

It introduces a novel doubly modulated Bose-Hubbard model that allows tuning of nearest-neighbour interactions and density-assisted tunnelling in optical lattices.

Findings

Doubly modulated fields induce nearest neighbour interactions.

Designing bi-chromatic driving patterns enables phase control.

Realization of Mott insulator and density wave phases.

Abstract

Optical lattice systems provide exceptional platforms for quantum simulation of many-body systems. We focus on the doubly modulated Bose-Hubbard model driven by both time-dependent on-site energy and interaction, and predict the emergence of the nearest neighbour interaction and density-assisted tunnelling. By specifically designing a bi-chromatic driving pattern for a one dimensional lattice, we demonstrate that the doubly modulated fields can be tuned to realize desired quantum phases, e.g. the Mott insulator phase with selective defects, and density wave phase.

Figures

Equations

$$\hat{H}(t)=\sum _{⟨jk⟩}{\hat{c}}_j^{†}{J}_{jk}{\hat{c}}_k+\frac{U(t)}{2}\sum _i{\hat{n}}_j({\hat{n}}_j-1)+F(t)\sum _j{ϵ}_j{\hat{n}}_j$$

$$U_I(t)=\exp \left(i\left(\sum _j{\theta }_j(t)n_j+\frac{\Gamma (t)}{2}\sum _j{n}_j(n_j-1)\right)\right),$$

$$\begin{array}{r}\tilde{H}(t)=\sum _{⟨jk⟩}{\hat{c}}_j^{†}{\hat{A}}_{jk}(t){\hat{c}}_k+\frac{U_0}{2}\sum _j{\hat{n}}_j({\hat{n}}_j-1),\end{array}$$

$${\hat{A}}_{jk}(t)=J_{jk}{e}^{i{\chi }_{jk}(t)}{e}^{i\Gamma (t)({\hat{n}}_j-{\hat{n}}_k)},$$

$$\left[c_j^{†}{c}_k,\left[{\hat{n}}_k({\hat{n}}_k-1),c_k^{†}{c}_j\right]\right]=4{\hat{n}}_j{\hat{n}}_k-2{\hat{n}}_j^2.$$

$$\begin{array}{r}{\hat{H}}_{eff}^{int}=\sum _{⟨jk⟩}\frac{U_{jk}^N}{2}{\hat{n}}_j{\hat{n}}_k+\sum _j\frac{U_j^O}{2}{\hat{n}}_j({\hat{n}}_j-1),\end{array}$$

$$\begin{array}{rl}{U}_{jk}^N& =-4\sum _{n\ne 0}\left[\frac{1}{\omega n}(g_{jk}^{-n}{t}_{kj}^n+t_{jk}^n{g}_{kj}^{-n})-\frac{U_0}{2n^2{\omega }^2}{g}_{jk}^n{g}_{kj}^{-n}\right],\\ U_j^O& =U_0-\sum _k{U}_{jk}^N,\end{array}$$

$$\begin{array}{rl}{g}_{jk}^n& =\frac{1}{T}{\int }_0^{T}dt{J}_{jk}{e}^{i{\chi }_{jk}(t)}{e}^{-in\omega t},\\ t_{jk}^n& =\frac{i}{T}{\int }_0^{T}dt{J}_{jk}{e}^{i{\chi }_{jk}(t)}\Gamma (t)e^{-in\omega t},\end{array}$$

$$U_f(t)=2U_d\cos (2\omega t)$$

$$F(t)=F_1\cos (\omega t)+2F_2\cos (2\omega t)$$

$${\hat{H}}_t=\sum _{jk}{c}_j^{†}{\hat{A}}_{jk}^0{c}_k,$$

$${\hat{A}}_{jk}^0=-j_0{J}_0({\tilde{F}}_1(j-k),{\tilde{F}}_2(j-k)+{\tilde{U}}_d({\hat{n}}_j-{\hat{n}}_k)).$$

$$c_j^{†}{\hat{A}}_{j,j+1}^0{c}_{j+1}|{\psi }_i⟩=h_{n_j,n_{j+1}}^L{c}_j^{†}{c}_{j+1}|{\psi }_i⟩,$$

$$h_{n_j,n_{j+1}}^L=-j_0{J}_0({\tilde{F}}_1,{\tilde{F}}_2+{\tilde{U}}_d(-n_j+n_{j+1}-1)),$$

$$h_{n_j,n_{j+1}}^R=-j_0{J}_0({\tilde{F}}_1,{\tilde{F}}_2+{\tilde{U}}_d(-n_j+n_{j+1}+1)).$$

$$\begin{array}{rl}{h}_{(0,1)}^L=h_{(1,0)}^R& =-j_0{J}_0({\tilde{F}}_1,{\tilde{F}}_2),\\ h_{(1,1)}^L=h_{(2,0)}^R& =-j_0{J}_0({\tilde{F}}_1,{\tilde{F}}_2-{\tilde{U}}_d),\\ h_{(1,1)}^R=h_{(0,2)}^L& =-j_0{J}_0({\tilde{F}}_1,{\tilde{F}}_2+{\tilde{U}}_d).\end{array}$$

$${\hat{A}}_{j,j+1}^0=-j_0\sum _{\nu =0}^{\infty }\frac{{\tilde{U}}_d^{\nu }({\hat{n}}_j-{\hat{n}}_{j+1}{)}^{\nu }}{\nu !}\frac{{\partial }^{\nu }}{\partial {\tilde{F}}_2^{\nu }}{J}_0({\tilde{F}}_1,{\tilde{F}}_2),$$

$$\begin{array}{r}{\hat{H}}_{eff}^{int}=V_e\sum _j{\hat{n}}_j{\hat{n}}_{j+1}+\frac{U_e}{2}\sum _j{\hat{n}}_j({\hat{n}}_j-1),\end{array}$$

$$\begin{array}{r}{V}_e=\frac{4j_0^2}{{\omega }^2}\sum _{n\ne 0}\left[-\frac{U_d}{n}\frac{\partial }{\partial {\tilde{F}}_2}+\frac{U_0}{2n^2}\right]J_n^2({\tilde{F}}_1,{\tilde{F}}_2),\end{array}$$

$$\begin{array}{r}{U}_e=U_0-2V_e,\end{array}$$

$${\hat{H}}_c=\sum _j{T}_c{c}_k^{†}{c}_k^{†}{c}_j{c}_j$$

$$T_c=-\frac{j_0^2}{{\omega }^2}\sum _{n\ne 0}{J}_n({\tilde{F}}_1,{\tilde{F}}_2)\left[\frac{U_d}{n}\frac{\partial }{\partial {\tilde{F}}_2}+\frac{U_0}{4n^2}\right]J_{-n}({\tilde{F}}_1,{\tilde{F}}_2).$$

$${\hat{H}}_s=\sum _j{T}_s{c}_{j+1}^{†}{c}_{j-1}^{†}{c}_j{c}_j$$

$${\hat{H}}_a=\sum _j{T}_c{c}_j^{†}\left({\hat{n}}_j-4{\hat{n}}_{j+1}+{\hat{n}}_{j+2}\right)c_{j+2}.$$

$$\begin{array}{rl}{T}_c{c}_j^{†}& \left({\hat{n}}_j-4{\hat{n}}_{j+1}+{\hat{n}}_{j+2}\right)c_{j+2}|{\psi }_i⟩\\ & =T_c(n_j-4n_{j+1}+n_{j+2}-1)c_j^{†}{c}_{j+2}|{\psi }_i⟩,\end{array}$$

$$\begin{array}{rl}{H}_e& =\sum _{⟨jk⟩}{\hat{c}}_j^{†}{J}_{jk}^{(e)}{\hat{c}}_k+\sum _{⟨jk⟩}{T}_{jk}^{(e)}{\hat{c}}_j^{†}({\hat{n}}_j-{\hat{n}}_k){\hat{c}}_k\\ & +\frac{U_e}{2}\sum _j{\hat{n}}_j({\hat{n}}_j-1)+\frac{V_e}{2}\sum _{⟨jk⟩}{\hat{n}}_j{\hat{n}}_k,\end{array}$$

$$\frac{V_e}{U_e}\simeq \frac{4j_0^2}{{\omega }^2}\sum _{n\ne 0}\left[\frac{1}{2n^2}-\frac{1}{n}\frac{U_d}{U_0}\frac{\partial }{\partial {\tilde{F}}_2}\right]J_n^2({\tilde{F}}_1,{\tilde{F}}_2).$$

$$\begin{array}{r}\frac{V_e}{h_{(0,1)}^L}\approx \frac{4j_0{U}_d}{{\omega }^2{J}_0({\tilde{F}}_1,{\tilde{F}}_2)}\sum _{n\ne 0}\frac{1}{n}\frac{\partial }{\partial {\tilde{F}}_2}{J}_n^2({\tilde{F}}_1,{\tilde{F}}_2),\end{array}$$

$$\begin{array}{rl}\hat{H}=& -j_0\sum _j{\hat{c}}_j^{†}{\hat{c}}_{j+1}+T_e{\hat{c}}_j^{†}({\hat{n}}_j-{\hat{n}}_{j+1}){\hat{c}}_{j+1}+h.c.\\ & +\frac{U_e}{2}\sum _j{\hat{n}}_j({\hat{n}}_j-1)+V_e\sum _j{\hat{n}}_j{\hat{n}}_{j+1},\end{array}$$

$$\begin{array}{rl}{h}_{(0,1)}^L& =h_{(1,0)}^R=-j_0,\\ h_{(1,1)}^L& =-j_0+T_e,\\ h_{(1,1)}^R& =-j_0-T_e.\end{array}$$

Full text

I Introduction

The rapid progress in accurate control of matter by means of laser fields enables the realization of quantum simulations QS in various controllable systems, such as optical lattices with atomic OLBHM or molecular gases Dipolar , and trapped ions RTLGT . One of the central tools for the realization of desired Hamiltonians is Floquet engineering HF ; SD that gives access to the simulation of physical processes that would be prohibitively difficult to realize by static means NS ; Mag1 . With the help of rapid periodic shaking of optical lattices it is possible to tune the amplitudes of quantum tunnelling between nearest neighbours  DC and beyond that, to create long-range tunnelling processes SD and modify topological properties TC . Periodically modulated external magnetic fields can be employed in the vicinity of Feshbach resonances FR3 ; FR4 for the dynamical control of on-site interactions in optical lattices. The tunability of such interactions enabled the observation of the phase transition from superfluid to Mott insulator SIT , and it can also be exploited for the realisation of density assisted tunnelling processes EU ; CH .

So far, most activities in driven optical lattices are concerned with the control of tunnelling properties or on-site interactions. Only recently, also control of nearest neighbour interactions with great importance for the realization of entangled many-body ground states EBHM ; QP ; Dipolar ; LR1 ; Entangle or far-from-equilibrium quantum many-body dynamics cor ; LR2 came into the spotlight. In principle, effective interactions beyond on-site interactions can be realised in terms of regular lattice shaking, but those effective interactions will typically be much smaller than the intrinsic on-site interactions. In order to overcome this, we consider, in addition to a periodic shaking, also a periodic modulation of the on-site interaction. We derive a complete classification of effective processes in this system and find driving profiles that allow us to selectively enhance or suppress such processes. With the explicit example of a density-wave we demonstrate what type of physics BHM1 ; BHMSFMI ; CDW1 ; HI ; EBHM can be simulated in terms of such a doubly-modulated Bose-Hubbard system.

II The Doubly Modulated Model

A system of spinless bosons occupying the lowest band in an optical lattice can be described by the Bose Hubbard Hamiltonian

[TABLE]

where $\hat{c}_{j}(\hat{c}_{j})^{\dagger}$ annihilates (creates) a boson at site $j$, $\hat{n}_{j}=\hat{c}_{j}^{\dagger}\hat{c}_{j}$ is the occupation number operator. $U(t)$ is the amplitude of the on-site interaction, $F(t)$ is the amplitude of tilt or shaking, $\epsilon_{j}$ is given by the lattice geometry and $J_{jk}$ is the bare tunnelling rate between nearest neighbours, i.e. the summation $\sum_{\langle jk\rangle}$ is performed over adjacent sites only. We will consider the doubly modulated version of the Bose-Hubbard Model(BHM) EU , in which both on-site energy and the on-site interactions have a periodic time-dependence with period $T$. In practice, these time-dependences result from a periodic tilt or acceleration of the lattice, and a modulation of a magnetic field inducing a Feshbach resonance FR4 ; FR3 . We will consider a vanishing average on-site energy, i.e. $\int_{t}^{t+T}F(\tau)d\tau=0$, but allow for a finite average interaction $U_{0}=\int_{t}^{t+T}U(\tau)d\tau/T$, and denote the periodic deviation from the average by $U_{f}(t)$.

For a periodically driven system, Floquet theory Flo1 ; Flo2 ; Flo3 permits to describe the dynamics in terms of an effective time-independent Hamiltonian $\hat{H}_{eff}$ that induces the propagation over one period of driving. If the driving frequency $\omega$ exceeds all system energy scales the effective Hamiltonian $\hat{H}_{eff}$ can be computed in a perturbative manner Mag1 ; HF as sketched in Appendix A. Since strong driving is conflicting with this condition, it is necessary to consider the system Hamiltonian in a frame in which the strong driving amplitude does not contribute to the magnitude of the Hamiltonian, i.e. the energy scales of the system. It is necessary for the shaking $F(t)$ to be strong (i.e. $F(t)\gtrsim\omega$), but the modulation of the interaction does not need to exceed the driving frequency. A possible choice for the unitary defining a suitable frame would thus be $\exp(i\sum_{j}\theta_{j}(t)n_{j})$ with $\theta_{j}(t)=\epsilon_{j}\int_{0}^{t}d\tau\ F(\tau)$ as frequently employed in regular shaken lattices HF . The unitary

[TABLE]

with $\Gamma(t)=\int_{0}^{t}d\tau\ U_{f}(\tau)$ will however ease the analysis and result in more compact expressions. Strictly speaking, the frame defined by $\hat{U}_{I}(t)$ would enable a treatment even if the amplitude of the time-dependent interaction exceeds the driving frequency $\omega$. Since this, however would result in new effective many-body interactions, we will limit our analysis to sufficiently weak interactions, and apply a corresponding series expansion where appropriate.

The transformed Hamiltonian $\tilde{H}(t)=U_{I}\hat{H}(t)U^{\dagger}_{I}-iU_{I}\dot{U}^{\dagger}_{I}$ in the frame defined by $U_{I}$ reads

[TABLE]

and admits a perturbative expansion even if the amplitude $F(t)$ exceeds the driving frequency $\omega=2\pi/T$. The system’s rich physics is contained in the time-dependent operators TC ; SIT

[TABLE]

defined in terms of the displacements $\textbf{r}_{jk}=\textbf{r}_{j}-\textbf{r}_{k}$, and the phase $\chi_{{jk}}(t)=\theta_{k}(t)-\theta_{j}(t)$, that a particle acquires when tunnelling from site $k$ to site $j$.

At the first glance, $\tilde{H}(t)$ given in Eq. (3) looks like a standard Bose-Hubbard Hamiltonian. The first term of $\tilde{H}(t)$, however, does not describe only single-particle tunnelling from site $k$ to site $j$, but the term $\hat{A}_{jk}^{0}$ may be interpreted as density assisted tunnelling, i.e. the amplitude for a tunnelling process between sites $j$ and $k$ depends explicitly on the particle number difference between those two sites. The physics described by $\tilde{H}(t)$ is thus substantially richer than that of the bare Bose-Hubbard Hamiltonian.

III Effective Processes

Since the perturbative construction of the effective Hamiltonian is entirely analytic, it can be done without specifying an explicit time-dependence of the system. We will therefore leave the driving profiles $\Gamma(t)$ and $\theta_{j}(t)$ unspecified and determine specific time-dependencies only later in Sec. III.

Standard perturbative expansions provide the effective Hamiltonian as series in powers of $J_{jk}/\omega$ and $U_{0}/\omega$, and an expansion including terms up to $1/\omega^{2}$ is appropriate for our purposes. Within this approximation, the effective Hamiltonian still contains general dependence on the amplitudes of shaking and the modulation of on-site interactions. Since the latter needs to be assumed to be sufficiently small, we can perform the corresponding series expansion and neglect terms beyond order $1/\omega^{2}$. Wherever, such an expansion is not helpful, we will keep the general dependence, but bear in mind that sufficiently weak interactions are assumed.

The lowest order ($1/\omega^{0}$) effective Hamiltonian is a regular Bose-Hubbard Hamiltonian with modified tunnelling and interaction rate HF ; SIT ; Colloquium . In first order, there is density assisted tunnelling (in terms of the operator $c_{j}^{\dagger}(\hat{n}_{j}-\hat{n}_{k})c_{k}$) as a physically distinct new process. Of particular interest for our present purposes is the second order ($1/\omega^{2}$) nearest neighbour interaction $\hat{n}_{j}\hat{n}_{k}$ resulting from a particle tunnelling from site $j$ to site $k$, interacting with the particles at site $k$ and tunnelling back to site $j$, as reflected by the commutation relation

[TABLE]

As discussed in more detail in Appendix. A the contributions to the effective Hamiltonian that are of purely interacting nature, result in

[TABLE]

where the summation $\sum_{\langle jk\rangle}$ is performed over neighbouring sites only with the interaction constants

[TABLE]

defined in terms of the Fourier components

[TABLE]

that depend on the driving functions $\Gamma(t)$ and $\theta_{j}(t)$ (via $\chi_{jk}(t)=\theta_{k}(t)-\theta_{j}(t)$).

In addition to these interaction terms and processes contained in the bare system, the effective Hamiltonian also contains additional tunnelling processes. The exact form of these processes and their amplitudes depend on the driving profiles similarly to the dependence of the interaction terms given in Eqs.(7) and (8). The general case is discussed in more detail in Appendix A.

In the following, we will exploit the freedom to chose driving profiles in order to realise an effective system with desired properties in a one-dimensional chain. The profiles can be chosen from a continuous set, but we will see that a monochromatically modulated on-site energy

[TABLE]

and a bi-chromatic shaking resulting in

[TABLE]

with $\epsilon_{j}=-j$ is all that is required. In order to simplify notation, we will also use the dimensionless quantities $\tilde{U}_{d}=U_{d}/\omega$, $\tilde{F}_{1}=F_{1}/\omega$ and $\tilde{F}_{2}=F_{2}/\omega$ in the following.

III.1 Elementary processes

Before discussing means to control the system dynamics in terms of the driving parameters $U_{d}$, $F_{1}$ and $F_{2}$ we will first give an overview over all processes found in the effective Hamiltonian up to second order. Since we employ bi-chromatic driving, the parameters in the effective Hamiltonian can no longer be expressed in terms of regular Bessel functions $J_{n}(x)$ of a single variable, but we will encounter two-dimensional Bessel functions; those can be expressed in terms of one-dimensional Bessel functions via the relation $J_{n}(x,y)=\sum_{s=-\infty}^{+\infty}J_{n-2s}(x)J_{s}(y)$, and for all practical purposes the summation can be limited to a finite range.

III.1.1 Nearest-neighbour density assisted tunnelling

The Hamiltonian

[TABLE]

given as the temporal average of the tunnelling term in Eq.4, depends on the driving parameters via the relation

[TABLE]

Because of the term $\hat{n}_{j}-\hat{n}_{k}$, the prefactors $\hat{A}^{0}_{jk}$ should not be taken literally as a rate but a proper rate can be specified for given initial and final states only EU . In a one dimensional system with the initial Fock state $\ket{\psi_{i}}=\ket{n_{1}\ldots n_{N}}$ one obtains

[TABLE]

with the prefactor

[TABLE]

that can be interpreted as the rate for a particle to tunnel from site $j+1$ to $j$, i.e. to the left. The equivalent rate for a particle to tunnel from site $j$ to $j+1$ reads

[TABLE]

States with low occupations per site will be particularly relevant later-on, and in those cases, the above rates reduce to the simple expressions

[TABLE]

In regular tunnelling processes described by the bare Bose-Hubbard Hamiltonian, the rates to tunnel left or right can only differ by a phase factor, but if the factors $h_{j,k}^{L/R}$ are interpreted as a rate, then one needs to accept that rates to tunnel left or right can also differ in their magnitude, as it generally is the case for $h_{(2,0)}^{R}$ and $h_{(0,2)}$EU .

In order to highlight the explicit dependence of particle densities on the tunnelling processes, it is appropriate to expand $\hat{A}^{0}_{jk}$ (Eq. 12) in powers of $U_{d}$. The expansion

[TABLE]

then yields explicit rates in terms of derivatives of the scalar Bessel function $J_{0}(\tilde{F}_{1},\tilde{F}_{2})$. The lowest order term describes regular tunnelling, and the higher order terms correspond to density assisted tunnelling processes.

III.1.2 Nearest-neighbour interactions

In the one dimensional chain the interaction terms of the effective Hamiltonian given in Eq. (6) reduces to

[TABLE]

with the nearest neighbour interaction constant

[TABLE]

and the on-site interaction constant

[TABLE]

which differs from the constant $U_{0}$ of the static system by the contribution to the two nearest neighbour interaction processes.

III.1.3 Co-tunnelling

Co-tunnelling processes with two particles tunnelling from one site $j$ to the same neighbour site $k$, can be described by

[TABLE]

with $k=j+1$ or $k=j-1$. The tunnelling rate can be expressed in Bessel functions as

[TABLE]

In contrast to the density assisted tunnelling discussed above, there is no preferred direction, and the rate for particles to tunnel to the left coincides with the rate for particles to tunnel to the right.

III.1.4 Split-tunnelling

Another emerging process appearing in second order can be interpreted as split-tunnelling, where particles on the same site $j$ tunnel to two distinct adjacent sites, i.e. site $j-1$ and $j+1$ in the one-dimensional case. This is captured by the Hamiltonian

[TABLE]

with the split tunnelling rate $T_{s}$ that is determined by the interaction constant $V_{e}$ given in Eq. 19 by the relation $T_{s}=-V_{e}/8$.

III.1.5 Next-nearest-neighbour density assisted tunnelling

In second order ($1/\omega^{2}$), there is nearest-neighbour density assisted tunnelling described by the Hamiltonian

[TABLE]

The amplitude for this process coincides indeed with the amplitude for co-tunnelling given above in Eq. (22)

Similarly to the discussion in Sec. III.1.1, the amplitude $T_{c}$ should not be taken literally as a rate, but rates for tunnelling processes depend on the occupation of all involved sites. For the initial Fock state $\ket{\psi_{i}}=\ket{n_{1}\ldots n_{N}}$, one obtains the rates

[TABLE]

such that $T_{c}(n_{j}-4n_{j+1}+n_{j+2}-1)$ can be interpreted as a rate, and, in contrast to the nearest-neighbour tunnelling discussed above, the rate for tunnelling to the right coincides with the rate to tunnel to the left.

III.2 Control

We will be interested in the dynamics induced by nearest neighbour interactions and density assisted tunnelling in addition to the regular tunnelling and on-site interaction captured by the Bose-Hubbard Hamiltonian.

We focus on the resulting effective model

[TABLE]

which features the rich physics of different strongly correlated phases of matter, as we show in Sec. IV.

In order to arrive at this model, we need to ensure that split-tunnelling, co-tunnelling and next-nearest neighbour tunnelling vanish or become negligible. Since the rate $T_{s}$ of split-tunnelling depends on the strength $V_{e}$ via the relation $T_{s}=-V_{e}/8$, it is not possible to find driving patterns resulting in vanishing $T_{s}$ but finite $V_{e}$. The factor $1/8$, however, ensures that the rate of split tunnelling is an order of magnitude smaller than the nearest neighbour interaction strength, such that we neglect it in the following.

In order to have a substantial impact of the nearest-neighbour interaction, it is important that it is not outweighed by the on-site interaction. Consistent with the perturbative expansion, the ratio between the interaction strengths $U_{e}$ and $V_{e}$ following Eqs. (19) and (20) can be approximated as

[TABLE]

Generically, the nearest-neighbour interaction will thus be substantially weaker than the on-site interaction due to the prefactor ${j_{0}^{2}}/{\omega^{2}}$. It is however possible to chose the static interaction constant $U_{0}$ to be much smaller than the amplitude $\omega U_{d}$ of the oscillating interaction, and thus arrive at a sizeable ratio between $V_{e}$ and $U_{e}$. On the other hand, it is necessary to choose a non-zero $F_{1}$ to make sure the second term in Eq. 27 remains finite, as discussed in more detail in Appendix B.

A finite driving amplitude $F_{2}$ is not necessary for the realisation of nearest neighbour interactions, but it is required for the realisation of an extended Bose Hubbard model with the parity symmetry breaking with $h_{(1,1)}^{R}\neq h_{(1,1)}^{L}$. In order to arrive at strongly correlated ground states, the ratio between the nearest neighbour interaction and the nearest neighbour tunnelling rates needs to be sufficiently large. Following Eqs. 16 and 19 the ratio of these two rates can be approximated as

[TABLE]

if $U_{0}$ is much smaller than $U_{d}$. The factor $j_{0}U_{d}/\omega^{2}$ is generically small, but the terms involving Bessel functions can be tuned within a wide range, and especially a substantial ratio can be realised with a small value of $J_{0}(\tilde{F}_{1},\tilde{F}_{2})$.

With these general observations in mind, we can now discuss the amplitudes for the effective processes and their dependence on the tuneable parameters $U_{d}$, $F_{1}$ and $F_{2}$ defined above in Eqs. (9) and (10). Fig. 1 depicts the amplitudes (in terms of the bare tunnelling rate $j_{0}=1$) for several effective processes for driving frequency $\omega=6j_{0}$, a static on-site interaction $U_{0}=0.3j_{0}$ and an AC on-site interaction amplitude $U_{d}=-4j_{0}$. The amplitude $F_{2}$ adopts the value $F_{2}=22j_{0}$ and $F_{1}$ ranges between [math] and $30j_{0}$. The rate $V_{e}$ of effective nearest-neighbour interaction (purple) is negligible for small values of $F_{1}$ as discussed above in Sec. III.2. In the parameter window $5j_{0}\lesssim F_{1}\lesssim 22j_{0}$, however, the nearest-neighbour interaction is by no means negligible as compared to the other processes; in particular around the sweet-spot $F_{1}\simeq 14j_{0}$ the rate of nearest-neighbour interaction is approximately twice as large as the prefactor of on-site interaction.

In addition, tunnelling processes can be suppressed with suitable choices of $F_{1}$. In particular for $F_{1}\simeq 15.5j_{0}$, the tunnelling process with rate $h_{(1,1)}^{R}$ is forbidden, but the remaining channels still exist with the ratio $|V_{e}/h_{(1,0)}^{L/R}|\approx 5$. Another interesting point is around $F_{1}\approx 14$ such that $\frac{\partial}{\partial\tilde{F}_{2}}J_{0}(\tilde{F}_{1},\tilde{F}_{2})\approx 0$, where three density-dependent tunnelling rates coincide at $-j_{0}J_{0}(\tilde{F}_{1},\tilde{F}_{2})$ and the parity symmetry remains preserved. The co-tunnelling rate $T_{c}$, depicted in black in Fig. 1, is negligible for all values of $F_{1}$.

As the main result of our paper we have thus shown that the combination of shaking and modulated on-site interactions significantly extends the achievable effective Hamiltonian, and the individual terms can be enhanced or suppressed simply by tuning the driving parameters. After neglecting suppressible undesired processes, the effective Hamiltonian up to the order of $\omega^{-2}$ can be written as Eq. 26.

IV Quantum phases

Having discussed the possibility to tune the rates of various processes, we can finally give an idea of the type of physics that can be explored with the processes described above. To this end, we consider the Hamiltonian

[TABLE]

comprised of a bare Bose-Hubbard Hamiltonian, and additional density assisted tunnelling and nearest-neighbour interaction.

In order to understand the ground state properties of this Hamiltonian, it is helpful to consider the dynamics induced by $\hat{H}$ with the Mott insulator state $|\Psi\rangle=\bigotimes_{i=1}^{N}c_{i}^{\dagger}|0\rangle$ as initial state. The transition amplitudes, defined in Eq. 16 – expressed in terms of $j_{0}$ and $T_{e}$ – read

[TABLE]

In particular for $T_{e}\simeq j_{0}$, the tunneling rates $h_{(1,1)}^{L}$ vanishes whereas the amplitude $h_{(1,1)}^{R}$ becomes particularly large. The system will thus evolve from the Mott insulator state $|\Psi\rangle$ towards the state

[TABLE]

The resulting imbalanced doublon-holon defect pair, with a single doubly occupied site and adjacent un-occupied site on its left, leads to violation of parity symmetry. The existence of defect pairs reduces the nearest-neighbour interaction in the system Hamiltonian (Eq. (29)), but it increases the on-site interaction energy. One might thus expect that the system ground state shows signatures of a density wave for a suitable balance of the several processes in the system.

We introduce the average number of two types of defect pairs in 02 or 20 order in the ground state as $N_{02/20}$ to assess it more quantitatively. The difference

[TABLE]

for defect pairs indicates a preferred direction of tunnelling.

The creation of a pronounced imbalance between the two different types of defect pairs around the the tunnelling amplitude $T_{e}\simeq j_{0}$ can be seen in Fig. 2, where $\Delta N$ is depicted based on numerically exact diagonalisation of the extended Bose Hubbard Hamiltonian (Eq. (29)) for $8$ particles on $8$ lattice sites with periodic boundary conditions to reduce finite size effect. Given the strong on-site interaction with $U_{e}=7j_{0}$, the system is in a Mott insulating state for sufficiently low values of $T_{e}$ and $V_{e}$. The blue area shows the expected preference of the system for defects resulting for tunnelling to the right, and the number of defect pairs is expected to grow with increasing nearest-neighbour interactions.

In order to analyse the phase transition from the Mott insulator to a density wave more accurately, we focus on the order parameter EBHM

[TABLE]

where $\delta n_{j}=n_{j}-\langle n\rangle$ characterizes the deviation of the atom number on site $j$ from its average. Due to periodic boundary conditions, the system has translation symmetry; the order parameter $R_{DW}$ is independent of $j$, and the prefactor $(-1)^{r}$ ensures that $R_{DW}$ is non-negative. For a density wave, the order parameter remains finite for large distances, but it vanishes for the Mott insulator. In practice, we will chose $r=5$ which is the largest possible length within $10$ sites.

This order parameter is plotted as function of density assisted tunnelling rate $T_{e}$ and nearest neighbour interaction constant $V_{e}$ in Fig. 3. Given the finite system size, there is not a sharp transition from Mott insulator (blue) to a density wave (red), but the threshold values $R_{DW}=0.25\pm 0.05$ and $R_{DW}=0.01\pm 0.005$ depicted as $V_{d}$ (black) and $V_{c}$ (green) provide an estimate for the transition. The establishment of a density wave is clearly supported by strong nearest neighbour interactions, whereas density-assisted tunnelling tends to destroy the density wave like behaviour. The top figure corresponds to an on-site interaction of strength $U_{e}=10j_{0}$ whereas the bottom figure corresponds to the weaker on-site interaction with $U_{e}=5j_{0}$. The comparison of both figures confirms that also strong on-site interactions tend to break density wave like behaviour as expected.

V Summary

Many experiments over the last decade highlighted the technological advances towards the control that is required for the realization of a quantum simulator. After those developments, it is now possible to not only control single-particle tunnelling processes, but to also introduce artificial interactions. This control over interactions pushes quantum simulations from proof-of-concept experiments on single-body physics to the exploration of many-body physics which would be strongly limited in terms of classical computational approaches.

The possibility to modulate both on-site energies and tunnelling processes as considered here offers quite some potential for the realisation of effective processes beyond single-body physics. Those processes help to realise correlated many-body phases like density waves as discussed here, but the range of accessibly many-body ground states is substantially wider than we would be able to explore here, e.g. the Haldane Bose insulator with unconventional string order HI . In addition, in two-dimensional lattices, where the interplay of different tunnelling processes can lead to the formation of states with topological order, there is an abundance of conceivable reactions of many-body system to different interaction effects. Exploring those in actual quantum simulations would allow us to actually exceed the limitation of our classical computers for questions of fundamental science.

In this article we focus on the engineered nearest neighbour interactions, a non-trivial extension of this work will be realizing notable longer ranged interactions. As seen from Eq. 5, the induced interactions are not limited to nearest neighbouring sites, as long as tunnellings beyond nearest neighbours are permitted, e.g. for quantum gas systems in shallow lattices. Consequently more exotic strongly correlated phases are probably within reachNS .

As we demonstrated here, the driving schemes to realize such quantum simulations do not require sophisticated shaped pulses, but simple bi-chromatic driving already gives us access to a broad range of physical processes.

VI Acknowledgement

We are indebted to Frederic Sauvage for stimulating discussions.

Appendix A Effective Hamiltonian

Physical origins and derivations of effective processes illustrated in Sec. III.1 will be discussed here. In the high-frequency regime in which the driving frequency $\omega$ exceeds all other relevant scales in the system Hamiltonian, the time evolution can be approximately captured by the time-independent effective Hamiltonian which can be obtained perturbatively in powers of $1/\omega$. For the present purposes, terms up to $1/\omega^{2}$ will be sufficient. In first order ($\propto 1/\omega$), two-step tunnelling results from a particle tunnelling from site $k$ to site $m$ and moving to site $j$ immediately after, as reflected by the commutation relation

[TABLE]

The process in which the particle tunnels back to site $k$ (i.e. $j=k$) results in a modification of the on-site energy.

In addition, two-step density assisted tunnelling($1/\omega^{2}$) appear resulting from a particle tunnelling from site $k$ to site $m$, interacting with the particles at site $m$ and keep moving to site $j$, as reflected by the commutation relation

[TABLE]

In fact, besides interacting with particle at site $m$, the interaction can happen at site $j,k$ as well, therefore the occupation number dependence also takes additional form of $c_{j}^{\dagger}n_{k}c_{k},c_{j}^{\dagger}n_{j}c_{k}$. Finally split-tunnelling results from two particle first interacting with other particles at the same site $m$, then tunnelling to distinct nearest neighbour sites $j$ and $k$ respectively as

[TABLE]

The co-tunnellings is almost the same as split-tunnelling apart from that two particles ending at same site as $c_{j}^{\dagger}c_{j}^{\dagger}c_{m}c_{m}$. Finally the engineered interactions have been discussed in Sec. III. In order to get the amplitudes for each process, in practice we use Magnus expansion to compute the lowest three orders effective Hamiltonians, and within each of them we neglect all terms with higher order than $1/\omega^{2}$.

Lowest order effective Hamiltonian

The lowest effective Hamiltonian is the temporal average of the time-dependent Hamiltonian (Eq. (3)), and reads

[TABLE]

where the density assisted tunnelling can be expressed as

[TABLE]

In the one dimensional case with the drivings defined in Eq. 9 10, it reduces to

[TABLE]

and the dependence on density difference in terms of $1/\omega$ can be derived from Eq. 17.

First order effective Hamiltonian

The construction of the first order effective Hamiltonian is eased by expressing the Hamiltonian $H(t)$ in terms of its Fourier components $\hat{H}_{n}=\sum_{\langle jk\rangle}\hat{c}_{j}^{\dagger}\hat{A}_{jk}^{n}\hat{c}_{k},$ with

[TABLE]

for $n\neq 0$. Up to first order, the operators $\hat{A}_{jk}^{n}$ can be approximated as $\hat{A}_{jk}^{n}\approx g_{jk}^{n}+(\hat{n}_{j}-\hat{n}_{k})t_{jk}^{n}$ with the Fourier components

[TABLE]

is sufficient. For the one dimensional case, we have

[TABLE]

where $k=j\pm 1$ for nearest neighbour tunnelling. The first order effective Hamiltonian including terms up to second order thus reads

[TABLE]

with the coefficients

[TABLE]

The first term of Eq. 40, $\sum_{\langle\langle jk\rangle\rangle,j=k}c_{j}^{\dagger}B_{jk}c_{k}$ indicates the two-step tunnelling and the modification to on-site energies with rates $B_{jk}$. As proven in the Appendix.A in  SD , both of them vanish regardless of the drivings for the one dimensional case. The second expression of Eq. 40 with rate $T_{jmk}$ illustrates the co-tunnelling and split-tunnelling processes. The operator $\sum_{\langle\langle jk\rangle\rangle}c_{j}^{\dagger}\hat{S}_{jk}c_{k}$ captures the two-step density-assisted tunnelling. The last two give us the induced nearest neighbour interactions and the modification to the on-site interaction.

Second order effective Hamiltonian

The second order processes ($1/\omega^{2}$) also exist in the second order effective Hamiltonian $H_{eff}^{2}$,which can be expressed HF as

[TABLE]

Since $\omega^{2}$ appears in the dominator, all contributions with terms $\hat{H}_{n}\propto 1/\omega$ are negligible and we can make the following approximation

[TABLE]

Such an approximation leads to the crucial simplification to compute $\hat{H}_{eff}^{2}$ by noticing that the commutator appearing in Eq. 43

[TABLE]

with $B_{jk}^{pq}=\sum_{l}g_{jl}^{p}g_{lk}^{q}-g_{jl}^{q}g_{lk}^{p}$, vanishes in the one dimensional case for arbitrary $p,q$. For instance, when $k=j+2$, the tunnelling rate reduces to

[TABLE]

which turns to be zero according to Eq. 39. Therefore the non-vanishing contribution to $\hat{H}^{2}_{eff}$ is merely from the presence of the constant on-site interaction $U_{0}$ as

[TABLE]

with the interacting Hamiltonian $\hat{H}_{int}=\frac{U_{0}}{2}\sum_{j}\hat{n}_{j}(\hat{n}_{j}-1)$. No new process will be introduced, but the amplitudes given in Eq. 42 will be modified as

[TABLE]

with

[TABLE]

where $S_{jmk}^{1/2},T_{jmk}$ are defined in Eq. 42. For the one dimensional case, those reduce to the co-tunnelling rate

[TABLE]

the split tunnelling rate

[TABLE]

and the interaction constant. $V_{e}=-8T_{s}$ According to Eq. 39, one can express in terms of Bessel functions and retrieve expression given in Sec. III.

Appendix B Special case with zero $F_{1}$

In order to achieve a notable ratio of $V_{e}/U_{e}$ as considered in Sec. III.2, non-vanishing $F_{1}$ is demanded, otherwise the contribution proportional to $U_{d}$ to the nearest neighbour interactions(Eq. 19) vanishes. According to Eq. 19, the first term involving $U_{d}$ reads

[TABLE]

All contributions with odd $n$ vanishes since the Bessel function $J_{n}(0,\tilde{F}_{2})$ is zero for all values of $\tilde{F}_{2}$ according to the symmetry relation of the two dimensional Bessel function Bessel . $V_{d}$ reduces to

[TABLE]

where the two dimensional Bessel function can be rewritten as a one dimensional Bessel function as $J_{-2k}(0,\tilde{F}_{2})=J_{k}(\tilde{F}_{2})$. Finally the property $J_{-k}(\tilde{F}_{2})=(-1)^{k}J_{k}(\tilde{F}_{2})$Bessel leads to the expression

[TABLE]

where the numerator $1-(-1)^{2k}$ cancels. Therefore a finite $\tilde{F}_{1}$ is crucial to generate $V_{d}$ and furthermore the notable nearest neighbour interactions.