Theory of interacting fermions in shaken square optical lattice

TL;DR

This paper develops a theoretical framework for weakly interacting fermions in shaken square optical lattices, revealing complex Fermi surfaces and unconventional pairing symmetries induced by orbital hybridization and periodic driving.

Contribution

It introduces a four-band effective Hamiltonian for shaken lattices, analytically derives momentum-dependent interactions, and predicts an $s+d$-wave pairing phase in driven fermionic systems.

Findings

Hybridized $s$-band develops multiple minima.

Effective interactions gain momentum dependence.

Predicted $s+d$-wave pairing symmetry.

Abstract

We develop a theory of weakly interacting fermionic atoms in shaken optical lattices based on the orbital mixing in the presence of time-periodic modulations. Specifically, we focus on fermionic atoms in circularly shaken square lattice with near resonance frequencies, i.e., tuned close to the energy separation between $s$-band and the $p$-bands. First, we derive a time-independent four-band effective Hamiltonian in the non-interacting limit. Diagonalization of the effective Hamiltonian yields a quasi-energy spectrum consistent with the full numerical Floquet solution that includes all higher bands. In particular, we find that the hybridized $s$-band develops multiple minima and therefore non-trivial Fermi surfaces at different fillings. We then obtain the effective interactions for atoms in the hybridized $s$-band analytically and show that they acquire momentum dependence on the Fermi surface even though the bare interaction is contact-like. We apply the theory to find the phase diagram of fermions with weak attractive interactions and demonstrate that the pairing symmetry is $s+d$-wave. Our theory is valid for a range of shaking frequencies near resonance, and it can be generalized to other phases of interacting fermions in shaken lattices.

Tables

Table 1. Table 1: Numerical values of the parameters in the effective Hamiltonian H eff subscript 𝐻 eff H_{\mathrm{eff}} for shaken square lattice for different lattice depth V 0 subscript 𝑉 0 V_{0} in units of recoil energy E R subscript 𝐸 𝑅 E_{R} . e n subscript 𝑒 𝑛 e_{n} and t n subscript 𝑡 𝑛 t_{n} are onsite energies and nearest neighbor hopping consistent with Ref.  Li and Liu ( 2016 ) . d 0 subscript 𝑑 0 d_{0} and d 1 subscript 𝑑 1 d_{1} are the onsite and nearest neighbor inter-orbital coupling strength, respectively.

$V_0$	$e_s$	$e_p$	$t_s$	$t_p$	$d_0$	$d_1$
3	-0.06	2.57	-0.11	0.50	2.07	-0.44
5	-0.56	2.76	-0.07	0.42	2.63	-0.28
10	-2.12	2.94	-0.02	0.24	3.47	-0.07
20	-5.80	1.98	-0.00	0.06	4.37	-0.01

Equations

$$V_{\mathrm{lat}}(\mathbf{x})=V_0\left[\cos (\frac{2\pi x}{{\lambda }_L})+\cos (\frac{2\pi y}{{\lambda }_L})\right].$$

$${\mathbf{x}}_0(\tau )={\mathbf{x}}_0(\tau +T).$$

$${\mathbf{x}}_0(\tau )=s_0[\sin \omega \tau ,\sin (\omega \tau +\vartheta )],$$

$${\mathcal{L}}_0$$

$$H_0(\mathbf{x},\tau )=\frac{{\mathbf{p}}^2}{2m}+V_{\mathrm{lat}}(\mathbf{x}+{\mathbf{x}}_0(\tau ))$$

$$\left\{{\psi }_{\sigma }(\mathbf{x},\tau ),{\psi }_{{\sigma }^{\prime }}^{†}({\mathbf{x}}^{\prime },\tau )\right\}={\delta }_{\sigma {\sigma }^{\prime }}\delta (\mathbf{x}-{\mathbf{x}}^{\prime }).$$

$$U(\mathbf{x}-{\mathbf{x}}^{\prime },\tau -{\tau }^{\prime })=g\delta (\mathbf{x}-{\mathbf{x}}^{\prime })\delta (\tau -{\tau }^{\prime }).$$

$${\mathcal{L}}_1$$

$$H_0^{\prime }=\frac{{\mathbf{p}}^2}{2m}+V_{\text{lat}}(\mathbf{x})+{\dot{\mathbf{x}}}_0(\tau )\cdot \mathbf{p}$$

$$H_0^{\prime \prime }=\frac{[\mathbf{p}+m{\dot{\mathbf{x}}}_0(\tau ){]}^2}{2m}+V_{\mathrm{lat}}(\mathbf{x}).$$

$$\mathbf{A}=-m{\dot{\mathbf{x}}}_0$$

$$\mathcal{U}(T)$$

$${\psi }_{\mathbf{k}}(\mathbf{x},\tau )=e^{i\mathbf{k}\cdot \mathbf{x}}\sum _{\mathbf{G}}{\Psi }_{\mathbf{k}}(\mathbf{G},\tau ){\mathrm{e}}^{\mathrm{i}\mathbf{G}\cdot \mathbf{x}}$$

$$H_0^{\prime \prime }(\mathbf{\ell },{\mathbf{\ell }}^{\prime };\tau )={\left[\mathbf{k}+2\mathbf{\ell }+\mathbf{A}(\tau )\right]}^2{\delta }_{\mathbf{\ell },{\mathbf{\ell }}^{\prime }}+V_{\mathbf{\ell },{\mathbf{\ell }}^{\prime }},$$

$$\mathbf{A}(\tau )=\beta \left[\cos \tau ,\cos (\tau +\vartheta )\right]$$

$$\beta =(\omega /E_R)(s_0{k}_L).$$

$$\mathcal{U}(T)=\prod _{{\tau }_i=0}^{2\pi }\exp \left\{-\frac{\mathrm{i}}{\omega }{H}_0^{\prime \prime }(\mathbf{\ell },{\mathbf{\ell }}^{\prime };{\tau }_i)\right\}.$$

$${\mathcal{E}}_n(\mathbf{k})=-\frac{1}{T}\mathrm{Im}\log {\Lambda }_n(\mathbf{k})$$

$$\Omega \equiv \max E_p-\min E_s,$$

$$\psi (\mathbf{x},\tau )=\sum _{i,n}{\psi }_{i,n}(\tau ){\mathcal{W}}_n(\mathbf{x}-{\mathbf{R}}_i).$$

$$\mathrm{i}{\partial }_{\tau }{\psi }_{i,n}(\tau )=\sum _{i^{\prime },n^{\prime }}[H_0^{\prime }(\tau ){]}_{i{i}^{\prime },n{n}^{\prime }}{\psi }_{i^{\prime },n^{\prime }}(\tau ),$$

$$[H_{\mathrm{static}}{]}_{i{i}^{\prime }}$$

