Stabilization of prethermal Floquet steady states in a periodically driven dissipative Bose-Hubbard model

TL;DR

This paper investigates how dissipation can stabilize prethermal Floquet steady states in a driven Bose-Hubbard model, preventing heating and enabling stable energy absorption from external drives.

Contribution

It demonstrates that dissipation suppresses parametric instabilities, leading to stable steady states in a periodically driven Bose-Hubbard system, a novel insight into controlling heating in quantum many-body systems.

Findings

Dissipation suppresses parametric instabilities.

Stable steady states form when dissipation exceeds on-site interactions.

Driven systems can emit energy to reservoirs, avoiding heating.

Abstract

We discuss the effect of dissipation on heating which occurs in periodically driven quantum many body systems. We especially focus on a periodically driven Bose-Hubbard model coupled to an energy and particle reservoir. Without dissipation, this model is known to undergo parametric instabilities which can be considered as an initial stage of heating. By taking the weak on-site interaction limit as well as the weak system-reservoir coupling limit, we find that parametric instabilities are suppressed if the dissipation is stronger than the on-site interaction strength and stable steady states appear. Our results demonstrate that periodically-driven systems can emit energy, which is absorbed from external drivings, to the reservoir so that they can avoid heating.

Equations

$$H_S(t)=$$

$$+\sum _j(\mathbf{E}(t)\cdot {\mathbf{r}}_j{n}_j-\mu n_j),$$

$$H_R=\sum _k{\omega }_k{A}_k^{†}{A}_k$$

$$\displaystyle H_{I}=\sum_{k,j}\bigg{[}\frac{\lambda_{k,j}}{\sqrt{L}}A^{\dagger}_{k}B_{j}+H.c.\bigg{]},$$

$$A_k={\Psi }_k+a_k$$

$$B_j={\Phi }_j+b_j,$$

$$i\frac{d}{dt}{A}_k(t)={\omega }_k{A}_k(t)+\sum _j\frac{{\lambda }_{k,j}}{\sqrt{L}}{B}_j(t)$$

$$\begin{array}{rl}i\frac{d}{dt}{B}_j(t)=& -t_0\sum _{⟨i,j⟩}{B}_i(t)+(\mathbf{E}(t)\cdot {\mathbf{r}}_j-\mu )B_j(t)\\ & +U{B}_j^{†}(t)B_j(t)B_j(t)+\sum _k\frac{{\lambda }_{k,j}^{\ast }}{\sqrt{L}}{A}_k(t).\end{array}$$

$$B_j^{†}{B}_j{B}_j\simeq$$

$$+2(|{\Phi }_j{|}^2+⟨b_j^{†}{b}_j⟩)b_j+({\Phi }_j^2+⟨b_j{b}_j⟩)b_j^{†}.$$

$$A_k(t)=A_k(0)e^{-i{\omega }_{k}t}-i\sum _j{\int }_0^{t}d{t}^{\prime }\frac{{\lambda }_{k,j}}{\sqrt{L}}{e}^{-i{\omega }_k{t}^{\prime }}{B}_j(t-t^{\prime }).$$

$$\displaystyle\begin{aligned} i\frac{d}{dt}\Phi_{j}(t)=&-t_{0}\sum_{\braket{i,j}}e^{i\bm{\mathcal{A}}(t)\cdot(\bm{r}_{j}-\bm{r}_{i})}\Phi_{i}(t)-\mu\Phi_{j}(t)\\ &+U\big{[}(n_{j}(t)+2\tilde{n}_{j}(t))\Phi_{j}(t)+\tilde{m}_{j}(t)\Phi^{\ast}_{j}(t)\big{]}\\ &-i\sum_{j^{\prime}}\int_{0}^{t}dt^{\prime}\Gamma_{jj^{\prime}}(t,t^{\prime})\Phi_{j^{\prime}}(t-t^{\prime})+i\Xi_{j}(t)\end{aligned}$$

$$\begin{array}{rl}i\frac{d}{dt}{b}_j(t)=& -t_0\sum _{⟨i,j⟩}{e}^{i\mathcal{A}(t)\cdot ({\mathbf{r}}_j-{\mathbf{r}}_i)}{b}_i(t)-\mu b_j(t)\\ & +2U(n_j(t)+{\tilde{n}}_j(t))b_j(t)\\ & +U(m_j(t)+{\tilde{m}}_j(t))b_j^{†}(t)\\ & -i\sum _{j^{\prime }}{\int }_0^{t}d{t}^{\prime }{\Gamma }_{j{j}^{\prime }}(t,t^{\prime })b_{j^{\prime }}(t-t^{\prime })+i{\xi }_j(t),\end{array}$$

$$n_j(t)$$

$$m_j(t)$$

$${\tilde{n}}_j(t)$$

$${\tilde{m}}_j(t)$$

$${\Phi }_{\mathbf{p}}(t)$$

$$b_{\mathbf{p}}(t)$$

$$E_{\mathbf{p}}$$

$${\theta }_{\mathbf{p}}(t)$$

$${\Phi }_{\mathbf{p}}^{(I)}(t)=e^{i{E}_{\mathbf{p}}t+i{\theta }_{\mathbf{p}}(t)}{\Phi }_{\mathbf{p}}(t)$$

$$b_{\mathbf{p}}^{(I)}(t)=e^{i{E}_{\mathbf{p}}t+i{\theta }_{\mathbf{p}}(t)}{b}_{\mathbf{p}}(t).$$