$$+\sum _{\pm }{\delta }_{i,i^{\prime }\pm \hat{x}}\mathrm{Diag}(t_s,t_p,t_s,t_p)$$

$$+\sum _{\pm }{\delta }_{i,i^{\prime }\pm \hat{y}}\mathrm{Diag}(t_s,t_s,t_p,t_p).$$

$$e_n=\int d\mathbf{x}{\mathcal{W}}_n^{\ast }(\mathbf{x})H_0^s(\mathbf{x}){\mathcal{W}}_n(\mathbf{x}),$$

$$t_s$$

$$t_p$$

$$[V(\tau ){]}_{i{i}^{\prime }}$$

$$+\sum _{\pm }{\delta }_{i,i^{\prime }\pm \hat{x}}{a}_x(\tau )V_1^x$$

$$+\sum _{\pm }{\delta }_{i,i^{\prime }\pm \hat{y}}{a}_y(\tau )V_1^y.$$

$$a_x(\tau )=s_0\omega \cos (\omega \tau ),a_y(\tau )=s_0\omega \cos (\omega \tau +\vartheta ).$$

$$V_0^x$$

$$d_0=\int d\mathbf{x}{\mathcal{W}}_s^{\ast }(\mathbf{x}){\hat{p}}_x{\mathcal{W}}_{p_x}(\mathbf{x}).$$

Full text

I Introduction

Manipulating quantum many-body systems by time-periodic driving, known as Floquet engineering, has recently emerged as a powerful way to experimentally control the band structure of ultracold atoms in optical lattices. One particular useful approach is lattice shaking, namely moving the entire lattice along a certain prescribed periodic trajectory in space by tuning the phases of interfering laser beams that give rise to the optical lattice potential. With a properly designed shaking protocol, the time-dependent many-body system may resemble a time-independent system described by an effective Hamiltonian with desired properties.

Floquet engineering by lattice shaking has been successfully demonstrated in a number of experiments. Dynamical control of tunneling and band flattening in optical lattice was reported in Madison et al. (1998); Lignier et al. (2007). Superfluid-Mott insulator transition induced by shaking has also been observed Zenesini et al. (2009). Lattice shaking can give rise to artificial gauge fields for neutral atoms, complementary to other approaches based on for example the Raman coupling scheme or laser insisted tunneling (for review, see for instance Holthaus (2016); Goldman et al. (2014); Eckardt (2016)). Along this line, increasingly sophisticated shaking schemes were implemented to realize, for instance, frustrated magnetism on triangular latticeStruck et al. (2011), tunable gauge fields including staggered flux Struck et al. (2012); Hauke et al. (2012), the $XY$-spin model Struck et al. (2013), and the Haldane model on honeycomb lattice Jotzu et al. (2014). Shaken one-dimensional lattice with mixed Bloch bands and the resulting spectrum with double minima was used to simulate ferromagnetism Parker et al. (2013) and to test the universality relations near quantum phase transitions Clark et al. (2016).

Stimulated by these experiments, theoretical work has explored and clarified various regimes of shaking. Shaken lattice is intrinsically a multi-band problem, where the relevant energy scales are the band widths, energy gaps, and the shaking frequency $\omega$. For deep optical lattices and when $\omega$ is much larger than the width of the lowest $s$-band (but much smaller than the excitation gap), a single band picture is possible and the system can be described by a static Hamiltonian with an effective tunneling amplitude tunable in both magnitude and sign Eckardt et al. (2005). In this frequency regime, if the shaking protocol is adjusted to have certain symmetries, artificial gauge fields can be realized Sørensen et al. (2005); Hauke et al. (2012); Goldman and Dalibard (2014). An extensively studied shaken system is the two-dimensional hexagonal lattice in the tight binding limit. When the shaking frequency matches the energy difference between the sublattice sites, unconventional pairing is predicted to occur for fermions Zhang et al. (2015). Tuning the shaking frequency to match the band-gap at the Brillouin zone center leads to the so-called moat band Sedrakyan et al. (2015). For larger $\omega$, the higher Bloch bands have to be incorporated and band mixing may drastically modify the band structure. For example, when $\omega$ is tuned to half of the energy gap between $s$-band and $p$-band, two-photon processes dominate and give rise to topologically nontrivial bands in both one dimension Zheng and Zhai (2014) and two dimensions Zhang and Zhou (2014). Finally, when $\omega$ is on the order of the gap between the $s$- and $p$-bands, single photon processes dominate and lead to dispersion spectra with multiple minima implying interesting many-body phenomena Zheng et al. (2015).

In contrast to the single particle properties, understanding the interaction effects for fermions in shaken lattices remains a theoretical challenge. The problem is complicated by the time dependence and the hybridization of multiple Bloch bands. A key point is that the effective interactions between the atoms are modified by shaking and must be derived and analyzed along with the band structures. In this paper, we formulate an effective theory to capture the essential physics of interacting fermions on shaken lattices. For concreteness, we will focus on the two-dimensional square lattice with circular shaking and near resonance frequencies, i.e., with $\omega$ tuned close to the energy separation between $s$- and $p$-bands at the Brillouin zone center. Our theory is not restricted to the tight binding limit and various values of the optical lattice depth will be considered.

This paper is organized as follows. First we describe the problem and summarize our main results in Section II. The single particle spectrum for the shaken lattice is solved numerically in Sec. III using the standard Floquet analysis. The main goal here is to obtain the quasienergy spectrum accurately, which will be used to benchmark our approximation schemes. For example, the hybridized $s$-band is shown to develop four minima as opposed to a single minimum at the Brillouin zone center. By comparing the numerical spectrum with simple folding construction, we show that it is sufficient to keep only a few lowest Bloch bands. This observation motivates our subsequent analytical theory. In Sec. IV, we derive a four-band effective Hamiltonian $H_{\mathrm{eff}}$ for the single-particle shaken system by using the rotating wave approximation (RWA) that is consistent with the numerical Floquet solution. Diagonalizing $H_{\mathrm{eff}}$, we determine the non-trivial Fermi surface geometries of fermionic atoms populating the hybridized $s$-band as function of lattice filling fraction. Next, we take interactions into account in Sec.V and derive the effective interaction $V_{\mathrm{eff}}(\mathbf{k},\mathbf{k}^{\prime})$ for fermions on the Fermi surface of the hybridized $s$-band. In particular, we show that it develops interesting momentum dependence and no longer point-like. In Sec. VI, we apply the effective model to fermions with weak attractive interactions in circularly shaken square lattice. We investigate the pairing symmetry and transition temperature for different filling fractions. We show that the order parameter can have $s+d$-wave symmetry. The theoretical framework is applied in Sec. VII to study red-detuned near-resonance shaking, where the bands directly overlap. We conclude with remarks on the implications of our work in Sec. VIII.

II The problem and main results

We are interested in formulating a theory of interacting fermionic ultracold atoms in an optical lattice potential that varies periodically in time. Consider for example a square optical lattice given by the potential

[TABLE]

Here $\mathbf{x}=(x,y)$ is the space coordinate in two dimensions (2D), $\lambda_{L}$ is the wavelength of the lasers forming the optical lattice, $V_{0}$ is the lattice depth and we assume tight confinement in the $z$-direction. For a shaken optical lattice, the lattice potential $V_{\mathrm{lat}}(\mathbf{x})$ is replaced by $V_{\mathrm{lat}}(\mathbf{x}+\mathbf{x}_{0}(\tau))$ and becomes a periodic function of time $\tau$, since the origin of the lattice, $\mathbf{x}_{0}(\tau)$, moves along a prescribed loop in space and returns to its starting point after one shaking period $T$,

[TABLE]

Different choices of $\mathbf{x}_{0}(\tau)$ are referred to as different shaking protocols. For example, a general shaking protocol corresponds to the choice

[TABLE]

where $s_{0}$ is the amplitude of shaking, $\omega\equiv 2\pi/T$ is the shaking frequency and $\vartheta$ is the phase difference between $x$ and $y$ directions. In this paper we take $\vartheta=\pi/2$ which correspondd to circular shaking.

The many-body system of fermionic atoms loaded in such lattices is described by an action $\mathcal{S}=\int d\tau[\mathscr{L}_{0}+\mathscr{L}_{1}]$ where the single particle part of the Lagrangian $\mathscr{L}_{0}$ has the form

[TABLE]

Here $\sigma=\uparrow,\downarrow$ is the spin index and we take $\hbar=1$. For shaken optical lattices,

[TABLE]

with $\mathbf{p}=-i\nabla$. Following the discussion above regarding $V_{\mathrm{lat}}$, the single particle Hamiltonian $H_{0}$ is periodic both in time $H_{0}(\mathbf{x},\tau)=H_{0}(\mathbf{x},\tau+T)$ and in space $H_{0}(\mathbf{x},\tau)=H_{0}(\mathbf{x}+\mathbf{R}_{i},\tau)$ where $\mathbf{R}_{i}$ are the lattice vectors and $T$ is the shaking period. Note that the single-particle time-dependent Schrodinger equation can be obtained as $\mathrm{i}\partial_{\tau}\psi=H_{0}\psi$.

The fermionic field operators obey the equal-time anti-commutation relation

[TABLE]

For a dilute gas of ultracold alkali atoms, the interactions are local both in time and space such that the two particle interaction potential takes the form

[TABLE]

Here the interaction strength $g$ is related to the low energy s-wave scattering length. The interaction part of the Lagrangian, $\mathscr{L}_{1}$, can be written as

[TABLE]

Our strategy is to study the single particle physics of $\mathscr{L}_{0}$ first and then incorporate interaction $\mathscr{L}_{1}$ later. Even without $\mathscr{L}_{1}$, it is a challenge to analyze the time-dependent Hamiltonian $H_{0}(\mathbf{x},\tau)$ directly. To make progress, we recast $H_{0}$ in slightly different forms that are convenient for subsequent analytic or numeric treatments. This is achieved by performing a gauge transformation to co-moving frame, $\psi^{\prime}=e^{-i\mathbf{x}_{0}\cdot\mathbf{p}}\psi$, in which $H_{0}$ becomes

[TABLE]

with $\dot{\mathbf{x}}_{0}(\tau)=\partial\mathbf{x}_{0}/\partial\tau$. Note the result is valid for arbitrary shaking protocol. In this form, the driving appears as a time-dependent perturbation $\dot{\mathbf{x}}_{0}(\tau)\cdot\mathbf{p}$ to the static problem $\mathbf{p}^{2}/{2m}+V_{\text{lat}}(\mathbf{x})$. The last term in Eq. (9) can be combined with the first term by completing the square. The resulting time-dependent term proportional to $\dot{\mathbf{x}}^{2}_{0}$ can be removed via another gauge transformation Holthaus (2016). In this case, Hamiltonian becomes

[TABLE]

This form suggests that the lattice shaking is equivalent to the presence of a time-dependent vector potential

[TABLE]

and the corresponding force field is given by $\mathbf{E}(\tau)=-\partial\mathbf{A}(\tau)/\partial\tau=m\ddot{\mathbf{x}}_{0}$. The presence of a vector field $\mathbf{A}$ may drastically modify the band dispersions. We stress that the Hamiltonians given in Eqs. (5), (9) and (10) are equivalent. For the numerics in section III, we use Eq. (10) which is consistent with Ref. Holthaus (2016). For analytical manipulations in section IV and after, we use Eq. (9) instead to treat time dependence separately.

Our key results for the readers who wish to skip the technical details are summarized as follows. For near resonance shaking, the quasienergy spectrum can be captured by a static, four-band effective Hamiltonian $H_{\mathrm{eff}}$ given in Eq. (54) where the band mixing is described by the off-diagonal matrix elements. $H_{\mathrm{eff}}$ can be diagonalized analytically to yield the dispersion of the hybridized bands. For example, the hybridized $s$-band dispersion is given by $\epsilon_{1}(\mathbf{k})=\epsilon_{++}(\mathbf{k})$ in Eq. (70). The effective interactions for fermions on the hybridized $s$-band is given by $V_{\mathrm{eff}}(\mathbf{k},\mathbf{k}^{\prime})$ in Eq. (83). Based on the effective band structure and interactions, we solve the pairing problem of fermions with weak attractive interaction in the hybridized $s$-band to get the phase diagram and the symmetry of the order parameter as shown in Fig. 4.

III Quasienergy spectrum from Floquet analysis

We first numerically calculate the single particle spectrum of time periodic Hamiltonian $H^{\prime\prime}_{0}$ in Eq. (10). This problem is previously considered in Refs. Miao et al. (2015) and Po and Zhou (2015) for bosonic systems using Wannier expansion up to three and four orbitals, respectively. Here we adopt an approach based on Bloch expansion that can include all higher bands to desired numerical accuracy.

Floquet operator $\mathscr{U}(T)$ is defined as the time evolution operator over one shaking period $T$. It can be written as the following time ordered exponential,

[TABLE]

We expand the wave functions using the Bloch theorem

[TABLE]