$$i\frac{d}{dt}{\Phi }_{\mathbf{q}}^{(I)}(t)=u(t){\Phi }_{\mathbf{q}}^{(I)}(t)+v(t){\Phi }_{\mathbf{q}}^{(I)\ast }(t)-i{\int }_0^{t}d{t}^{\prime }{\Gamma }_{\mathbf{q}\mathbf{q}}(t^{\prime }){\Phi }_{\mathbf{q}}^{(I)}(t-t^{\prime })+i{\Xi }_{\mathbf{q}}^{(I)}(t)$$

$$i\frac{d}{dt}{b}_{\mathbf{p}}^{(I)}(t)=\tilde{u}(t)b_{\mathbf{p}}^{(I)}(t)+{\tilde{v}}_{\mathbf{p}}(t)b_{2\mathbf{q}-\mathbf{p}}^{(I)†}(t)-i\sum _{{\mathbf{p}}^{\prime }\in \mathcal{S}(E_{\mathbf{p}})}{\int }_0^{t}d{t}^{\prime }{\Gamma }_{\mathbf{p}{\mathbf{p}}^{\prime }}(t^{\prime })b_{{\mathbf{p}}^{\prime }}^{(I)}(t-t^{\prime })+i{\xi }_{\mathbf{p}}^{(I)}(t),$$

$$\displaystyle v(t)=U\big{[}\tilde{m}_{\bm{q}}(t)+\sum_{\bm{p}\in\mathcal{R}}(h_{\bm{p},\bm{q}}^{(1)})^{\ast}\tilde{m}_{\bm{p}}(t)e^{-2i(E_{\bm{p}}-E_{\bm{q}}-\Omega/2)t}\big{]}$$

$${\tilde{v}}_{\mathbf{p}}(t)$$

$$n_0(t)$$

$$m_0(t)$$

$$\tilde{n}(t)$$

$${\tilde{m}}_{\mathbf{p}}(t)$$

$${\Phi }_{\mathbf{q}}^{(I)}(t)$$

$$b_{\mathbf{q}}^{(I)}(t)$$

$$b_{\mathbf{p}\ne \mathbf{q}}^{(I)}(t)$$

$$\begin{array}{rl}i\frac{d}{dt}{\Phi }_{\mathbf{q}}^{(I)}(t)& =g{\Phi }_{\mathbf{q}}^{(I)\ast }(t)\\ & -i{\int }_0^{t}d{t}^{\prime }{\Gamma }_{\mathbf{q}\mathbf{q}}(t^{\prime }){\Phi }_{\mathbf{q}}^{(I)}(t-t^{\prime })+i{\Xi }_{\mathbf{q}}^{(I)}(t)\end{array}$$

$$\displaystyle\begin{aligned} i\frac{d}{dt}b^{(I)}_{\bm{q}}(t)&=n_{0}Ub^{(I)}_{\bm{q}}(t)+\big{[}m_{0}U+g\big{]}b^{(I)\dagger}_{\bm{q}}(t)\\ &-i\int_{0}^{t}dt^{\prime}\Gamma_{\bm{q}\bm{q}}(t^{\prime})b^{(I)}_{\bm{q}}(t-t^{\prime})+i\xi^{(I)}_{\bm{q}}(t)\end{aligned}$$

$$\begin{array}{rl}i\frac{d}{dt}{b}_{\mathbf{p}\ne \mathbf{q}}^{(I)}(t)& ={\tilde{g}}_{\mathbf{p}}{b}_{2\mathbf{q}-\mathbf{p}}^{(I)†}(t)+i{\xi }_{\mathbf{p}}^{(I)}(t)\\ & -i\sum _{{\mathbf{p}}^{\prime }\in \mathcal{S}(E_{\mathbf{p}})}{\int }_0^{t}d{t}^{\prime }{\Gamma }_{\mathbf{p}{\mathbf{p}}^{\prime }}(t^{\prime })b_{{\mathbf{p}}^{\prime }}^{(I)}(t-t^{\prime }).\end{array}$$

$$\sum_{\bm{p}^{\prime}\in\mathcal{R}}\big{(}\delta_{\bm{p},\bm{p}^{\prime}}z+\Gamma_{\bm{p}\bm{p}^{\prime}}(z)\big{)}b^{(I)}_{\bm{p}^{\prime}}(z)+i\tilde{g}_{\bm{p}}b^{(I)\dagger}_{2\bm{q}-\bm{p}}(z)=b_{\bm{p}}+\xi^{(I)}_{\bm{p}}(z).$$

$$M(z){\mathbf{b}}^{(I)}(z)=\mathbf{b}+\mathbf{\xi }(z),$$

$${\mathbf{b}}^{(I)}(t)$$

$$\displaystyle+\sum_{k}\sum_{n\in\mathbb{Z}}\mathop{\rm Res}\limits_{z=\pm i(\omega_{k}-\epsilon_{\bm{p}}^{(n)})}\big{[}M^{-1}(z)\bm{\xi}(z)e^{zt}\big{]},$$

Full text

I Introduction

Long time behavior of periodically driven quantum systems has attracted particular attention due to the success of the Floquet engineering of mattersLignier et al. (2007); Struck et al. (2012); Aidelsburger et al. (2013, 2011); Miyake et al. (2013); Aidelsburger et al. (2015); Wang et al. (2013); Ishikawa et al. (2014); Sie et al. (2015). Dynamics of such system is described by the Floquet effective Hamiltonian which is defined by $U(T)=\exp[-iH_{{\rm eff}}T]$, where $U(t)$ is the time evolution operator of the system and $T$ is the period of the external field. This Floquet effective Hamiltonian is usually obtained by truncating the Magnus expansion and the truncated Floquet effective Hamiltonian can have non-trivial properties such as topologically non-trivial band structures Oka and Aoki (2009); Lindner, Netanel H. and Refael, Gil and Galitski, Victor (2011); Kitagawa et al. (2011) which the original Hamiltonian does not have.

Recent studies have carefully investigated essential features of the Floquet effective Hamiltonian of quantum isolated non-integrable systems and revealed that the system absorbs energy from the external field and steady states realized after a long time are the trivial infinite temperature stateD’Alessio and Rigol (2014); Lazarides et al. (2014); Ponte et al. (2015). Further studies have also shown that the energy absorption rate is small and there exists quasi-steady states (i.e. their life time is finite) which have finite temperature properties of the truncated Floquet effective HamiltonianAbanin et al. (2015, 2017); Mori et al. (2016); Kuwahara et al. (2016); Canovi et al. (2016); Lindner et al. (2016); Weidinger and Knap (2016); Chandran and Sondhi (2016). These facts suggest that there remains a possibility that the life time of quasi-steady states in isolated systems becomes infinite in dissipative systems because the system can emit the energy to the reservoir. Solid state systems irradiated by intense laser lights are essentially a quantum dissipative system and in this sense, understanding the effect of dissipation on the heating in periodically driven quantum systems is a fundamentally important issue.

There has been also intensive research on periodically driven dissipative quantum systemsZerbe and Hänggi (1995); Kohler et al. (1997); Breuer et al. (2000); Kohn (2001); Hone et al. (2009); Ketzmerick and Wustmann (2010); Iadecola et al. (2013); Langemeyer and Holthaus (2014); Iadecola and Chamon (2015); Iadecola et al. (2015); Shirai et al. (2015); Dehghani et al. (2014); Liu (2015); Seetharam et al. (2015); Shirai et al. (2016); Iwahori and Kawakami (2016). Recent studies have revealed that long time asymptotic states of periodically driven dissipative quantum systems are effectively described by the Gibbs distribution of the Floquet effective Hamiltonian in the high frequency regime of the external fieldSeetharam et al. (2015); Shirai et al. (2016); Iwahori and Kawakami (2016). However, these studies are restricted to integrable systems or numerical calculations for a small non-integrable spin chain. Moreover, most of previous studies concentrated on the case where the system reservoir coupling is taken to be infinitesimal compared with any other energy scales of the system, therefore it is not clear how the heating which is caused by periodic drivings is affected by the dissipation.

Motivated by these backgrounds, we investigate the effect of dissipation on the heating in the periodically driven quantum systems. We consider a periodically driven Bose-Hubbard model coupled to a reservoir and the on-site interaction is treated by the Hartree-Fock-Bogoliubov approximation. If there is no dissipation, this system is known to exhibit parametric instabilities with an infinitesimal on-site interactionCreffield (2009); Bukov et al. (2015); Lellouch et al. (2016) and recent studies have revealed that the heating before the quasi steady state is governed by the parametric instabilitiesChandran and Sondhi (2016); Weidinger and Knap (2016). By taking the weak on-site interaction limit as well as the weak system-reservoir coupling limit, we find that the Floquet steady state is stabilized if the dissipation is stronger than the on-site interaction strength. Our results demonstrate that heating which is inevitable in the periodically driven isolated quantum systems can be suppressed by introducing dissipation and we can obtain the low energy properties of the Floquet effective Hamiltonian in the prethermal Floquet steady state if the dissipation is strong enough to overcome the heating rate.

In the next section, we show the model considered in this paper and derive the equation of motion of the system within the Hartree-Fock-Bogoliubov approximation. In Sec. III, the equation of motion in the weak system-reservoir coupling and the weak on-site interaction limit is described, then we discuss the stability of the system and properties of the steady state in Sec. IV. A summary of our results is presented in Sec. VI.

II Setup

We study a periodically driven Bose-Hubbard model on a $d$-dimensional cubic lattice which couples to a reservoir. The Hamiltonian of the total system is given by $H_{\text{tot}}(t)=H_{S}(t)+H_{R}+H_{I}$, where $H_{R}$ and $H_{S}(t)$ are the Hamiltonian of the reservoir and the system, respectively, and $H_{I}$ describes the interaction between them. The Hamiltonian of the system is the Bose-Hubbard model driven by a time periodic force $\bm{E}(t)$:

[TABLE]

where $B_{j}$ is the boson annihilation operator for particle on site $\bm{r}_{j}$, $n_{j}=B^{\dagger}_{j}B_{j}$ is the number operator, and $\braket{i,j}$ represents a pair of nearest neighbor sites. The number of lattice site is taken as $N^{d}$ and $\mu$ is the chemical potential. The external force has a period of $2\pi/\Omega$ and the time average of the external force is zero : $\bm{E}(t+2\pi/\Omega)=\bm{E}(t),\ \int_{0}^{2\pi/\Omega}dt\bm{E}(t)=\bm{0}$. The Hamiltonian of reservoir and the interaction between the system and the reservoir are taken as

[TABLE]

where $A_{k}$ describes the boson in the reservoir and $\lambda_{k,j}$ determines the coupling between the system and reservoir. We consider the continuum limit $\sum_{k}\frac{2\pi}{L}F(\omega_{k})=\int_{-\infty}^{\infty}d\omega D(\omega)F(\omega)$, where $D(\epsilon)$ is the density of states of the reservoir and $F(\epsilon)$ is an arbitrary function. The system and the reservoir exchange particles and energy because $H_{I}$ is a bilinear form of the Bose operators.