where $\Psi_{\mathbf{k}}(\mathbf{G},\tau)$ are expansion coefficients for crystal momentum $\mathbf{k}\equiv(k_{x},k_{y})$ and reciprocal lattice vectors are given by $\mathbf{G}=\frac{2\pi}{\lambda_{L}}\boldsymbol{\ell}$ with intergers $\boldsymbol{\ell}=(\ell_{x},\ell_{y})$. In numerical calculations, we take $\ell_{x,y}=-N_{b},\cdots,N_{b}$ where the cutoff $N_{b}=4$ corresponds to the inclusion of $(2N_{b}+1)^{2}=81$ orbitals. We checked that increasing $N_{b}$ further does not change the results. To express the resulting Hamiltonian matrix in a simple form, we measure lattice momentum in units of the recoil momentum $k_{L}=\pi/\lambda_{L}$, energy and shaking frequency in units of the recoil energy $E_{R}=k_{L}^{2}/2m$: $k\rightarrow k/k_{L}$, $H_{0}\rightarrow H_{0}/E_{R}$ and $\omega/E_{R}\rightarrow\omega$. Then Eq. (10) becomes

[TABLE]

where the $\mathbf{k}$ dependence of $H^{\prime\prime}_{0}$ is suppressed for brevity. For the square optical lattice potential given in Eq. (1), the matrix elements can be calculated as $V_{\boldsymbol{\ell},\boldsymbol{\ell}^{\prime}}=(V_{0}/2E_{R})\left[\delta_{\ell_{x}^{\prime},\ell_{x}+1}+\delta_{\ell_{x}^{\prime},\ell_{x}-1}\right]\delta_{\ell_{y}^{\prime},\ell_{y}}+(x\leftrightarrow y)$ and the vector potential coming from shaking has the form

[TABLE]

where the dimensionless shaking strength is defined as

[TABLE]

The time ordered product is calculated by dividing the time evolution into many small slices $\{\tau_{i}\}$,

[TABLE]

For a given $\mathbf{k}$ and discrete $\tau_{i}$, we compute the corresponding matrix exponentials in Eq. (14). Then by taking their product, we obtain the Floquet operator. The eigenvalues of the Floquet operator, $\mathscr{U}(T)v_{n}(\mathbf{k})=\Lambda_{n}(\mathbf{k})v_{n}(\mathbf{k})$ has the form $\Lambda_{n}(\mathbf{k})=\mathrm{e}^{-\mathrm{i}T\mathscr{E}_{n}(\mathbf{k})}$ which defines the quasi-energy spectrum by the relation

[TABLE]

where $n$ is the Floquet band index. Note that the quasienergy spectrum is periodic, i.e., an energy level at $\mathscr{E}_{n}(\mathbf{k})$ is identical to $\mathscr{E}_{n}(\mathbf{k})+\omega$. The region $\mathscr{E}_{n}\in[0,\omega]$ is called the quasienergy Brillouin zone (QeBZ), analogous to the quasimomentum Brillouin zone of a solid.

The quasienergy band structure depends on two independent energy scales: the lattice depth in terms of recoil energy $V_{0}/E_{R}$ and the shaking frequency $\omega/E_{R}$ which couples different energy sectors with a strength on the order of $\beta$. We have explored several regimes in our numerics and found that one of the most interesting regimes is near resonance shaking, i.e., when $\omega$ is comparable to $\Omega$ defined by

[TABLE]

where the bottom of the $s$-band is $\min E_{s}$ and the top of two-fold degenerate the $p$-bands is $\max E_{p}$. For such frequencies, shaking strongly couples the $s$- and $p$-bands. This is illustrated in Fig. 1A, where $\omega$ is set to $\omega=1.05\Omega$. The bands in Fig. 1A are labeled using standard notation, e.g., the $s$, $p_{x}$, $p_{y}$ and $d$-bands, in line with Refs. Po and Zhou (2015); Miao et al. (2015). Note that the orbital symmetries like $s$, $p_{x,y}$ and $d$ are mixed in the shaken system so formally it is not possible to label bands by their symmetries. However, mixed orbitals still retain a dominant symmetry after shaking as long as shaking is not perturbatively strong. Based on this reasoning, we call the bands after shaking hybridized $s$-band, hybridized $p$-bands etc.

The gross features of the quasienergy spectrum of the shaken lattice can be captured by a simple folding construction: take the static bands given in Fig. 1A and fold them into the same QeBZ, $\mathscr{E}_{n}(\mathbf{k})\rightarrow\mathrm{mod}_{\omega}[E_{n}(\mathbf{k})]$, one then obtains a rough caricature of the quasienergy spectrum shown in Fig. 1B. We observe that, after folding, the $s$, $p_{x}$, $p_{y}$ and one of the $d$-bands, $d_{xy}$, come close in energy in regions indicated by the blue circles around the $\Gamma$ point. For the shaking frequency chosen here, even though there is no direct level crossing inside these circles and the $s$-band seems to be isolated, the mutual influence of these bands turns out to be important. There are however level crossings between higher bands, e.g., between the $p$- and $d$-bands.

Fig. 1C shows the full numerical solution of the quasienergy band structure. By comparing it to Fig. 1B, we can identify the relatively flat band within the small gray region as the hybridized $s$-band, which is plotted separately for clarity in Fig. 1D with a finer energy resolution. We notice that the level repulsion between the $s$-band and $p$-bands, particularly around the $\Gamma$ point, pushes up the bottom of the $s$-band located at the $\Gamma$ point. As a result, the minimum of the $s$-band moves away from $\Gamma$ to four $\mathbf{k}$ points on the diagonal $\Gamma-\mathrm{M}$ line. The four minima can also be seen in the full dispersion shown in Fig. 1E for three different shaking amplitudes. It clearly demonstrates that the $s$-band is strongly modified by lattice shaking. Similar effect has recently been observed in experiments with bosons in one-dimensional shaken lattice in Ref. Parker et al. (2013).

There are also other dramatic consequences of shaking to the quasienergy spectrum. For example, many of the level crossings of the $p$-bands and $d$-bands in Fig. 1B become avoided crossings in Fig. 1C. Moreover, the $d$-bands are also modified by coupling to higher bands such as the $f$-bands (not shown). Hereafter we are mainly interested in what happens to the $s$-band in the presence of non-perturbative shaking. For this purpose, our numerical findings here suggest that it is sufficient to consider the lowest four orbitals; $s$, $p_{x}$, $p_{y}$, and $d_{xy}$.

IV Effective Four-band model

In this section, we derive analytically an effective four-band Hamiltonian that can accurately describe the quasienergy spectrum up to the $d_{xy}$ band, consistent with the numerical Floquet solution. Based on previous section, a truncation up to four lowest Bloch bands is sufficient to understand the interplay of multiple orbitals in the shaking problem, as far as the modifications to the $s$-band and $p$-bands are concerned. The effective Hamiltonian will serve as the basis to study interaction effects in the next section.

To analyze the Hamiltonian $H^{\prime}_{0}$ in Eq. (9), we first carry out Wannier expansion for the wavefunction,

[TABLE]

Here $\mathcal{W}_{n}(\mathbf{x}-\mathbf{R}_{i})$ are Wannier functions localized at site $i$ for the static square optical lattice. We will truncate the orbital index $n$ and only keep four orbitals, $n=s,p_{x},p_{y},d_{xy}$. The truncation error can be assessed by comparing to the full numerical results of the previous section. In this basis, the Schrödinger equation becomes