We consider the BEC state and separate out the condensate part as

[TABLE]

where $\Psi_{k}=\braket{A_{k}}$ and $\Phi_{j}=\braket{B_{j}}$ describe macroscopic Bose fields. The initial state of the reservoir is taken as the equilibrium state with the temperature $T$, that is, $\braket{a^{\dagger}_{k}a_{k}}=f_{T}(\omega_{k})$, where $f_{T}(\omega)=1/(e^{\omega/T}-1)$ is the Bose distribution function (setting $k_{B}=1$ and $\hbar=1$). We take the chemical potential of the reservoir as zero and thereby the energy of the reservoir is positive $\omega_{k}\geq 0$. The initial state of the macroscopic Bose field in the reservoir is taken as $\Psi_{k}=\Psi\delta_{k,k_{0}}$, where $k_{0}$ specifies the smallest energy mode of the reservoir ($\text{min}_{k}\omega_{k}=\omega_{k_{0}}$ ).

The Heisenberg equation of the total system is written as

[TABLE]

We treat the on-site interaction term by the Hartree-Fock-Bogoliubov approximationGriffin (1996):

[TABLE]

The Hartree-Fock-Bogoliubov approximation cannot describe the scattering of quasi-particles correctly. However, within the weak on-site interaction limit which we will consider later, this process is irrelevantBukov et al. (2015); Lellouch et al. (2016). Eq. (4) is integrated as

[TABLE]

Substituting Eq. (7) into Eq. (5), we can obtain the closed form of the equation of motion for the system. We here also apply a gauge transformation $B_{j}(t)\rightarrow B_{j}(t)\exp[-i\bm{\mathcal{A}}(t)\cdot\bm{r}_{j}]$, where $\bm{\mathcal{A}}(t)=\int_{0}^{t}dt^{\prime}\bm{E}(t^{\prime})$. Then, we have the following equation of motion:

[TABLE]

where $\Gamma_{jj^{\prime}}(t,t^{\prime})=\sum_{k}\lambda^{\ast}_{k,j}(t)\lambda_{k,j^{\prime}}(t-t^{\prime})e^{-i\omega_{k}t}/L$ and $\Xi_{j}(t)=-i\sum_{k}\lambda_{k,j}^{\ast}(t)e^{-i\omega_{k}t}\Psi_{k}(0)/\sqrt{L},\ \xi_{j}(t)=-i\sum_{k}\lambda_{k,j}^{\ast}(t)e^{-i\omega_{k}t}a_{k}(0)/\sqrt{L}$. $\lambda_{k,j}(t)$ is defined by $\lambda_{k,j}(t)=\lambda_{k,j}\exp[-i\bm{\mathcal{A}}(t)\cdot\bm{r}_{j}]$. We have also introduced local densities:

[TABLE]

The last line in Eq. (8) and Eq. (9) describes a dissipation to the reservoir and a noise from the reservoir.

III Equation of motion in the weak coupling limit

We solve Eqs. (8) and (9) in both weak on-site interaction and weak system-reservoir coupling limit. We here write down the equation of motion in the interaction picture. Without the system-reservoir coupling and the on-site interaction, the solutions of Eqs. (8) and (9) are $\Phi_{\bm{p}}(t)=\exp[-iE_{\bm{p}}t-i\theta_{\bm{p}}(t)]\Phi_{\bm{p}}(0)$ and $b_{\bm{p}}(t)=\exp[-iE_{\bm{p}}t-i\theta_{\bm{p}}(t)]b_{\bm{p}}(0)$, where $\Phi_{\bm{p}}$ and $b_{\bm{p}}$ are the Bloch state:

[TABLE]

The quasi-energy of the system $E_{\bm{p}}$ and the time dependent phase $\theta_{\bm{p}}(t)$ are

[TABLE]

where $\mathcal{J}_{\alpha}^{(n)}$ is the Fourier coefficient of $\exp[-i\mathcal{A}_{\alpha}(t)]$ and $\mathcal{J}_{\alpha}^{(0)}=|\mathcal{J}_{\alpha}^{(0)}|\exp[-iq_{\alpha}]$. If the time reversal symmetry of the Floquet effective Hamiltonian is broken, we have $q_{\alpha}\neq 0,\pi$Sacha et al. (2012); Struck et al. (2012). $b_{\bm{p}}(t)$ is an annihilation operator of a Floquet state.

We define the interaction picture in a usual way Bilitewski and Cooper (2015):

[TABLE]

We take the initial state of the macroscopic Bose field as $\Phi_{\bm{p}}(0)=\Phi\delta_{\bm{p},\bm{q}}$, where $\bm{q}$ labels the smallest energy Bloch state of the system (i.e. $\text{min}_{\bm{p}}E_{\bm{p}}=E_{\bm{q}}$). We also assume that the initial state of the system does not have anomalous average, that is, $\braket{b_{\bm{p}}(0)b_{\bm{k}}(0)}=0$. In the following, the frequency of the external force is taken in the range of $\Delta<\Omega<2\Delta$, where $\Delta=\text{max}_{\bm{p},\bm{k}}|E_{\bm{p}}-E_{\bm{k}}|$ is the single particle band width of the system. Therefore, there is no resonance in the single particle level. We define sets $\mathcal{S}(\omega):=\{\bm{p}\in\mathcal{K}^{d}|E_{\bm{p}}-E_{\bm{q}}=\omega\}$ and $\mathcal{R}=\{\bm{p}\in\mathcal{K}^{d}|E_{\bm{p}}-E_{\bm{q}}=\Omega/2+O(U)\}$, where $\mathcal{K}^{d}$ is the discrete $d$-dimensional quasi-momentum space. $\mathcal{S}(\omega)$ describes a set of states whose quasi-energy is $\omega$. In the weak system-reservoir coupling and the weak on-site interaction limit, $\mathcal{R}=\emptyset$ or there exists a $\omega=E_{\bm{p}^{\ast}}\ (\bm{p}^{\ast}\in\mathcal{K}^{d})$ that satisfies $\mathcal{R}=\mathcal{S}(\omega)\neq\emptyset$. The parametric instability occurs in $\bm{p}\in\mathcal{R}$.

Let us now consider the weak system-reservoir coupling limit (van Hove limitHove (1955)) as done in many of previous papersShirai et al. (2015); Liu (2015); Iadecola and Chamon (2015); Shirai et al. (2016); Iwahori and Kawakami (2016). Here, we also take the weak on-site interaction limit. This treatment enables us to discuss the effect of dissipation on the heating caused by the on-site interaction. The meaning of the weak system-reservoir coupling limit and the weak on-site interaction limit is as follows: We define the dissipation rate $\gamma=\max_{\bm{p}}\gamma_{\bm{p}}$ and the energy spacing $\delta=\text{min}_{E_{\bm{p}}\neq E_{\bm{k}}}|E_{\bm{p}}-E_{\bm{k}}|$, where $\gamma_{\bm{p}}=\int_{0}^{\infty}dt\sum_{n,m,\bm{p}^{\prime}}|\Gamma^{nm}_{\bm{p}\bm{p}^{\prime}}(t)|$ is the dissipation rate of the Bloch state $\bm{p}$ and $\Gamma^{nm}_{\bm{p}\bm{p}^{\prime}}(t)$ is defined in Eq (33) in Appendix. Then, take $\lim_{\gamma t,\ Ut\rightarrow\infty}\lim_{\gamma/\delta,\ U/\delta\rightarrow 0}$. In this limit, we obtain the following equation of motions in the interaction picture (for details of the derivation, see Appendix):

[TABLE]

where $u(t)=U(n_{0}(t)+\tilde{n}(t))$ and $\tilde{u}(t)=2U(n_{0}(t)+\tilde{n}(t))$,

[TABLE]

We have here defined the mean fields:

[TABLE]

$h^{(n)}_{\bm{p},\bm{k}}$ is the Fourier coefficient of $\exp[i(\theta_{\bm{p}}(t)+\theta_{2\bm{q}-\bm{p}}(t)-\theta_{\bm{k}}(t)-\theta_{2\bm{q}-\bm{k}}(t))]$. The definition of all the coefficients is summarized in Appendix. The macroscopic Bose fields whose quasi-momentum does not equal to $\bm{q}$ are zero at any time. The first term in Eq. (14) describes the normal mean field and only results in the energy shift. The coefficients in the second term $v(t)$ and $\tilde{v}_{\bm{p}}(t)$ are written by using the square of the macroscopic Bose field (not the absolute value) and the anomalous average of the Bose operators. This term describes a process of the pair creation or annihilation of the Bogoliubov phonon. Note that there are two processes of the pair creation (annihilation) of Bogoliubov phonon. One is the ordinary pair creation (annihilation) which the photon does not affect, the other is the pair creation (annihilation) in which a photon is absorbed (emitted). The latter process is a cause of the parametric instability in this system Bukov et al. (2015); Lellouch et al. (2016). When the time reversal symmetry of the Floquet effective Hamiltonian is not broken, the photon mediated pair creation and annihilation process does not occur in the weak on-site interaction limitLellouch et al. (2016) (i.e. $h^{(1)}_{\bm{p},\bm{k}}=0$ in our expression). We can break the time reversal symmetry of the Floquet effective Hamiltonian by applying a multichromatic forceSacha et al. (2012); Struck et al. (2012), and therefore we consider the external force which breaks the time reversal symmetry of the Floquet effective Hamiltonian.

IV Stability of the Floquet state

We assume that the system is in the steady state of which the amplitude of the macroscopic Bose field $n_{0}(t)=|\Phi_{\bm{q}}|^{2}/N^{d}$ and the number of the excitations $\tilde{n}=\sum_{\bm{p}}\braket{b^{\dagger}_{\bm{p}}(t)b_{\bm{p}}(t)}/N^{d}$ are time independent. We take the energy of the macroscopic Bose field in the reservoir as $\omega_{k_{0}}=E_{\bm{q}}+(n_{0}+2\tilde{n})U$ and the frequency of the external force as $\Omega=2(E_{\bm{p}}-E_{\bm{q}}+n_{0}U)\ (\bm{p}\in\mathcal{R})$. If there is a detuning, we can eliminate the detuning by taking a gauge transformation of $\Phi_{\bm{q}}$ and $b_{\bm{p}}$ respectively.

We apply the following gauge transformation:

[TABLE]

Then, we have the equation of motion as

[TABLE]