[TABLE]

where $H^{\prime}_{0}$ is a $4\times 4$ matrix in orbital space. Since $H^{\prime}_{0}$ naturally splits into a static part $H_{\mathrm{static}}=\mathbf{p}^{2}/{2m}+V_{\text{lat}}(\mathbf{x})$ describing the unshaken lattice and a time-dependent part $V(\tau)=\dot{\mathbf{x}}_{0}(\tau)\cdot\mathbf{p}$ describing shaking, we evaluate their matrix elements in turn as follows.

The familiar $H_{\mathrm{static}}$ contains onsite and the nearest-neighbor hopping terms

[TABLE]

The onsite energy for each band is defined as

[TABLE]

and the hopping integrals for the $s$- and $p$-orbitals are given by

[TABLE]

where we have taken the lattice spacing to be one.

Similarly, we can split the time periodic shaking term $V(\tau)$ into onsite ($V_{0}$) and nearest neighbor ($V_{1}$) coupling terms,

[TABLE]

Here $\mathbf{a}(\tau)=\partial_{\tau}\mathbf{x}_{0}(\tau)$, i.e.,

[TABLE]

The symmetry of $V(\tau)$ dictates that for the onsite terms, only $s-p_{x,y}$ and $p_{x,y}-d_{xy}$ couplings are allowed,

[TABLE]

Such shaking induced band mixing is characterized by the coupling strength

[TABLE]

The nearest neighbor coupling terms have a similar matrix structure,

[TABLE]

Here shaking induces transitions between two orbitals on two neighboring sites with coupling strength

[TABLE]

Note that there are also diagonal terms given by

[TABLE]

The matrix elements of $V(\tau)$ obtained here Eq. (25-48) are crucial for our subsequent analysis.

To obtain a time-independent effective Hamiltonian, we use the RWA $\psi_{n}\rightarrow[U_{R}]_{nn^{\prime}}\psi_{i,n^{\prime}}$ for $\psi$ in Eq. (21) for a given site $i$. The transformation matrix is given by

[TABLE]

where the choice of the exponentials in the matrix $U_{R}$ is motivated by our numerical results in Fig. 1. Specifically, the band gap between $s$ and $p$-bands is of the same order of the band gap between $p$ and $d$-bands. Thus near resonance shaking couples the $s$ and $p$-bands, and also the $p$ and $d$-bands. $U_{R}$ accounts the interplay between these four orbitals by shifting $s$-band by energy $2\omega$ and $p$-bands by $\omega$ such that all four levels are within the same energy window. The second step of RWA is to drop remaining rapidly oscillating terms in $U^{\dagger}_{R}H^{\prime}_{0}U_{R}$. In particular, we find the diagonal terms $\mathrm{i}t_{s,p}^{\prime}$ in $V_{1}$ are removed by RWA. After Fourier transformation to momentum space, the resulting effective Hamiltonian for the shaken system takes a clean form,

[TABLE]

Here diagonal elements are the bare energies of four bands (see the dashed curves in Fig. 2)

[TABLE]

where $\varepsilon_{s}(k_{\mu})={e}_{s}+2{t}_{s}\cos(k_{\mu})$ and $\varepsilon_{p}(k_{\mu})={e}_{p}+2{t}_{p}\cos(k_{\mu})$, $\mu=x,y$. The off-diagonal terms are inter-orbital couplings induced by shaking,

[TABLE]

For given $V_{0}/E_{R}$, we can calculate parameters $e_{n}$, $t_{s,p}$, $d_{0,1}$ numerically from the Wannier functions constructed from the Bloch waves. For reference we provided a few typical values of these parameters in Table 1.

Eq. (54) is one of the main results of this paper. We can further rewrite $H_{\mathrm{eff}}$ as the sum of two direct products of the form

[TABLE]

with $\mathbb{1}_{\mu}$ the two-by-two unit matrix and $H_{\mu}$ given by

[TABLE]

where $\zeta_{\mu}=[\varepsilon_{s}(k_{\mu})+\omega+\varepsilon_{p}(k_{\mu})]/2$ and $h_{\mu}=[\varepsilon_{s}(k_{\mu})+\omega-\varepsilon_{p}(k_{\mu})]/2$. $H_{\mathrm{eff}}$ is diagonalized by a unitary transformation

[TABLE]

where

[TABLE]

The two angles here, $\theta_{\mu}$ and $\varphi_{\mu}$, are defined by

[TABLE]

Notice that $D_{y}$ couples to $h_{x}$ and $D_{x}$ couples to $h_{y}$ in the definitions of $\theta_{\mu}$ and $\varphi_{\mu}$. The four eigenvalues of $H_{\mathrm{eff}}$ are found as

[TABLE]

where $\kappa,\kappa^{\prime}=\pm$. For example, the hybridized $s$-band, modified from the bare $s$-band due to mixing with other bands by shaking, is given by

[TABLE]

Let us introduce the fermion creation operators in the basis of $H_{\mathrm{eff}}$, $\Phi^{\dagger}=(\phi^{\dagger}_{1},\phi^{\dagger}_{2},\phi^{\dagger}_{3},\phi^{\dagger}_{4})$, with $\phi^{\dagger}_{1}$ corresponding to the hybridized $s$-band $\epsilon_{1}(\mathbf{k})$ for example (the momentum and spin index are suppressed for brevity). They are related to the creation operators in the original basis $\Psi^{\dagger}\equiv(\psi_{s}^{\dagger},\psi_{p_{x}}^{\dagger},\psi_{p_{y}}^{\dagger},\psi_{d_{xy}}^{\dagger})$ by the unitary transformation $\Psi_{\mathbf{k}}=\mathcal{U}(\mathbf{k})\Phi_{\mathbf{k}}$ or $\psi_{n\mathbf{k}}=\mathcal{U}_{nm}(\mathbf{k})\phi_{m\mathbf{k}}$ where $\mathcal{U}(\mathbf{k})$ is given in Eq. (64). This relation will become important in the next section.

The band structure of $H_{\mathrm{eff}}$ is illustrated in Fig. 2 for lattice depth $V_{0}/E_{R}=5$. The top row compares the four energy bands (solid lines) described by $\epsilon_{\kappa\kappa^{\prime}}(\mathbf{k})$ in Eq. (70) with the bare band dispersions (dashed lines) given in Eq. (55). One can see that the level repulsion between the $s$-, $p$- and $d$-bands pushes the $s$-band up around the $\Gamma$ point. The dispersion of the hybridized $s$-band in the entire 2D Brillouin zone and the corresponding density of states are shown in the middle row of Fig. 2. The four band minima are clearly seen here. The spectrum obtained here is in good agreement with the numerical solution in section III. When the hybridized $s$-band is gradually filled with fermions, the resulting Fermi surfaces undergoes a non-trivial evolution as shown in the bottom row of Fig. 2. The Fermi surface topology change found here highlights the capability of lattice shaking in engineering the band structures.