We have also transformed $m_{0}$ and $\tilde{m}_{\bm{p}}$, $\Gamma_{\bm{p}\bm{k}}(t)$, $\Xi^{(I)}_{\bm{q}}(t)$, $\xi^{(I)}_{\bm{p}}(t)$ in the same way as Eq. (17). $g$ and $\tilde{g}_{\bm{p}}$ are defined as $g=U[\tilde{m}_{\bm{q}}+\sum_{\bm{p}\in\mathcal{R}}(h^{(1)}_{\bm{p},\bm{q}})^{\ast}\tilde{m}_{\bm{p}}]$ and $\tilde{g}_{\bm{p}\in\mathcal{R}}=U[(m_{0}+\tilde{m}_{\bm{q}})h^{(1)}_{\bm{p},\bm{q}}+\sum_{\bm{k}\in\mathcal{R}}h^{(0)}_{\bm{p},\bm{k}}\tilde{m}_{\bm{k}}]$, $\tilde{g}_{\bm{p}\not\in{R}}=0$. We here also assume that the anomalous average $\tilde{m}_{\bm{p}}$ and the macroscopic Bose field $m_{0}$ are time independent. We will see that this assumption is also satisfied if $n_{0}$ and $\tilde{n}$ are time independent. Then, we can solve Eqs. (18) by the Laplace transformation. Here we solve the equation of motion for $b^{(I)}_{\bm{p}\in\mathcal{R}}$. The other equation of motion can be solved in the same way. By applying the Laplace transformation, the equation of motion is transformed as

[TABLE]

Functions in the frequency space are specified by their argument $z$ while those in the time space by $t$. We write the number of elements in $\mathcal{R}$ as $n$ and write each element as $\bm{p}_{i}\in\mathcal{R}\ (i=1,2,\cdots,n)$. We define $\bm{b}^{(I)}(z)=(b^{(I)}_{\bm{p}_{1}}(z),\cdots,b^{(I)}_{\bm{p}_{n}}(z),b^{(I)\dagger}_{\bm{p}_{1}}(z),\cdots,b^{(I)\dagger}_{\bm{p}_{n}}(z))^{\text{T}}$ and $\bm{b}=(b_{\bm{p}_{1}},\cdots,b_{\bm{p}_{n}},b^{\dagger}_{\bm{p}_{1}},\cdots,b^{\dagger}_{\bm{p}_{n}})^{\text{T}}$, $\bm{\xi}(z)=(\xi^{(I)}_{\bm{p}_{1}}(z),\cdots,\xi^{(I)}_{\bm{p}_{n}}(z),\xi^{(I)\dagger}_{\bm{p}_{1}}(z),\cdots,\xi^{(I)\dagger}_{\bm{p}_{n}}(z))^{\text{T}}$. Then, we rewrite Eq. (19) in the matrix form:

[TABLE]

where the matrix $M(z)$ is defined through Eq. (19). In this expression, the solution is obtained in the following way:

[TABLE]

where $z_{l}$ is the solution of $\det M(z)=0$ and $\epsilon_{\bm{p}}^{(n)}=E_{\bm{p}}+2(n_{0}+\tilde{n})U+n\Omega$. We here introduce the energy cutoff of the system-reservoir coupling $\omega_{c}$ Shirai et al. (2016); Iwahori and Kawakami (2016) and also consider the limit of $\gamma/\omega_{c},\ U/\omega_{c}\rightarrow 0$. When we solve $\det M(z)=0$, we can take $\Gamma_{\bm{p}\bm{p}^{\prime}}(z)$ as $\Gamma_{\bm{p}\bm{p}^{\prime}}(0)$ because the solution of $\det M(z)=0$ are $z=O(\gamma)+O(U)$ and we have taken the limit of $\gamma/\omega_{c},\ U/\omega_{c}\rightarrow 0$. Therefore, $\det M(z)=0$ reduces to a $2n$th polynomial equation and there are $2n$ solutions. The second term in Eq. (21) describes the steady state. If the the first term decays (i.e. ${\rm Re}[z_{l}]<0,\ \forall l$), we obtain the steady state solution in the long time limit $z_{l}t\rightarrow\infty$ as

[TABLE]

Without dissipation (i.e. $U/\gamma\rightarrow\infty$), we have $z_{l}=\pm|\tilde{g}|$, and thereby we cannot obtain the steady state due to the occurrence of the parametric instabilityBukov et al. (2015); Lellouch et al. (2016). In this regime, a number of excitations (Bogoliubov phonons whose momenta are in $\mathcal{R}$) are created by the external force and the mean field analysis is not applicable in the long time steady stateBukov et al. (2015); Lellouch et al. (2016). On the other hand, we have ${\rm Re}[z_{l}]<0$ without on-site interaction (i.e. $U/\gamma\rightarrow 0$) and the steady state is written in the same way as our previous results Iwahori and Kawakami (2016). Therefore, there is a point where the stability and instability changes in $U/\gamma\in\mathbb{R}_{>0}$ and the parametric instabilities can be suppressed by the dissipation. We calculate this value in Sec. V. We can prove the same thing for the other equation of motion.

In the steady state, the time evolution of the operator $b^{(I)}_{\bm{p}}(t)$ is written in the following way:

[TABLE]

where each coefficient is determined by $M^{-1}(z)$ and $z_{l}$. $b^{(I)}_{\bm{p}}(t)$ is written in terms of $a^{\dagger}_{k}$ as well as $a_{k}$. This is because there is a coupling between $b^{(I)}_{\bm{p}}(t)$ and $b^{(I)\dagger}_{2\bm{q}-\bm{p}}(t)$. This coupling is induced by the on-site interaction $U$, and thereby the second term is $O(U/\gamma)$. By using Eq. (23), the occupation number of the Floquet state and the anomalous average can be obtained as