The analytical form of the effective Hamiltonian $H_{\mathrm{eff}}$ and the resulting spectrum clarify the physics of shaken square optical lattice. It captures succinctly how the relevant orbitals, either on the same site or two neighboring sites, are coupled by shaking. The simplicity achieved is partly due to our choice of a convenient gauge, where $H^{\prime}_{0}$ splits into $H_{\mathrm{static}}$ and $V(\tau)=\dot{\mathbf{x}}_{0}(\tau)\cdot\mathbf{p}$. It is also derived from the symmetries of the Wannier functions and $V(\tau)$, leading to for example only two independent coupling strength $d_{0}$ and $d_{1}$. These results from the treatment of circularly shaken square lattice may be useful for the study of other lattice geometries and shaking protocols.

V Effective Interactions

The peculiar Fermi surfaces found in the previous section suggests that interaction may drive interesting many-body instabilities for Fermi gases in shaken lattices. The effective interactions for two fermions on the Fermi surface of the hybridized bands will differ from the bare interactions. In this section, we outline a procedure to derive these effective interactions $V_{\mathrm{eff}}$ and then work out its explicit expression for fermions on the hybridized $s$-band. We only consider weak interactions of spin-$1/2$ fermions in the sense that the energy scale related to interactions is assumed to be much smaller than the shaking frequency and the band-width. In other words, we will treat interaction as a weak perturbation to $H_{\mathrm{eff}}$ in Eq. (54).

We first expand the contact interaction in Eq. (8) in the Wannier basis using Eq. (20). It then takes the Hubbard-like onsite form in the multi-orbital basis,

[TABLE]

where the time dependence of $\psi$ is suppressed for brevity and $U_{nmm^{\prime}n^{\prime}}\equiv g\int d\mathbf{x}\mathcal{W}_{n}^{*}(\mathbf{x})\mathcal{W}_{m}^{*}(\mathbf{x})\mathcal{W}_{m^{\prime}}(\mathbf{x})\mathcal{W}_{n^{\prime}}(\mathbf{x})$. In accordance with the previous section, we only keep the orbitals $n=s,p_{x},p_{y},d_{xy}$. Also we will approximate the Wannier functions with local harmonic oscillator eigenstates to evaluate the integrals $U_{nmm^{\prime}n^{\prime}}$. This is justified for deep lattices and it simplifies the algebra greatly. In fact, as we will show below, all the onsite interactions in different orbital channels can be expressed in terms of the $s$-orbital interaction constant

[TABLE]

Next we address the question what happens to the interaction term $V_{I}$ during the RWA. To answer this question in a transparent way, we split the terms in $V_{I}$ into three distinct channels

[TABLE]

The first term is the density-density interaction given by

[TABLE]

where $\vec{\mathrm{n}}_{\sigma}^{\top}=[\mathrm{n}_{\sigma s},\mathrm{n}_{\sigma x},n_{\sigma y},\mathrm{n}_{\sigma d}]$, the superscript $\top$ denotes matrix transposition and the density operator is defined as $\mathrm{n}_{\sigma n}=\psi_{\sigma n}^{\dagger}\psi_{\sigma n}$ with the site index $i$ dropped for brevity. The elements of the constant matrix

[TABLE]

are obtained by evaluating the overlap integrals using the approximate Wannier functions. Now consider the effect of RWA, $[\psi_{\sigma s},\psi_{\sigma x},\psi_{\sigma y},\psi_{\sigma d}]\overset{\text{\tiny RWA}}{\rightarrow}[\mathrm{e}^{-2\mathrm{i}\omega\tau}\psi_{\sigma s},\mathrm{e}^{-\mathrm{i}\omega\tau}\psi_{\sigma x},\mathrm{e}^{-\mathrm{i}\omega\tau}\psi_{\sigma y},\psi_{\sigma d}]$, on these terms. One can see that for every $\psi_{\sigma n}$ operator, there is a corresponding $\psi_{\sigma n}^{\dagger}$. Therefore no time dependent terms like $\mathrm{e}^{-\mathrm{i}\omega\tau}$ and $\mathrm{e}^{-2\mathrm{i}\omega\tau}$ will remain, thus $V_{\mathrm{density}}$ is invariant under RWA.

The second term is the orbital exchange interaction

[TABLE]

where

[TABLE]

with $\vec{\mathrm{S}}^{+}=[\mathrm{S}^{+}_{s},\mathrm{S}^{+}_{x},\mathrm{S}^{+}_{y},\mathrm{S}^{+}_{d}]$, and $\vec{\mathrm{S}}^{-}=[\mathrm{S}^{-}_{s},\mathrm{S}^{-}_{x},\mathrm{S}^{-}_{y},\mathrm{S}^{-}_{d}]^{\top}$. The raising and lowering operators are defined as usual $\mathrm{S}_{n}^{\pm}=\mathrm{S}_{n}^{1}\pm\mathrm{i}\mathrm{S}_{n}^{2}$, with $\mathrm{S}_{n}^{\mu}=\frac{1}{2}\gamma^{\mu}_{\sigma\sigma^{\prime}}\psi_{\sigma n}^{\dagger}\psi_{\sigma^{\prime}n}$ and $\boldsymbol{\gamma}=(\gamma^{1},\gamma^{2},\gamma^{3})$ are the Pauli matrices. Similar to density-density interactions, one can see that there is a $\psi_{\sigma n}^{\dagger}$ for every $\psi_{\sigma^{\prime}n}$ operator. Thus, $V_{\mathrm{ex}}$ is also invariant under RWA.

Finally, the last term in Eq. (73) describes pair transfers between the orbitals,

[TABLE]

where $\vec{\mathrm{P}}^{\dagger}=[\mathrm{P}_{s}^{\dagger},\mathrm{P}_{x}^{\dagger},\mathrm{P}_{y}^{\dagger},\mathrm{P}_{d}^{\dagger}]$, $\mathrm{P}_{n}^{\dagger}=\psi_{\uparrow n}^{\dagger}\psi_{\downarrow n}^{\dagger}$ is the pair creation operator for orbital $n$, and $\hat{G}_{\mathrm{pt}}$ is identical to $\hat{G}_{\mathrm{ex}}$ above. The pair transfer between the $p_{x}$ and $p_{y}$ orbitals is invariant under RWA since they have the same exponential time dependence. However, for pair transfers between $s$ and $p$, $s$ and $d$ as well as $p$ and $d$ orbitals, time dependent exponentials will remain after the RWA. Since such fast oscillating terms are subsequently ignored in RWA, $\hat{G}_{\mathrm{pt}}$ becomes simplified,