[TABLE]

where $w^{(n)}_{\bm{p}}\geq 0$ and $\tilde{w}^{(n)}_{\bm{p}}\in\mathbb{C}$ determine the contribution from the $n$th Floquet sideband $E_{\bm{p}}+n\Omega$. $c_{\bm{p}}\geq 0$ and $\tilde{c}_{\bm{p}}\in\mathbb{C}$ are independent of the temperature of the reservoir. The behavior of the first term in Eq. (24) is the same as many of the previous resultsSeetharam et al. (2015); Iadecola and Chamon (2015); Iadecola et al. (2015); Iwahori and Kawakami (2016). The second term comes from the on-site interaction effect and describes small number of excitations created by the external force. Even if we take the zero temperature limit, the second term remains finite while the first term goes to zero. This means that the reservoir can prevent the system from heating but there remains small number of excitations in the system. We will show in the next section, $c_{\bm{p}}$ increases as the on-site interaction becomes stronger.

The occupation number of the Floquet states $\bm{p}\not\in\mathcal{R}$ can be obtained in the same way:

[TABLE]

By utilizing our previous resultsIwahori and Kawakami (2016), when the energy cutoff of the system-reservoir coupling $\omega_{c}$ is smaller than the frequency of the external force $\Omega$, we have $w^{(0)}_{\bm{p}}\simeq 1$ and $w^{(n\neq 0)}_{\bm{p}}\ll 1$ for states $E_{\bm{p}}<\omega_{c}$. Therefore the occupation number of the Floquet state whose quasi-energy is smaller than the energy cutoff of the system-reservoir coupling is effectively written by the Bose distribution function:

[TABLE]

The other time correlation functions have the same structure. As we can see, there are no excitations in states $\bm{p}\not\in\mathcal{R}$. This is because the origin of the excitation $c_{\bm{p}}$ is the coupling between $b_{\bm{p}}$ and $b^{\dagger}_{2\bm{q}-\bm{p}}$ via a photon assist (i.e. pair creation/annihilation of Bogoliubov phonon with one photon absorption/emission) in this case. If we consider a stronger on-site interaction than this case, the scattering of the quasi-particles are induced. This effect causes the diffusion of excitations to all the states and may enable us to define the effective temperature which is obtained in the previous papersBreuer et al. (2000); Ketzmerick and Wustmann (2010); Shirai et al. (2016).

From these observations, we conclude that we can prevent periodically driven quantum systems from heating by introducing dissipation which overwhelms the heating rate. In the realized steady state, we obtain low energy properties of the Floquet effective Hamiltonian. Recent studies on periodically driven quantum isolated systems have revealed that there are two phases of heatingChandran and Sondhi (2016); Weidinger and Knap (2016): The initial stage of heating is governed by the parametric instabilities. Then, the system reaches a quasi-steady state whose energy absorption rate is small. After the quasi-steady state, the system starts to absorb energy again and goes to the infinite temperature state. This heating is governed by the scattering of the excitations which are created in the initial stage of heating. Comparing these previous results of isolated systems, our results demonstrate that the initial stage of heating is prevented and quantum phases which are expected to be realized from the Floquet effective Hamiltonian will be stabilized in the prethermal Floquet steady stateBukov et al. (2015); Lellouch et al. (2016). Note that our results do not directly mean that we can avoid the heating to the infinite temperature state. However, from the discussion above, our results suggest that dissipation will also suppress the second stage of heating which is caused by the diffusion of the excitations to all the states because a number of excitations are suppressed by the dissipation and the dissipation channels exist in all the states.

V Example

Here, we show an example of the above discussion. We will discuss the stability of the steady state and calculate the number of excitations defined in Eq. (24). We consider a one-dimensional system and the external force is taken as $E(t)=E_{1}\cos(\Omega t)+E_{2}\cos(2\Omega t+\phi)$. In this case, the Fourier coefficient of $\exp[-i\mathcal{A}(t)]$ is obtained asSacha et al. (2012)

[TABLE]

where $\zeta_{1}=E_{1}/\Omega$, $\zeta_{2}=E_{2}/\Omega$ and $J_{n}(\zeta)$ is the $n$th Bessel function. The momentum of the BEC is obtained by $q=-\arg(\mathcal{J}^{(0)})$. The system reservoir coupling which is defined in Sec. II is taken as $\lambda_{k,j}=\lambda F(\omega_{k})\delta_{j,0}$ and $D(\omega)F(\omega)^{2}=\theta(\omega)/(1+(\omega/\omega_{c})^{6})$, where $\theta(x)$ is the step function. We define the dissipation rate by $\gamma=|\lambda|^{2}$. For simplicity, we solve Eq. (18) within the Bogoliubov approximation. We only consider the equation of motion for the state in $\mathcal{R}$ because the parametric instability does not occur in the other equation of motion. Then, the equation of motion for the state $p\in\mathcal{R}$ is written as

[TABLE]

The sum in the third term $\sum_{p^{\prime}=p,2q-p}$ results from the degeneracy $E_{p}=E_{2q-p}$. For simplicity, we neglect the sum for $p^{\prime}=2q-p$ and the equation of motion is

[TABLE]

where $u=m_{0}Uh^{(1)}_{p,q}$. Formally, this treatment is justified by applying a perturbation which breaks degeneracy of $p,2q-p\in\mathcal{R}$ and keeps $E_{p}+E_{2q-p}-2E_{q}+2n_{0}U=\Omega$. Then, after taking the weak system-reservoir coupling and the weak on-site interaction limit, we switch off the perturbation. Eq. (29) can be solved in the same way as Eq. (23) and (24). The solution of $\det M(z)=0$, which determines the stability of the system is explicitly written as

[TABLE]

where $\Gamma_{p,p}=\Gamma_{p,p}(z=0)$ and ${\rm Re}[\Gamma_{p,p}],\ {\rm Re}[\Gamma_{2q-p,2q-p}]>0$. If ${\rm Re}[z_{1}]$ and ${\rm Re}[z_{2}]$ are both negative, the steady state is stable. Note that ${\rm Re}[z_{1}]$ is always negative, therefore we discuss only ${\rm Re}[z_{2}]$. The figure of ${\rm Re}[z_{2}]$ and $c_{\bm{p}}$ as a function of $U/\gamma$ is shown in Fig. 1 and 2. The parameters are taken as $(\zeta_{1},\zeta_{2},\phi,N,t_{0}/\Omega,\omega_{c}/\Omega,n_{0})=(1.7965,1,\pi/2,101,0.5213,0.5,1)$. ${\rm Re}[z_{2}]$ gradually increases as $U/\gamma$ increases and at $U/\gamma\simeq 0.82$, ${\rm Re}[z_{2}]$ crosses zero and the steady state becomes unstable. The number of excitations which are created by the external force also monotonically increases as the on-site interaction becomes strong. At the stable/unstable transition point, the number of excitations $c_{\bm{p}}$ diverges. Similar monotonic behaviors of ${\rm Re}[z_{2}]$ and $c_{\bm{p}}$ can be observed in other parameter sets, but the stable/unstable transition point is different depending on the parameters employed. The stable/unstable transition point is determined by the dissipation rate of states $p\in\mathcal{R}$ ($\Gamma_{p,p}$ and $\Gamma_{2q-p,2q-p}$). Therefore, the transition point especially depends on the details of the system reservoir coupling $\lambda_{k,j}$, which determines the dissipation rate of each Bloch state $p\in\mathcal{K}$.

VI Summary

To summarize, we have investigated steady states of a periodically driven dissipative Bose-Hubbard model especially focusing on whether the dissipation can prevent the system from heating. We have revealed that the parametric instability which is considered as an initial stage of the heatingChandran and Sondhi (2016); Weidinger and Knap (2016) can be prevented if the dissipation is stronger than the on-site interaction and stabilizes the steady states. Our findings demonstrate that the low temperature properties of the Floquet effective Hamiltonian are obtained in the prethermal Floquet steady state if the dissipation overwhelms the heating rate. Our results also elucidate that the excitations which are created in the initial stage of heating are suppressed. This fact suggests that the second stage of heating which is governed by the scattering of the excitationsWeidinger and Knap (2016) will be also suppressed by introducing dissipation, and thus the heating to the infinite temperature state will be avoided. In order to check this suggestion, we have to go beyond the weak system-reservoir coupling and the weak on-site interaction limit. This is left for a future study.

Acknowledgements.

This work was partly supported by a Grand-in-Aid for Scientific Research on Innovative Areas (JSPS KAKENHI Grant No. JP15H05855) and also JSPS KAKENHI (No. JP16K05501).

Appendix A Derivation of eq. (14)

We derive the equation of motion in the weak on-site interaction and the weak system reservoir coupling limit. In the interaction picture defined by Eq (13), the equation of motion for the Bloch state is obtained in the following way:

[TABLE]

where $n_{\bm{k},\bm{k}^{\prime}}$, $m_{\bm{k},\bm{k}^{\prime}}(t)$, $\tilde{n}_{\bm{k},\bm{k}^{\prime}}(t)$, $\tilde{m}_{\bm{k},\bm{k}^{\prime}}(t)$ are the mean fields in the interaction picture,

[TABLE]

The dissipation $\Gamma^{nm}_{\bm{p}\bm{p}^{\prime}}(t)$ and the noise $\Xi^{(I)}_{\bm{p}}(t)$, $\xi^{(I)}_{\bm{p}}(t)$ are defined by

[TABLE]

$h^{(m)}_{\bm{p},\bm{p}^{\prime},\bm{k},\bm{k}^{\prime}}$ is the Fourier coefficient of $\exp[i(\theta_{\bm{p}}(t)+\theta_{\bm{k}}(t)-\theta_{\bm{p}^{\prime}}(t)-\theta_{\bm{k}^{\prime}}(t))]$. As we mentioned in the main text, the frequency of the external force is taken in the regime of $\Delta<\Omega<2\Delta$, that is, there is no resonance in the single particle level. In the weak system reservoir coupling and the weak on-site interaction limit, we can show that the only $n=m,\ E_{\bm{p}}=E_{\bm{p}^{\prime}}$ term in the dissipation term contribute to the dynamics. We can also show similar things for the on-site interaction term. Then, the equation of motion in the interaction picture reads

[TABLE]

$u(t)$, $v(t)$, $\tilde{u}$, $\tilde{v}_{\bm{p}}(t)$ are defined in the main text (see Eq. (15)). The dissipation term is $\Gamma_{\bm{p}\bm{p}^{\prime}}(t)=\sum_{n\in\mathbb{Z}}\Gamma^{nn}_{\bm{p}\bm{p}^{\prime}}(t)\exp[i(E_{\bm{p}}+n\Omega)t]$.