[TABLE]

It follows that after the RWA, the (time-independent) effective interaction in the momentum space takes the following form,

[TABLE]

Here $\psi$ no longer depends on $\tau$, and $U^{\prime}_{nmm^{\prime}n^{\prime}}$ differs from $U_{nmm^{\prime}n^{\prime}}$ by absence of all pair transfer terms but the one in between $p_{x}$ and $p_{y}$ orbitals, in accordance with Eq. (79). $U^{\prime}$ can be straightforwardly constructed from the $\hat{G}$ matrices above [Eqs. (75), (77) and (79)] and its various terms will not be tabulated here.

The last step is to rewrite $V^{\prime}_{I}$ in terms of the field operators in the basis of $H_{\mathrm{eff}}$. This is achieved by the unitary transformation, $\psi_{\sigma n\mathbf{k}}=\mathcal{U}_{nm}(\mathbf{k})\phi_{\sigma m\mathbf{k}}$, with $\mathcal{U}(\mathbf{k})$ given in Eq. (64). As an example, let us focus on the effective interactions for two fermions of opposite momenta on the hybridized $s$-band, denoted with $V^{s}_{I}$ below. For this purpose, we can project $V^{\prime}_{I}$ onto the $n=1$ band by substituting $\psi_{\sigma m\mathbf{k}}=\mathcal{U}_{m1}\phi_{\sigma 1\mathbf{k}}$ and $\psi^{\dagger}_{\sigma m\mathbf{k}}=\phi^{\dagger}_{\sigma 1\mathbf{k}}\mathcal{U}^{\dagger}_{1m\mathbf{k}}$ into the expression for $V^{\prime}_{I}$ and collecting the relevant terms to get

[TABLE]

where we have factored out onsite $s$-band interaction constant $U$ for convenience and defined the interaction vertex $V_{\mathrm{eff}}(\mathbf{k},\mathbf{k}^{\prime})$ as

[TABLE]

This expression can be simplified using the matrix elements of $\mathcal{U}(\mathbf{k})$ in Eq. (64) with some algebra as

[TABLE]

The first two terms are given by

[TABLE]

in terms of

[TABLE]

where $\theta_{\mu}$ is defined in Eq. (68). And the pair transfer term takes the following form,

[TABLE]

Note that we have considered the interactions between two particles with zero center of mass momentum above. It is straightforward to obtain more general interaction vertices of the form

[TABLE]

by a similar projection procedure. The result is rather lengthy and will not be given here.

The most interesting property of $V_{\mathrm{eff}}(\mathbf{k},\mathbf{k}^{\prime})$ is its nontrivial $\mathbf{k}$ dependence. This is in marked contrast with the bare interaction which is constant in $\mathbf{k}$ space. In analogy with the hybridized band dispersion discussed in the previous section, we may say that the interaction is strongly modified by shaking induced band mixing. Fig. 3 shows a few examples of $V_{\mathrm{eff}}(\mathbf{k},\mathbf{k}^{\prime})$ for different fillings, where $\mathbf{k}$ and $\mathbf{k}^{\prime}$ reside on the corresponding Fermi surfaces. For $n=0.55$, the variation in $V_{\mathrm{eff}}(\mathbf{k},\mathbf{k}^{\prime})$ reaches the order of 50%. Formally, the momentum dependence of $V_{\mathrm{eff}}(\mathbf{k},\mathbf{k}^{\prime})$ is inherited from the $\mathbf{k}$-dependence of the $\mathcal{U}$ matrix under the projection procedure. We can understand the variation of $V_{\mathrm{eff}}(\mathbf{k},\mathbf{k}^{\prime})$ qualitatively as follows. The $p$-orbitals hybridize with the $s$-orbital in a non-uniform way around the Fermi surface. In particular, the mixing is stronger near the $\Gamma\mathrm{X}$ line, where the ridge of the $p$-band come close to the $s$-band (see Fig. 2), than that along the $\Gamma\mathrm{M}$ line. Since the onsite interaction constant of the bare $p$-bands is smaller than the $s$-band due to the reduction in the overlap integrals, regions on Fermi surface with more mixture of $p$-orbitals have smaller effective interaction. Thus, we expect the effective interaction reaches maximum around, e.g., $\varphi_{\mathbf{k}}\sim\varphi_{\mathbf{k}^{\prime}}\sim\pi/4$ (where $\varphi_{\mathbf{k}}$ is the polar angle on the Fermi surface) along the $\Gamma\mathrm{M}$ line, consistent with the numerical results in Fig. 3.

VI Pairing of Fermions in shaken square lattice

To summarize, we have arrived at the following effective Hamiltonian for weakly interacting fermions in the hybridized $s$-band of the shaken lattice,

[TABLE]

Here we have dropped the orbital index $n=1$, $\phi_{\sigma 1\mathbf{k}}\rightarrow\phi_{\sigma\mathbf{k}}$, and defined $\xi_{\mathbf{k}}=\epsilon_{1}(\mathbf{k})-\mu$ with chemical potential $\mu$. The band dispersion $\epsilon_{1}(\mathbf{k})$ is given in Eq. (70) and the reduced effective interaction $V_{\mathrm{eff}}(\mathbf{k},\mathbf{k}^{\prime})$ is given in Eq. (83). We only consider weak, attractive interactions with $U<0$. The dominant instability is Cooper pairing, and it is justified to only include $V_{\mathrm{eff}}(\mathbf{k},\mathbf{k}^{\prime})$, i.e., the scattering of fermions with opposite spins and opposite momenta. For strong coupling or repulsive interactions, one needs to consider more general interaction vertices $V(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3})$ and other relevant instabilities.

Following the standard BCS theory, we define the pairing order parameter

[TABLE]

and perform mean-field decoupling of the quartic interaction terms in Eq. (89). The resulting quadratic Hamiltonian is diagonalized using Bogoliubov transformation, and the gap equation becomes

[TABLE]

where $E_{\mathbf{k}}=\sqrt{\xi_{\mathbf{k}}^{2}+|\Delta(\mathbf{k})|^{2}}$ is the quasiparticle excitation spectrum. For temperatures close to the critical temperature $T_{c}$, the magnitude of the gap is small. The gap equation can be linearized to become an eigenvalue problem

[TABLE]

Here $\varphi_{\mathbf{k}^{\prime}}$ is the angle of the momenta $\mathbf{k}^{\prime}$ with respect to the center of the Fermi surface. The largest eigenvalue $\lambda$ yields $T_{c}$ through the non-linear equation

[TABLE]

And the corresponding eigenvector gives the orbital symmetry of the pairing order parameter. Here $N_{\mathrm{sites}}$ is the number of lattice sites. For cases with multiple Fermi surfaces (see bottom row of Fig. 2 for $n=0.15$ for example), Eq. (92) is solved separately for each Fermi surface such that the Fermi surface with the largest eigenvalue determines the leading unstable surface.

Due to the momentum dependence of $V_{\mathrm{eff}}(\mathbf{k},\mathbf{k}^{\prime})$, the solution to Eq. (92) will yield an order parameter $\Delta(\mathbf{k})$ that depends on the angular location $\varphi$ of $\mathbf{k}$ on the Fermi surface. In other words, $\Delta(\mathbf{k})$ is in general anisotropic and includes higher harmonics,

[TABLE]

Fig. 4 shows the phase diagram of $H_{\mathrm{BCS}}$ as functions of filling and temperature. In the gray shaded region, we find the pairing is predominantly of $s+d$-wave symmetry, i.e., with a small $\Delta_{d}$ component and higher harmonics can be neglected. This is illustrated for three different fillings marked as $A$, $B$ and $C$ in the top row of Fig. 4. The corresponding Fermi surfaces are shown in the middle row, and the order parameters $\Delta(\varphi)$ are shown in the bottom row. The reduction in pairing amplitude along the $\Gamma\mathrm{X}$ line is consistent with the weaker effective interaction there found earlier in section V. Thus the $s+d$-wave pairing is a direct result of shaken induced anisotropy of the effective interactions. For other fillings outside the shaded region, the pairing symmetry is the usual $s$-wave.

VII Red-detuned Near-resonance Shaking

So far we have focused on shaking frequencies near resonance, but blue detuned from the band separation (not the band gap) of the $s$- and $p$-bands at the Brillouin zone center, $\Omega=\max E_{p}-\min E_{s}$. For example, we have set $\omega/\Omega=1.01$ and $1.05$. Now we move on to discuss red detuned shaking frequencies $\omega<\Omega$, e.g., $\omega/\Omega=0.85$, and compare them to the blue detuned case. For red detuned frequencies, the folded $s$- and $p$-bands directly cross each other (see the dashed lines in the top panel of Fig. 5), and hybridize strongly near these crossing points located away from the Brillouin zone center. We emphasize that the theory developed in the previous sections are valid for all near resonance shaking frequencies, $\omega\sim\Omega$. The expression for $H_{\mathrm{eff}}$, the band dispersion in Eq. (70), and the effective interactions in Eq. (83) can be directly applied to red detuned cases without any change.

The calculation for $\omega/\Omega=0.85$ proceeds the same way as before, and the results are summarized in Fig. 5. A main difference from the blue detuned case is the orbital character of the hybridized $s$-band, shown at the top row of Fig. 2 with solid line. It is predominantly of $d_{xy}$-orbital character at the Brillouin zone center $\Gamma$, becomes more $p$-orbital like on the side of the Brillouin zone around $\mathrm{X}$, and remains $s$-orbital like at the Brillouin zone corner $\mathrm{M}$. The shaking induced band hybridization is strongest near the crossing points of the folded bare bands (dashed lines). As shown in the middle panel of Fig. 5, the dispersion of the hybridized $s$-band also has four minima. Compared to the blue detuned case, the location of the energy minima, characterized by momentum vector $\mathbf{Q}$, is further away from $\Gamma$ even though a smaller shaking amplitude $\beta=0.05$ is used. The shapes of the Fermi surfaces are quite different as indicated by the contour lines of $\epsilon({k_{x},k_{y}})$. For example, at filling $n=0.45$, the Fermi surface is centered at the Brillouin zone corner $M$.

In the red detuned regime, the anisotropy of effective interactions on the Fermi surface is much more pronounced. Take again $n=0.45$ as an example. The effective interaction in the pairing channel $V(\varphi,\varphi^{\prime})$ varies by as much as seven folds. This is because the hybridized $s$-band contains a significant contribution from the $d$-band, the bare interaction of which is smaller than those of the $s$- and $p$-bands. As a result, $V$ is smaller for $\mathbf{k}$ points closer to $\Gamma$, e.g. the red dots in Fig. 5, and thus of more $d$-orbital character. Due to the anisotropy of $V$, the pairing order parameter has $s+d$-wave symmetry, as shown in the bottom panel of Fig. 5. It differs slightly from the one given in Fig. 4 in two aspects. First, the $d$-wave component of the order parameter is increased. Secondly, $\Delta(\varphi)$ reaches minimum at $\varphi=\pi/4$ where the effective interaction is weak (recall again that in this case, the center of Fermi surface is at $M$).

VIII Concluding remarks

Lattice shaking provides a valuable tool to engineer the band structures and effective interactions for cold atoms in optical lattice beyond the reach of the static optical lattices. The key physics at play is the mixing or hybridization of different orbital bands induced by lattice shaking. Presently, the new phases arising from interactions in these mixed bands remain largely unexplored. Our work constitutes a first step towards a quantitative understanding of the interaction effects for fermionic atoms in shaken square optical lattice. We derived a four-band effective Hamiltonian to clarify the matrix elements for inter-orbital couplings and obtained the analytical expressions for the eigenenergy bands, in good agreement with the full numerical Floquet analysis. We further derived the effective interactions on the hybridized $s$-band and explained the origin of its acquired momentum dependence. Applying the theory to spin 1/2 fermions with attractive interactions, we monitored the nontrivial evolution of the Fermi surface and found the symmetry of the pairing order parameter can be $s+d$ wave. The similarities and differences between the red and blue detuned shaking frequency are discussed using insights gained from the analytical understanding of $H_{\mathrm{eff}}$ and $V_{\mathrm{eff}}(\mathbf{k},\mathbf{k}^{\prime})$. These concrete examples support the presence of complex Fermi surfaces, anisotropic interactions, and interesting many-body phases for fermions in shaken lattices.

Our work can be generalized in several directions. For example, for repulsive interactions, the general effective interaction vertex can be obtained and subsequently applied to discuss competing many-body phases including spin density waves, superfluidity, and Pomeranchuk instability etc. In particular, interesting phenomena are expected when the Fermi surface are partially nested, i.e., with large segments connected by some common nesting wave vector. The technical procedures outlined here to derive the effective Hamiltonian will also be useful to treat other regimes of shaking frequencies, such as two-photon resonances known to give rise topologically nontrivial bands. Lastly, another interesting generalization is to analyze other lattice geometries, e.g., interacting fermions on the shaken honeycomb and checkerboard lattice.

Acknowledgements.

We thank C. Chin for helpful discussions. This work is supported by AFOSR Grant No. FA9550-16-1-0006 (A.K., E.Z., and W.V.L.), ARO Grant No. W911NF-11-1-0230 (A.K. and W.V.L.), NSF Grant No. PHY-1205504 (A.K. and E.Z.), and NSF of China Overseas Scholar Collaborative Program Grant No. 11429402 sponsored by Peking University (W.V.L.).