Synthesizing lattice structures in phase space

The quantum phenomenona in periodic structures have been of great interest since the beginning of modern solid state physics. Some phenomena which are difficult to observe in natural solid crystals, such as Bloch oscillations [2, 3] and Wannier–Stark ladders [4], have been observed in superlattices [5, 6] and optical lattices [7, 8], which are artificial periodic structures in position space [9–16]. Meanwhile, novel periodic structures in time domain (time crystals [17–19]) and in momentum space (superradiance lattice [20]), have been proposed to be new platforms to study the rich physics in periodic structures. In a recent paper [1], we introduced the idea of phase space crystals, i.e., a lattice structure in phase space created by breaking the continuous phase rotational symmetry via a driving field. Phase space crystals have the key advantage of being conveniently tunable in experiments through changes in the driving field and may provide a new platform to simulate condensed matter phenomena.

In our previous work of phase space crystals [1], we used the model of ultracold atoms trapped in a time-dependent power-law potential, i.e., $\sim {x}^{n}\mathrm{cos}({\omega }_{d}t)$, to illustrate our idea. However, the power-law trapping [21–25] is technically difficult to be realized in experiments. In this work, we present a realistic model, i.e, the ultracold atom in driven optical lattice to realize phase space crystals. In this model, the power-law driving is replaced by a cosine-type driving, i.e., $\sim \mathrm{cos}({kx}+{\omega }_{d}t)$. The model proposed here synthesizes a more complex phase space lattice structure than that produced in [1]. We analyze the quasienergy band structure and identify the artificial (magnetic) gauge field in phase space, which is a result of the noncommutative geometry of phase space. We investigate the dissipative dynamics and study the thermal properties which can be observed in experiments.

The article is organized as follows. In section 2, we describe our model of ultracold atoms in driven optical lattice and introduce the Floquet Hamiltonian under rotating wave approximation(RWA). We scale various experimental parameters into three dimensionless parameters, i.e., the effective Planck constant λ, the scaled detuning and the scaled driving strength μ. In the case of zero detuning $\epsilon =0$, the properties of our system only depend on the effective Planck constant λ. In section 3, we first analyze the symmetries of phase space lattice produced by our model. Then, we introduce the quasinumber theory, which is an application of Bloch theorem in solid state theory to phase space crystals. In section 4, we calculate the quasienergy band structure based on the tight-binding model. We calculate the quantum tunnelling rate by developing WKB theory. We find the tunnelling rate is a complex number which can be explained by the noncommutative geometry of phase space. In section 5, we investigate the dissipative dynamics by the method of master equation at finite temperature. To test our theory, we study the thermal properties which can be observed in the experiments. In section 6, we first discuss a gap to gapless transition of quasienergy bands. Then we discuss the possibility to realize phase space crystals in circuit-QED system. In section 7, we summarize our work and give an outlook of future work.

The model we propose here can be realized by ultracold atoms trapped in a time-dependent optical lattice. The Hamiltonian is given by

$$\begin{eqnarray}&&H(t)=\displaystyle \frac{{p}^{2}}{2m}+\displaystyle \frac{1}{2}m{\omega }^{2}{x}^{2}+2A\mathrm{cos}({kx}+{\omega }_{d}t).\end{eqnarray} \tag{ 1 }$$

Here, the parabolic term is the harmonic confinement potential of ultracold atoms, which can be created by a gaussian beam profile of a laser [26] or introduced by another external field. As sketched in figure 1, the characteristic length of the ground state in the confinement potential is $b\equiv 2\pi \sqrt{{\rm{\hslash }}/(m\omega )}$. Experimentally the optical lattice is created by the interference of two counter-propagating laser beams, which form an optical standing wave with period $d\equiv 2\pi /k$. The ultracold atoms are trapped by the interaction between the laser light field and the oscillating dipole moment of atoms induced by the laser light [27]. We can drive the optical lattice simply by tuning the phase difference of the two laser beams linearly as described by Hamiltonian (1). Effectively, this creates a propagating optical lattice with a velocity of ${\omega }_{d}/k$. An important parameter is $\lambda \equiv {(b/d)}^{2}={\rm{\hslash }}{k}^{2}/(m\omega )$, which defines the 'quantumness' of our system. The parameter λ is large in the quantum regime and goes to zero in the semiclassical limit. We emphasize that the optical potential is time-dependent and the confinement potential also plays an important role. Thus, our system does not have spatial periodicity and the Bloch theory in real space does not apply directly for the Hamiltonian (1).

Zoom In Zoom Out Reset image size

We are interested in the regime near the high-order resonant condition ${\omega }_{d}\approx n\omega$ with a large integer $n\gg 1$. For the duration of this paper we will use n = 30. The detuning $\delta \omega \equiv \omega -{\omega }_{d}/n$ is much smaller than the natural frequency ω. We perform a unitary transformation of the Hamiltonian $H(t)$ via the operator $\hat{U}={{\rm{e}}}^{{\rm{i}}({\omega }_{d}/n){\hat{a}}^{\dagger }\hat{a}t}$, where $\hat{a}$ is the annihilation operator of the oscillator. In the spirit of the rotating wave approximation (RWA), we drop the fast oscillating terms and arrive at the time-independent Hamiltonian (see the detailed derivations in appendix A)

$$\begin{eqnarray}&&\hat{g}=\lambda \epsilon ({\hat{a}}^{\dagger }\hat{a}+\displaystyle \frac{1}{2})+\mu \left[{\left(\displaystyle \frac{\lambda }{2}\right)}^{-\displaystyle \frac{n}{2}}{{\rm{e}}}^{-\displaystyle \frac{\lambda }{4}-{\rm{i}}\displaystyle \frac{n\pi }{2}}{\hat{a}}^{n}\ {L}_{{\hat{a}}^{\dagger }\hat{a}}^{-n}\left(\displaystyle \frac{\lambda }{2}\right)+{\rm{h.c.}}\right].\end{eqnarray} \tag{ 2 }$$

In the context of Floquet theory[28–33], $\hat{g}$ is called quasi-energy Hamiltonian , which has been scaled by the energy $m{(\omega /k)}^{2}={\rm{\hslash }}\omega /\lambda$. The eigenvalues of $\hat{g}$ are called quasienergies. The parameters $\epsilon \equiv \delta \omega /\omega$ and $\mu \equiv \lambda A/({\rm{\hslash }}\omega )$ are the dimensionless detuning and driving strength respectively. Functions ${L}_{{\hat{a}}^{\dagger }\hat{a}}^{-n}(\bullet )$ are the generalized Laguerre polynomials, as a function of the photon number ${\hat{a}}^{\dagger }\hat{a}| k\rangle =k| k\rangle$, where $| k\rangle$ are the Fock states.

3.1. Symmetries

In the following, we are particularly interested in the resonant condition, i.e., the detuning is zero $\epsilon =0$. Without loss of generality, we set the scaled driving strength to unity, i.e., μ = 1. In this case, the RWA Hamiltonian (2) has two new symmetries which are not visible in the original Hamiltonian (1). To visualize them, we replace the operator $\hat{a}$ by a complex number in the semiclassical limit and plot the quasienergy g in the phase space spanned by Re[a] and Im[a]. As displayed in figure 2(a) and figure 2(b), we first see the discrete angular symmetry$g(\theta )=g(\theta +2\pi /n)$. Additionally we have the chiral symmetry $g(\theta )=-g(\theta +\pi /n)$, which divides the whole lattice structure into two identical sublattices as indicated in figure 2(b) by the different colors.

Zoom In Zoom Out Reset image size

To describe the two symmetries in quantum mechanics, we define a unitary operator ${\hat{T}}_{\tau }={{\rm{e}}}^{-{\rm{i}}\tau {\hat{a}}^{\dagger }\hat{a}}$ with the properties ${\hat{T}}_{\tau }^{\dagger }\hat{a}{\hat{T}}_{\tau }=\hat{a}{{\rm{e}}}^{-{\rm{i}}\tau }$ and ${\hat{T}}_{\tau }^{\dagger }{\hat{a}}^{n}{\hat{T}}_{\tau }=\hat{a}{{\rm{e}}}^{-{\rm{i}}{n}\tau }$. Since ${\hat{T}}_{\tau }^{\dagger }{\hat{a}}^{\dagger }\hat{a}{\hat{T}}_{\tau }={\hat{a}}^{\dagger }\hat{a}$, it is not difficult to check that the RWA Hamiltonian (2) is invariant under this discrete transformation

$$\begin{eqnarray}&&{T}_{\tau }^{\dagger }\hat{g}{T}_{\tau }=\hat{g}\ \ \ \mathrm{for}\ \ \ \tau =\displaystyle \frac{2\pi }{n}.\end{eqnarray} \tag{ 3 }$$

We call this symmetry discrete phase translation symmetry. The chiral symmetry follows from the fact ${T}_{\tau }^{\dagger }\hat{g}{T}_{\tau }=-\hat{g}$ for $\tau =\pi /n$. The chiral symmetry suggests that the two sublattices are symmetric with respect to g = 0, except a phase shift $\theta \to \theta +\pi /n$. The angular symmetry indicates it is convenient to introduce the radial and angular operators $\hat{r}$ and $\hat{\theta }$ via

$$\begin{eqnarray}&&\hat{a}\equiv {{\rm{e}}}^{-{\rm{i}}\hat{\theta }}\hat{r}/\sqrt{2\lambda }\ \ \ \ \mathrm{and}\ \ \ \ \ {\hat{a}}^{\dagger }\equiv \hat{r}{{\rm{e}}}^{{\rm{i}}\hat{\theta }}/\sqrt{2\lambda }.\end{eqnarray} \tag{ 4 }$$

They obey the commutation relation

$$\begin{eqnarray}&&[{\hat{r}}^{2},{{\rm{e}}}^{{\rm{i}}\hat{\theta }}]=2\lambda {{\rm{e}}}^{{\rm{i}}\hat{\theta }}\end{eqnarray} \tag{ 5 }$$

where λ plays the role of a dimensionless effective Planck constant.

3.2. Phase space lattice

In the semiclassical limit $\lambda \to 0$, the quantum Hamiltonian $\hat{g}$ can be written in its classical form (see more details in appendix A)

$$\begin{eqnarray}&&g=\displaystyle \frac{1}{2}\epsilon {r}^{2}+2\mu {J}_{n}(r)\mathrm{cos}\left(n\theta -\displaystyle \frac{n\pi }{2}\right).\end{eqnarray} \tag{ 6 }$$

Here, we have used the asymptotic property of Laguerre polynomials, i.e.,

$$\begin{eqnarray*}&&\underset{k\to \infty }{\mathrm{lim}}{L}_{k}^{n}(x/k)={k}^{n}{{\rm{e}}}^{\displaystyle \frac{x}{2k}}{x}^{-n/2}{J}_{n}(2\sqrt{x}),\end{eqnarray*}$$

where ${J}_{n}(\bullet )$ is the Bessel function of order n. The angular periodicity comes from the cosine function in equation (6) while the radial structure is created by the Bessel function ${J}_{n}(r)$. A similar situation has recently been studied in voltage biased Josephson junctions [34, 35].

As shown in figure 2(b), the zero lines of quasienergy (i.e., $g(r,\theta )=0$) divide the whole phase space lattice into many small ' cells'. The center of each cell is a stable point corresponding to either a local minimum or a local maximum of $g(r,\theta )$ (see more discussions in section appendix C). The area inside the cell represents the basin of attraction for the stable state in the center. In figure 2(c), we plot the angular periodical structure of the quasienergy g for a fixed radius indicated by the dashed blue line in figure 2(b). In figure 2(d), we show the radial structure of the quasienergy g by plotting it along two angular directions $\theta =0$ and $\theta =\pi /n$. We see the quasienergy oscillates as a function of the radius r in the form of Bessel functions ${J}_{n}(r)$. We divide the whole lattice structure into 'loops', which correspond to ring-like areas in figure 2(b) between two radii satisfying ${J}_{n}(r)=0$. We label them from inside to outside by Roman numerals $I$, ${II}$, ${III}$ and so on as indicated in figure 2(d). The corresponding localized quantum levels inside the loops are labelled by ${Level}\ I$-1, ${Level}\ I$-2, ${Level}\ {II}$-1, etc.

3.3. Quasinumber theory

We diagonalize the quantum Hamiltonian (2) and study the properties of its spectrum. With zero detuning $\epsilon =0$ and driving μ = 1, the spectrum is only determined by the effective Planck constant λ. In figure 3(a) we show the structure of the quasienergy spectrum as function of the parameter $1/\lambda$. It is clear that the quasienergy spectrum is symmetric with respect to g = 0 because of the chiral symmetry. We also see that the gaps are opened for small λ and closed for sufficiently large λ. The transition happens around $\lambda \approx 5$. We will calculate the gaps using WKB theory and discuss the physical mechanism of gap closing in section 4. In figure 3(b) we show the gapless quasienergy spectrum for $\lambda =6$ and the band structure for $\lambda =4$. The band structure comes from the discrete phase translation symmetry.

Zoom In Zoom Out Reset image size

We introduce the quasinumber theory [1] according to the Bloch's theorem. Due to ${\hat{T}}_{\tau }^{\dagger }\hat{g}{\hat{T}}_{\tau }=\hat{g}$ for $\tau =2\pi /n$, the eigenstates ${\psi }_{m}(\theta )$ of the quasienergy Hamiltonian, $\hat{g}{\psi }_{m}(\theta )=g(m){\psi }_{m}(\theta )$, have the form

$$\begin{eqnarray*}&&{\psi }_{m}(\theta )={\varphi }_{m}(\theta ){{\rm{e}}}^{-{\rm{i}}{m}\theta },\end{eqnarray*}$$

with a periodic function ${\varphi }_{m}(\theta +\tau )={\varphi }_{m}(\theta )$. Here, the integer number m is called quasinumber, which is conjugate to the phase θ. It is an analogue of the quasimomentum $\overset{\rightharpoonup }{k}$ in a solid crystal. In figure 3(c), we plot the quasienergy band structure in the reduced Brillouin zone $m\tau \in [0,2\pi )$. We count the bands from the bottom and relabel the eigenstates ${\psi }_{m}(\theta )$ by ${\psi }_{l,m}(\theta )$, where the subscript $l=1,2,\ldots$ indicates the band that the eigenstate belongs to.

To visualize the quasinumber states, we plot the Q-function of state ${\psi }_{\mathrm{1,0}}(\theta )$ in figure 4(a). The Q-function is a quasi-probability distribution in phase space defined by [36]

$$\begin{eqnarray*}&&Q(\alpha ,{\alpha }^{*})\equiv \displaystyle \frac{1}{\pi }{| \langle \alpha | {\psi }_{l,m}(\theta )\rangle | }^{2},\end{eqnarray*}$$

where $| \alpha \rangle$ is the coherent state given by $\hat{a}| \alpha \rangle =\alpha | \alpha \rangle$ or $| \alpha \rangle ={{\rm{e}}}^{-{| \alpha | }^{2}/2}{\displaystyle \sum }_{k=0}^{\infty }{\alpha }^{k}/\sqrt{k!}| k\rangle$. The crystalline structure of Q-function in angular direction reflects the n-fold discrete phase translation symmetry. In figure 4(b), we plot the occupation number statistics of Fock states, i.e., ${{\rm{P}}}_{1,m}(k)={| \langle k| {\psi }_{1,m}(\theta )\rangle | }^{2}$, for quasinumber states with m = 0 and m = 15 in the first band l = 1. As we can see from the probability distribution, the quasinumber states are the superposition of Fock states with photon numbers being multiples of n.

Zoom In Zoom Out Reset image size

From the form of the Q-function we see that the eigenstates of the system ${\psi }_{m}(\theta )$ are delocalized states in phase space, which are superposition of localized states corresponding to the discrete energy levels as indicated in figure 2(d). We label these levels in the first loop by ${Level}\ I-1$, ${Level}\ I-2$ and those in the second loop ${Level}\ {II}-1$ etc. In the semiclassical limit, these quantum levels become classical orbits of iso-quasienergy contours represented by the boundaries of the colored elliptical areas inside each cell shown in figure 2(b). The shapes of these orbits vary in different loops as displayed on the top of figure 6(b).

The formation of quasienergy bands near the bottom can be understood in the frame of the tight-binding model. If we neglect quantum tunnelling, the n localized states in each loop are n degenerate states. If we consider quantum tunnelling, they are broadened and form bands. We can label the bands by the labels of corresponding localized levels, e.g., the bottom band of the whole quasienergy spectrum is ${Band}\ I-1$. We can describe the structure of the l-th tight-binding band approximately by

$$\begin{eqnarray}&&{g}_{l}(m)={E}_{l}-2| {J}_{l}| \mathrm{cos}(m\tau +{\delta }_{l}\tau ).\end{eqnarray} \tag{ 7 }$$

Here, E_l represents the center of the l-th band and the quasienergy of the corresponding localized level. The l-th bandwidth d_l is determined by the amplitude of dynamical tunnelling rate J_l, i.e., ${d}_{l}=4| {J}_{l}|$. From figure 3(c) we see that the bands are NOT symmetric with respect to the center of the Brillouin zone in general. We describe the asymmetry by an asymmetry factor${\delta }_{l}$. The asymmetry factor comes from the fact that the two dimensions of phase space are not commutative. We will calculate the gap, bandwidth and asymmetry factor by WKB theory below.

4.1. Dynamical quantum tunnelling in phase space

From the commutation relation (5), it can be shown that [1]

$$\begin{eqnarray*}&&\left[\displaystyle \frac{1}{2}{\hat{r}}^{2},\hat{\theta }\right]\approx {\rm{i}}\lambda \end{eqnarray*}$$

in the region of $r\gg 1$. We thus can view operators ${\hat{r}}^{2}/2$ and $\hat{\theta }$ as 'coordinate' and 'momentum' respectively, i.e.,

$$\begin{eqnarray*}&&\hat{\theta }\approx -{\rm{i}}\lambda {r}^{-1}\displaystyle \frac{\partial }{\partial r}.\end{eqnarray*}$$

In the semiclassical limit, the variables ${r}^{2}/2$ and θ define the phase space for our WKB calculation. In figure 5(a), we plot the quasienergy g in the range of $\theta \in [-2\pi /n,2\pi /n]$ in the phase space spanned by ${r}^{2}/2$ and θ. For a fixed g, all the branches of classical orbits are obtained from equation (6)

$$\begin{eqnarray}&&{\theta }_{\pm }(r)=\displaystyle \frac{\pi }{2}\pm \displaystyle \frac{1}{n}\mathrm{arccos}\left[\displaystyle \frac{g-\epsilon {r}^{2}/2}{2\mu {J}_{n}(r)}\right]+k\displaystyle \frac{2\pi }{n},\end{eqnarray} \tag{ 8 }$$

where k takes integers 0, 1, 2, $\cdot \cdot \cdot$, and $n-1$. Two real solutions ${\theta }_{\pm }(r)$ together represent one closed classical orbit. There are n identical orbital branches with a $2\pi /n$-shift of θ. From the condition

$$\begin{eqnarray*}&&\left|\displaystyle \frac{g-\epsilon {r}^{2}/2}{2\mu {J}_{n}(r)}\right|\lt 1,\end{eqnarray*}$$

we can determine the boundaries of classical motion for the fixed g as shown by the black closed curves in figure 5(a). Outside the boundaries, the solutions ${\theta }_{\pm }(r)$ have imaginary parts. In figure 5(a), we indicate the maximal boundaries of classical motion by r₁, r₂ and r₃. The region between r₂ and r₃ is the classically forbidden region for a fixed g. In the quantum regime, these states can tunnel into each other. In figure 5(a), we show how the two neighboring ${Level}\ I-1$ states tunnel into each other through the paths in phase space. The optimal tunnelling path is indicated by the white arrows, i.e., one ${Level}\ I-1$ state first tunnels into the nearest region in ${Loop}\ {II}$ across one saddle point (white dot) and then tunnels back to the neighboring ${Level}\ I-1$ state across another saddle point. There are also many other possible tunnelling paths which are indicated by yellow arrows in figure 5(a). But the contributions from these paths are exponentially small compared to the main tunnelling path (see more details in appendix B).

Zoom In Zoom Out Reset image size

4.2. Quasienergy levels and bandwidths

From the WKB theory, we know the phase space area enclosed by the classical orbit is quantized according to the so called Bohr-Sommerfeld quantization condition [37]

$$\begin{eqnarray}S(g) & = & {\displaystyle \int }_{{r}_{1}}^{{r}_{2}}{\theta }_{+}(r){r}{\rm{d}}{r}-{\displaystyle \int }_{{r}_{1}}^{{r}_{2}}{\theta }_{-}(r){r}{\rm{d}}{r}\\ & = & \displaystyle \frac{2}{n}{\displaystyle \int }_{{r}_{1}}^{{r}_{2}}\mathrm{arccos}\left[\displaystyle \frac{g-\epsilon {r}^{2}/2}{2\mu {J}_{n}(r)}\right]{r}{\rm{d}}{r}\\ & = & 2\pi \lambda \left(k+\displaystyle \frac{1}{2}\right),\end{eqnarray} \tag{ 9 }$$

where k takes nonnegative integers. From the above condition we can calculate the quasienergy levels. As shown in figure 5(b), the left subfigure shows several lowest levels calculated using the above quantization condition. We compare our WKB calculations (dashed lines) to the numerical simulations (solid lines) in figure 5(b). The agreements are satisfied. Noticeably, ${Level}\ I$-2 and ${Level}\ {II}$-1 cross each other near $1/\lambda =1.2$. The level crossing has significant effect on the bandwidths as we discuss below.

The width of the l-th band d_l is given by the tunnelling rate J_l, i.e., ${d}_{l}=4| {J}_{l}|$. The amplitude of J_l is given by the integral of the imaginary part of 'momentum' θ in the classical forbidden region ${r}_{2}\lt r\lt {r}_{3}$

$$\begin{eqnarray}&&| {J}_{l}| =\displaystyle \frac{\lambda }{2\pi }{{\left(\displaystyle \frac{{\rm{d}}{S}}{{\rm{d}}{g}}\right)}^{-1}|}_{g={g}_{l}}\mathrm{exp}\left(-\displaystyle \frac{2}{\lambda }{\displaystyle \int }_{{r}_{2}}^{{r}_{3}}\mathrm{Im}[\theta ]\ {r}{\rm{d}}{r}\right).\end{eqnarray} \tag{ 10 }$$

Here, $S(g)$ in the prefactor as function of g is given by the first equality of equation (9). In appendix B, we give a detailed description of the behavior of $\mathrm{Im}[\theta ]$ in the classical forbidden region. Here we just present our results. In figure 5(b) we show the bandwidths of ${Level}\ I$-1 and ${Level}\ I$-2 calculated by equation (10) and compare them to the numerical calculations. There is a cusp in the curve of ${Level}\ I$-2. This happens because of the crossing of ${Level}\ I$-2 and ${Level}\ {II}$-1 which significantly enhances the quantum tunnelling of ${Level}\ I$-2. In this case, we need to consider three interacting levels, i.e., two neighboring ${Level}\ I$-2 states and the medium state of ${Level}\ {II}$-1 as indicated by the closed orbits in figure 5(a). The Hamiltonian of three interacting levels (TIL) is described by the following 3 × 3 matrix

$$\begin{eqnarray}{H}_{{\rm{TIL}}}=\left(\begin{array}{ccc}{g}_{1} & {J}_{12} & {J}_{11}\\ {J}_{12} & {g}_{2} & {J}_{12}\\ {J}_{11} & {J}_{12} & {g}_{1}\end{array}\right).\end{eqnarray} \tag{ 11 }$$

Here g₁, g₂ represent the quasienergies of ${Level}\ I$-2 and ${Level}\ {II}$-1 respectively. Parameter J₁₁ represents the tunnelling rate between the two neighboring ${Level}\ I$-2 states. Parameter J₁₂ represents the tunnelling rate between the state of ${Level}\ I$-2 and the state of ${Level}\ {II}$-1. The tunnelling rate J₁₁ is given by equation (10) by taking $g={g}_{1}$, while the tunnelling rate J₁₂ is given by the following

$$\begin{eqnarray}&&{J}_{12}=\displaystyle \frac{\lambda }{2\pi }{{\left(\displaystyle \frac{{\rm{d}}{S}}{{\rm{d}}{g}}\right)}^{-1}|}_{g={g}_{2}}\mathrm{exp}\left(-\displaystyle \frac{1}{\lambda }{\displaystyle \int }_{{r}_{2}}^{{r}_{3}}\mathrm{Im}[\theta ]\ {r}{\rm{d}}{r}\right).\end{eqnarray} \tag{ 12 }$$

We can get the modified quasienergy levels by diagonalizing the matrix H_TIL. Half of the level spacing ${{\rm{\Delta }}}_{11}$ of the two modified ${Level}\ I$-2 levels gives the effective tunnelling rate between them, i.e., ${J}_{{\rm{eff}}}={{\rm{\Delta }}}_{11}/2$. Therefore, the correct bandwidth of ${Band}\ I$-2 is $4{J}_{{\rm{eff}}}=2{{\rm{\Delta }}}_{11}$. Figure 5(b) shows our theoretical calculations (solid lines) agree with the numerical calculations (dashed lines) very well.

4.3. Band asymmetry and artificial magnetic field

From figure 3(c), we see that the quasienergy bands are not symmetric with respect to the center of the reduced Brillouin zone. The asymmetry is described by the asymmetry factor ${\delta }_{l}$. In the frame of tight-binding approximation, the Bloch eigenstate ${\psi }_{{lm}}(\theta )$ is given by

$$\begin{eqnarray*}&&{\psi }_{{lm}}(\theta )=\displaystyle \frac{1}{\sqrt{n}}\displaystyle \sum _{q=0}^{n-1}{{\rm{e}}}^{{\rm{i}}{m}{q}\tau }{\hat{T}}_{\tau }^{q}{\phi }_{l}(\theta ),\end{eqnarray*}$$

where ${\phi }_{l}(\theta )$ is the wave functions of localized states (which is called Wannier state in solid state physics). The quantum tunnelling rate can be calculated by

$$\begin{eqnarray*}&&{J}_{l}=-\int {[{\hat{T}}_{\tau }{\phi }_{l}(\theta )]}^{*}\hat{g}{\phi }_{l}(\theta )d\theta .\end{eqnarray*}$$

The corresponding quasienergy spectrum of the l-th band then is given by

$$\begin{eqnarray*}&&{g}_{l}(m)={\displaystyle \int }_{0}^{2\pi }{\psi }_{{lm}}^{*}(\theta )\hat{g}{\psi }_{{lm}}(\theta ){\rm{d}}\theta \approx {E}_{l}-{J}_{l}{{\rm{e}}}^{{\rm{i}}{m}\tau }-{J}_{l}^{*}{{\rm{e}}}^{-{\rm{i}}{m}\tau }.\end{eqnarray*}$$

The band asymmetry comes from the fact that quantum tunnelling rate J_l in driven systems is generally a complex number [1, 38], i.e., ${J}_{l}=| {J}_{l}| {{\rm{e}}}^{-{\rm{i}}{\delta }_{l}\tau }$, and the phase parameter ${\delta }_{l}$ is exactly the asymmetry factor. We can calculate the phase ${\delta }_{l}$ using the WKB theory discussed above.

In fact, when r is approaching one of the roots ${r}^{(0)}$ with ${J}_{n}({r}^{(0)})=0$, from equation (8) we see the amplitude of 'momentum' θ goes to infinity $| \theta ({r}^{(0)})| \to \infty$. This means the WKB approximation breaks down near the root ${r}^{(0)}$ and we need a connecting condition. In the range of ${r}^{(0)}\gg 1$, we can expand the phase translation operator ${\hat{T}}_{\tau }={{\rm{e}}}^{-{\rm{i}}\tau {\hat{a}}^{\dagger }\hat{a}}$ by [1]

$$\begin{eqnarray*}&&{\hat{a}}^{\dagger }\hat{a}\approx {\lambda }^{-1}{({r}^{(0)})}^{2}/2+{\rm{i}}\partial /\partial \theta .\end{eqnarray*}$$

The connecting condition, i.e., the neighboring localized state of ${\phi }_{l}(\theta )$, is given by

$$\begin{eqnarray*}&&{\hat{T}}_{\tau }{\phi }_{l}(\theta )\approx {{\rm{e}}}^{-{\rm{i}}{\lambda }^{-1}{({r}^{(0)})}^{2}\tau /2}{\phi }_{l}(\theta +\tau ).\end{eqnarray*}$$

Thus we get the symmetry factor

$$\begin{eqnarray*}&&{\delta }_{l}={\delta }_{l}^{0}+\displaystyle \frac{1}{2\lambda }{({r}^{(0)})}^{2},\end{eqnarray*}$$

where ${\delta }_{l}^{0}$ is the residual asymmetry beyond WKB calculation and can be removed by redefining the phase translation operator ${\hat{T}}_{\tau }={{\rm{e}}}^{-{\rm{i}}\tau ({\hat{a}}^{\dagger }\hat{a}-{\delta }_{l}^{0})}$. The asymmetry factor ${\delta }_{l}$ is linearly dependent on the parameter $1/\lambda$ with the slope ${({r}^{(0)})}^{2}/2$ differing between bands. In appendix B, we give more detailed discussion on tunnelling paths and show more results about the linear relationship between ${\delta }_{l}$ versus $1/\lambda$ for different bands. If we count ${r}^{(0)}=0$ as the first root of ${J}_{n}(r)$, then the asymmetry factors of bands in the l-th ($l\geqslant 2$) loop are all given by the l-th ($l\geqslant 2$) root of the Bessel function. But the asymmetry factors of the bands in the first loop are determined by the second root of the Bessel function. The reason is that the localized states inside the first loop tunnel through its upper boundary while states in other loops tunnel through lower boundaries.

The fact that the tunnelling amplitudes are complex means there is an artificial magnetic field B_eff in phase space. Imagine we have a loop of atoms forming a one dimensional lattice in real space with magnetic field B across the loop. The magnetic field induces an additional phase to the tunnelling amplitude between neighbored atoms $J=| J| {{\rm{e}}}^{-{\rm{i}}\delta }$, where $\delta \propto B$ is called Peierls phase [39]. Comparing the Peierls phase to the asymmetry factor of the phase space lattice calculated above, we can identify there is an effective magnetic field ${B}_{{\rm{eff}}}\propto 1/\lambda$ in phase space. The coordinate system of a phase space lattice has a noncommutative geometry [40], which is fundamentally different from spatial lattices. It is this noncommutative phase space which creates an artificial magnetic field and is responsible for the asymmetry of the quasienergy band structure.

The above calculation of the quasienergy bandstructure does not consider the dissipative environment due to the quantum and thermal fluctuations. We use the master equation method to describe the dissipative evolution in experiments. Already previously it has been shown that a Lindblad type of master equation [33, 41–44] is sufficient as description,

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}=-\displaystyle \frac{{\rm{i}}}{\lambda }[\hat{g},\rho ]+\kappa (1+\bar{n}){ \mathcal D }[a]\rho +\kappa \bar{n}{ \mathcal D }[{a}^{\dagger }]\rho ,\end{eqnarray} \tag{ 13 }$$

where the time t is dimensionless and scaled by the natural frequency ω. The Lindblad superoperator is defined through ${ \mathcal D }[A]\rho \equiv A\rho {A}^{\dagger }-({A}^{\dagger }A\rho +\rho {A}^{\dagger }A)/2$. Bose distribution $\bar{n}={({{\rm{e}}}^{{\rm{\hslash }}{\omega }_{0}/{k}_{B}T}-1)}^{-1}$ represents temperature and dimensionless damping κ is scaled by ω. The type of Lindblad master equation is widely used in the field of cold atoms [45–47].

Based on the master equation (13), we calculate the density matrix of the stationary distribution, i.e., $\langle k| \rho | k\rangle$, in the basis of the Fock states $\{| k\rangle ,k=0,1,\cdot \cdot \cdot \}$. By the relationship of $k={r}^{2}/(2\lambda )$ from equation (4), we define the probability density $\rho (r)$

$$\begin{eqnarray*}&&\rho (r)\equiv \langle k| \rho | k\rangle \displaystyle \frac{{\rm{d}}{k}}{{\rm{d}}{r}}=\displaystyle \frac{r}{\lambda }\langle k| \rho | k\rangle .\end{eqnarray*}$$

In figure 6(a), we plot $\rho (r)$ for different temperatures $\bar{n}=0$ and $\bar{n}=0.1$. We see that $\rho (r)$ oscillates with radius r. The zero nodes of $\rho (r)$ correspond to the loop boundaries of phase space lattice shown in figure 2(b).

Zoom In Zoom Out Reset image size

The probabilities over the loops are not equally distributed because the quantum heating [48] of each loop is not the same. On the bottom of each loop, the stationary distribution can be described by an effective temperature ${\bar{n}}_{e}$. The localized ground state of each loop can be approximately described by a squeezed state with the squeezing factor u and the corresponding effective temperature is given by (see the detailed calculations in appendix D)

$$\begin{eqnarray*}&&{\bar{n}}_{e}={| u| }^{2}+\bar{n}(2{| u| }^{2}+1).\end{eqnarray*}$$

In our case, as we can see from figure 6(a), the peak of $\rho (r)$ is in the third loop. The reason is the effective temperature of the third loop is lower than the temperatures in other loops. In figure 6(b), we calculate the squeezing factor u and the effective temperature ${\bar{n}}_{e}$ for the first ten loops and compare them to fully numerical simulations. The agreement is very good. Another interesting fact is the squeezing factor u changes from a negative value to a positive value. This means the shape of the squeezed state in each loop is different as displayed by the classical orbits on the top of figure 6(b). The orbital shapes are taken from the plot in figure 2(b). The third orbit is very close to a round circle, which means the squeezing factor $u\approx 0$ and the resulting effective temperature ${\bar{n}}_{e}\approx \bar{n}$ is the lowest one of all the loops. The stationary distribution can be directly measured in the experiments [49].

6.1. Transition from gapped to non-gapped spectrum

From figure 3(a), we see that the quasienergy spectrum undergoes a transition from gapless state to gapped state as the effective Planck constant λ decreases. In figure 7, we plot the behavior of bottom gap Δ near the critical point ${\lambda }_{c}\approx 5$ for n = 20, 30 and 50, respectively. The results are obtained from numerical simulations and the quantity Δ is determined by the largest level spacing. The critical point ${\lambda }_{c}$ can be estimated by the condition that the bottom gap is filled by the widths of neighboring bands. Basically, the value of the critical point depends on the parameter n weakly. This transition may be explored in the study of heat transport in the ion traps [50, 51].

Zoom In Zoom Out Reset image size

6.2. Possibility in circuit-QED systems

The phase space lattice may be realized in circuit-QED systems, i.e., a superconducting cavity coupled to Josephson junctions with the Hamiltonian

$$\begin{eqnarray*}&&{H}_{{\rm{cQED}}}={\rm{\hslash }}\omega {a}^{\dagger }a+2{E}_{J}\mathrm{cos}(4\pi {\rm{e}}{{\rm{\hslash }}}^{-1}{\rm{\Phi }})\mathrm{cos}\varphi .\end{eqnarray*}$$

The Josephson junction can be driven by either a dc voltage [34, 35], which creates $\varphi ={\varphi }_{0}+{\omega }_{d}t$ with ${\omega }_{d}=2{eV}/{\rm{\hslash }}$, or a time-dependent magnetic flux [52]${\rm{\Phi }}={\omega }_{d}t/(4\pi {\rm{e}}{{\rm{\hslash }}}^{-1})$. The effective Planck constant in this case is given by $\lambda =8\pi \omega L/(h/{e}^{2})$, where L is the inductance of the circuit and $h/{e}^{2}\approx 25.8k{\rm{\Omega }}$ is the von Klitzing constant. The typical impedance $\omega L$ of circuit-QED systems using only geometrical inductors and capacitors, can not exceed the characteristic impedance of vacuum ${\mu }_{0}c\approx 376.73{\rm{\Omega }}$ [53], which means that we have $\lambda \lt 0.015$ in circuit-QED systems. However, there are several proposals to realized a super-inductance based on the design of Josephson junction arrays [53, 54] which can increase the impedance significantly up to 35 $k{\rm{\Omega }}$ resulting $\lambda \gt 1$. Thus, it is possible to realize phase space lattices in circuit-QED systems combined with a proper design of Josephson junction arrays.

We have studied a type of phase space crystals, which can be realized by untracold atoms driven by optical lattice. We discussed the physics near the high-order resonant condition, i.e., the driving frequency is close to a large multiple of the frequency of trapping potential. In the framework of Floquet theory, we showed the Hamiltonian under RWA has complex periodic lattice structure in phase space. We analyzed the discrete phase translation symmetry and chiral symmetry of the phase space lattice in the case of zero detuning. Based on the quasinumber theory, we studied the properties of the quasienergy bands and the quasinumber states. We calculated the quasienergy band structure using WKB theory based on the tight-binding model. We discussed the quantum tunnelling paths in phase space and found the quantum tunnelling rate is a complex number due to the noncommutative geometry of phase space. The bands can touch each other and enhance their bandwidths significantly. We investigated the dissipative dynamics of our phase space crystals due to quantum and thermal fluctuations. We found the density distribution $\rho (r)$ oscillates with the loops of phase space lattice. The whole distribution can not be described by a uniform temperature but by a series of effective temperatures in different loops. In the end, we discuss a transition from gapped to gapless quasienergy spectrum and the possibility to realize phase space crystals in a circuit-QED system.

We should clarify that the quasienergy band theory discussed here is based on the single particle picture like the band theory in the solid state physics. Future work may go beyond the single particle approximation and consider the interactions between atoms.

We acknowledge helpful discussions with Dr P Kotetes, Dr P Jin, Prof G Johansson, Prof G Schön, Prof F Marquart and Dr V Peano. We acknowledge financial support from Carl-Zeiss Stiftung. We acknowledge support by Deutsche Forschungsgemeinschaft and Open Access Publishing Fund of Karlsruhe Institute of Technology.

In this section, we give detailed derivation from the time-dependent Hamiltonian (1) to the RWA Hamiltonian (2) and the semiclassical Hamiltonian (6) in the main text. Now, we introduce $a,{a}^{\dagger }$ via $x=\sqrt{{\rm{\hslash }}/(2m\omega )}({a}^{\dagger }+a)$ and $p={\rm{i}}\sqrt{m{\rm{\hslash }}\omega /2}({a}^{\dagger }-a)$. By introducing parameter $\lambda \equiv {\rm{\hslash }}{k}^{2}/(m\omega )$, we map the Hamiltonian (1) to the following

$$\begin{eqnarray}&&H(t)={\rm{\hslash }}\omega {a}^{\dagger }a+2A\mathrm{cos}\left[\sqrt{\displaystyle \frac{\lambda }{2}}({a}^{\dagger }+a)+{\omega }_{d}t\right].\end{eqnarray} \tag{ A.1 }$$

We introduce the scaled coordinate and momentum operators $\hat{Q}=\sqrt{\displaystyle \frac{\lambda }{2}}({a}^{\dagger }+a)$ and $\hat{P}={\rm{i}}\sqrt{\displaystyle \frac{\lambda }{2}}({a}^{\dagger }-a)$ with the noncommutative relationship $[\hat{Q},\hat{P}]={\rm{i}}\lambda .$ We write Hamiltonian (A.1) in an alternative form

$$\begin{eqnarray}&&H(t)=\displaystyle \frac{1}{2\lambda }{\rm{\hslash }}\omega ({\hat{Q}}^{2}+{\hat{P}}^{2})+2A\mathrm{cos}(\hat{Q}+{\omega }_{d}t),\end{eqnarray} \tag{ A.2 }$$

Now, we employ an unitary operator $U={{\rm{e}}}^{{\rm{i}}{\rm{\hslash }}\displaystyle \frac{{\omega }_{d}}{n}{a}^{\dagger }{at}}$ to transform Hamiltonian (A.2) into a rotating frame with frequency ${\omega }_{d}/n$

$$\begin{eqnarray}{H}_{{\rm{RF}}} & = & {UH}(t){U}^{\dagger }-{\rm{i}}{U}{\dot{U}}^{\dagger }\\ & = & \displaystyle \frac{1}{2\lambda }{\rm{\hslash }}\delta \omega ({\hat{Q}}^{2}+{\hat{P}}^{2})+2{AU}\mathrm{cos}(\hat{Q}+{\omega }_{d}t){U}^{\dagger }\\ & = & \displaystyle \frac{1}{2\lambda }{\rm{\hslash }}\delta \omega ({\hat{Q}}^{2}+{\hat{P}}^{2})+({{A}{\rm{e}}}^{{\rm{i}}{\omega }_{d}t}{{U}{\rm{e}}}^{{\rm{i}}\hat{Q}}{U}^{\dagger }+{\rm{h.c.}})\\ & = & \displaystyle \frac{1}{2\lambda }{\rm{\hslash }}\delta \omega ({\hat{Q}}^{2}+{\hat{P}}^{2})+\{{{A}{\rm{e}}}^{{\rm{i}}{\omega }_{d}t}{{\rm{e}}}^{{\rm{i}}[\hat{Q}\mathrm{cos}({\omega }_{d}t/n)+\hat{P}\mathrm{sin}({\omega }_{d}t/n)]}+{\rm{h.c.}}\}\\ & \equiv & \displaystyle \frac{1}{2\lambda }{\rm{\hslash }}\delta \omega ({\hat{Q}}^{2}+{\hat{P}}^{2})+[{{A}{\rm{e}}}^{{\rm{i}}{\omega }_{d}t}M(\hat{Q},\hat{P})+{\rm{h.c.}}].\end{eqnarray} \tag{ A.3 }$$

Here, we define $M(\hat{Q},\hat{P})\equiv {{\rm{e}}}^{{\rm{i}}[\hat{Q}\mathrm{cos}({\omega }_{d}t/n)+\hat{P}\mathrm{sin}({\omega }_{d}t/n)]}$ and the detuning $\delta \omega \equiv {\omega }_{0}-{\omega }_{d}/n$. To calculate the matrix element of $M(\hat{Q},\hat{P})$, we define the displacement operator $D(\alpha ,{\alpha }^{*})$ by

$$\begin{eqnarray}&&D(\alpha ,{\alpha }^{*})\equiv \mathrm{exp}(\alpha {a}^{\dagger }-{\alpha }^{*}a)=\mathrm{exp}(\mathrm{Re}[\alpha ]({a}^{\dagger }-a)+i\ \mathrm{Im}[\alpha ]({a}^{\dagger }+a)).\end{eqnarray} \tag{ A.4 }$$

Since the operator $M(\hat{Q},\hat{P})$ can be written as

$$\begin{eqnarray}M(\hat{Q},\hat{P}) & \equiv & \mathrm{exp}[{\rm{i}}\ (\hat{Q}\mathrm{cos}({\omega }_{d}t/n)+\hat{P}\mathrm{sin}({\omega }_{d}t/n))]\\ & = & \mathrm{exp}\left[-\sqrt{\displaystyle \frac{\lambda }{2}}\mathrm{sin}({\omega }_{d}t/n)({a}^{\dagger }-a)+{\rm{i}}\sqrt{\displaystyle \frac{\lambda }{2}}\mathrm{cos}({\omega }_{d}t/n)({a}^{\dagger }+a)\right],\end{eqnarray} \tag{ A.5 }$$

we get the relationship between the parameter α of $D(\alpha ,{\alpha }^{*})$ and parameters of $M(\hat{Q},\hat{P})$

$$\begin{eqnarray}&&\alpha \equiv -\sqrt{\displaystyle \frac{\lambda }{2}}\mathrm{sin}({\omega }_{d}t/n)+{\rm{i}}\sqrt{\displaystyle \frac{\lambda }{2}}\mathrm{cos}({\omega }_{d}t/n)=\sqrt{\displaystyle \frac{\lambda }{2}}{{\rm{e}}}^{{\rm{i}}(\varphi +\pi /2)},\end{eqnarray} \tag{ A.6 }$$

with $\varphi ={\omega }_{d}t/n$. We further define the following notations

$$\begin{eqnarray}\mathrm{Coherent}\ \ \mathrm{state}:| \alpha \rangle & \equiv & {{\rm{e}}}^{-\displaystyle \frac{1}{2}{| \alpha | }^{2}}\displaystyle \sum _{k=0}^{\infty }\displaystyle \frac{{\alpha }^{k}}{\sqrt{k!}}| k\rangle ,\\ \langle \beta | \alpha \rangle & = & {{\rm{e}}}^{\alpha {\beta }^{*}-({| \alpha | }^{2}+{| \beta | }^{2})/2},\\ \mathrm{Displaced}\ \ \mathrm{Fock}\ \ \mathrm{state}:\ \ | \alpha ,k\rangle & \equiv & D(\alpha ,{\alpha }^{*})| k\rangle ,\\ | 0,k\rangle & = & | k\rangle .\end{eqnarray} \tag{ A.7 }$$

According to equation (3.11) in [55], we have

$$\begin{eqnarray}&&\langle \beta ,l| \alpha ,k\rangle =\langle \beta | \alpha \rangle \sqrt{\displaystyle \frac{l!}{k!}}{({\beta }^{*}-{\alpha }^{*})}^{k-l}{L}_{l}^{k-l}\{{| \beta -\alpha | }^{2}\}.\end{eqnarray} \tag{ A.8 }$$

Here, ${L}_{l}^{k-l}(\bullet )$ is the Laguerre polynomials. Let $\beta =0$, we have the exact form of matrix element of displacement operator $D(\alpha ,{\alpha }^{*})$

$$\begin{eqnarray}&&\langle l| \alpha ,k\rangle \equiv \langle l| D(\alpha ,{\alpha }^{*})| k\rangle ={{\rm{e}}}^{-{| \alpha | }^{2}/2+{\rm{i}}\pi (k-l)}{({\alpha }^{*})}^{k-l}\sqrt{\displaystyle \frac{l!}{k!}}{L}_{l}^{k-l}({| \alpha | }^{2}).\end{eqnarray} \tag{ A.9 }$$

Using the relationship (A.6) we get the explicit form of matrix elements of $M(\hat{Q},\hat{P})$

$$\begin{eqnarray}&&\langle l| M(\hat{Q},\hat{P})| k\rangle ={{\rm{e}}}^{-\lambda /4+{\rm{i}}(k-l)(\pi /2-{\omega }_{d}t/n)}\sqrt{\displaystyle \frac{l!}{k!}}\ {\left(\displaystyle \frac{\lambda }{2}\right)}^{\displaystyle \frac{k-l}{2}}{L}_{l}^{k-l}(\lambda /2).\end{eqnarray} \tag{ A.10 }$$

Thus, quantum Hamiltonian (A.3) is

$$\begin{eqnarray}{H}_{{\bf{RF}}} & = & \displaystyle \frac{1}{2\lambda }{\rm{\hslash }}\delta \omega ({\hat{Q}}^{2}+{\hat{P}}^{2})+[{{A}{\rm{e}}}^{{\rm{i}}{\omega }_{d}t}M(\hat{Q},\hat{P})+{\rm{h.c.}}]\\ & = & {\rm{\hslash }}\delta \omega \left({a}^{\dagger }a+\displaystyle \frac{1}{2}\right)+[{{A}{\rm{e}}}^{{\rm{i}}{\omega }_{d}t}M(\hat{Q},\hat{P})+{\rm{h.c.}}]\\ & = & {\rm{\hslash }}\delta \omega \left({a}^{\dagger }a+\displaystyle \frac{1}{2}\right)+A\left[\displaystyle \sum _{k,l}\langle l| M(\hat{Q},\hat{P})| k\rangle {{\rm{e}}}^{{\rm{i}}{\omega }_{d}t}| l\rangle \langle k| +{\rm{h.c.}}\right]\\ & = & {\rm{\hslash }}\delta \omega \left({a}^{\dagger }a+\displaystyle \frac{1}{2}\right)+A\left[\displaystyle \sum _{k,l}{{\rm{e}}}^{{\rm{i}}{\omega }_{d}t}{{\rm{e}}}^{-\lambda /4+{\rm{i}}(k-l)(\pi /2-{\omega }_{d}t/n)}\right.\\ & & \times \left.\sqrt{\displaystyle \frac{l!}{k!}}\ {\left(\displaystyle \frac{\lambda }{2}\right)}^{\displaystyle \frac{k-l}{2}}{L}_{l}^{k-l}(\lambda /2)| l\rangle \langle k| +{\rm{h.c.}}\right].\end{eqnarray} \tag{ A.11 }$$

Under RWA, we drop the fast oscillating terms ($k-l\ne n$) and get RWA Hamiltonian ($k-l=n$)

$$\begin{eqnarray} & & {H}_{{\rm{RWA}}}\\ & = & {\rm{\hslash }}\delta \omega \left({a}^{\dagger }a+\displaystyle \frac{1}{2}\right)+A\left[\displaystyle \sum _{l}{{\rm{e}}}^{-\lambda /4+{\rm{i}}{n}\pi /2}\right.\\ & & \times \left.\sqrt{\displaystyle \frac{l!}{(l+n)!}}\ {\left(\displaystyle \frac{\lambda }{2}\right)}^{\displaystyle \frac{n}{2}}{L}_{l}^{n}(\lambda /2)\ | l\rangle \langle l+n| +{\rm{h.c.}}\right]\\ & = & {\rm{\hslash }}\delta \omega \left({a}^{\dagger }a+\displaystyle \frac{1}{2}\right)+A\left[{{\rm{e}}}^{-\lambda /4+{\rm{i}}{n}\pi /2}{\left(\displaystyle \frac{\lambda }{2}\right)}^{\displaystyle \frac{n}{2}}\right.\\ & & \times \left.\displaystyle \sum _{l}| l\rangle \langle l+n| \sqrt{\displaystyle \frac{l!}{(l+n)!}}\ {L}_{l}^{n}(\lambda /2)+{\rm{h.c.}}\right]\\ & = & {\rm{\hslash }}\delta \omega \left({a}^{\dagger }a+\displaystyle \frac{1}{2}\right)+A\left[{{\rm{e}}}^{-\lambda /4-{\rm{i}}{n}\pi /2}{\left(\displaystyle \frac{\lambda }{2}\right)}^{-\displaystyle \frac{n}{2}}\right.\\ & & \times \left.\displaystyle \sum _{l}| l\rangle \langle l+n| \sqrt{\displaystyle \frac{(l+n)!}{l!}}\ {L}_{l+n}^{-n}(\lambda /2)+{\rm{h.c.}}\right]\\ & = & {\rm{\hslash }}\delta \omega ({a}^{\dagger }a+\displaystyle \frac{1}{2})+A\left[{{\rm{e}}}^{-\lambda /4-{\rm{i}}{n}\pi /2}{\left(\displaystyle \frac{\lambda }{2}\right)}^{-\displaystyle \frac{n}{2}}\right.\\ & & \times \left.\displaystyle \sum _{l}| l\rangle \langle l| {a}^{n}\ {L}_{{a}^{\dagger }a}^{-n}(\lambda /2)+{\rm{h.c.}}\right]\\ & = & {\rm{\hslash }}\delta \omega ({a}^{\dagger }a+\displaystyle \frac{1}{2})+A\left[{{\rm{e}}}^{-\lambda /4-{\rm{i}}{n}\pi /2}{\left(\displaystyle \frac{\lambda }{2}\right)}^{-\displaystyle \frac{n}{2}}{a}^{n}\ {L}_{{a}^{\dagger }a}^{-n}(\lambda /2)+{\rm{h.c.}}\right].\end{eqnarray} \tag{ A.12 }$$

Here we have used the relationship [55]${L}_{l}^{n}(x)/{L}_{l+n}^{-n}(x)={(-x)}^{-n}(l+n)!/l!$ for $x\gt 0$. We now scale the RWA Hamiltonian by ${\rm{\hslash }}\omega /\lambda$ and get the dimensionless Hamiltonian $\hat{g}$

$$\begin{eqnarray}\hat{g} & \equiv & \displaystyle \frac{\lambda }{{\rm{\hslash }}\omega }{H}_{{\rm{RWA}}}\\ & = & \lambda \epsilon \left({a}^{\dagger }a+\displaystyle \frac{1}{2}\right)+\mu \left[{{\rm{e}}}^{-\lambda /4-{\rm{i}}{n}\pi /2}{\left(\displaystyle \frac{\lambda }{2}\right)}^{-\displaystyle \frac{n}{2}}{a}^{n}\ {L}_{{a}^{\dagger }a}^{-n}(\lambda /2)+{\rm{h.c.}}\right],\end{eqnarray} \tag{ A.13 }$$

$\epsilon =\delta \omega /\omega$ and $\mu =\lambda A/({\rm{\hslash }}\omega )$ are the dimensionless detuning and driving strength, respectively.

Using the following asymptotic form of Laguerre polynomials [56, 57]

$$\begin{eqnarray}&&\underset{k\to \infty }{\mathrm{lim}}{L}_{k}^{\alpha }(x/k)={k}^{\alpha }{{\rm{e}}}^{\displaystyle \frac{x}{2k}}{x}^{-\alpha /2}{J}_{\alpha }(2\sqrt{x}),\end{eqnarray} \tag{ A.14 }$$

we have the following relationship in the limit of $k,l\gg | k-l|$ for a fixed k − l

$$\begin{eqnarray}{L}_{k}^{k-l}(\lambda /2) & \approx & {{\rm{e}}}^{\lambda /4}{k}^{k-l}{(k\lambda /2)}^{-(k-l)/2}{J}_{k-l}(2\sqrt{k\lambda /2})\\ & = & {{\rm{e}}}^{\lambda /4}{\left(\displaystyle \frac{\lambda }{2k}\right)}^{-\displaystyle \frac{1}{2}(k-l)}{J}_{k-l}(\sqrt{2k\lambda }).\end{eqnarray} \tag{ A.15 }$$

Thus, in the semiclassical limit, i.e., $k,l\to \infty$ and fixed k − l, equation (A.10) goes to the following

$$\begin{eqnarray}\langle l| M(\hat{Q},\hat{P})| k\rangle & = & {{\rm{e}}}^{-\lambda /4+{\rm{i}}(k-l)(\pi /2-\varphi )}\sqrt{\displaystyle \frac{l!}{k!}}\ {\left(\displaystyle \frac{\lambda }{2}\right)}^{\displaystyle \frac{k-l}{2}}{L}_{l}^{k-l}(\lambda /2)\\ & \approx & {{\rm{e}}}^{{\rm{i}}(k-l)(\pi /2-\varphi )}{k}^{(k-l)/2}\sqrt{\displaystyle \frac{l!}{k!}}{J}_{k-l}(\sqrt{2k\lambda })\\ & \approx & {{\rm{e}}}^{{\rm{i}}(k-l)(\pi /2-\varphi )}{J}_{k-l}(\sqrt{2k\lambda }).\end{eqnarray} \tag{ A.16 }$$

Here, we have used the limit relationship $\sqrt{\displaystyle \frac{l!}{k!}}\ {k}^{(k-l)/2}\to 1$. Therefore, we have the RWA Hamiltonian (A.13) in the Fock representation $\hat{g}={\displaystyle \sum }_{k=0,l=0}^{\infty }f(k,l)| k\rangle \langle l|$ with

$$\begin{eqnarray}&&f(k,l)\approx \lambda \epsilon (k+\displaystyle \frac{1}{2}){\delta }_{k,l}+\mu \left[{{\rm{e}}}^{{\rm{i}}{n}\displaystyle \frac{\pi }{2}}{J}_{k-l}(\sqrt{\lambda (k+l+1)}){\delta }_{l-k,n}+{\rm{h.c.}}\right].\end{eqnarray} \tag{ A.17 }$$

We define the radial and angular operators $\hat{r}$ and $\hat{\theta }$ by $a={{\rm{e}}}^{-{\rm{i}}\theta }\hat{r}/\sqrt{2\lambda }$ and ${a}^{\dagger }=\hat{r}{{\rm{e}}}^{{\rm{i}}\theta }/\sqrt{2\lambda }$. In the Fock representation, the operator ${{\rm{e}}}^{{\rm{i}}\hat{\theta }}$ is defined by

$$\begin{eqnarray}&&{{\rm{e}}}^{{\rm{i}}\hat{\theta }}=\displaystyle \sum _{k=0}^{\infty }| k\rangle \langle k+1| ,\ \ \ \ \mathrm{and}\ \ \ \ \ {{\rm{e}}}^{-{\rm{i}}\hat{\theta }}=\displaystyle \sum _{k=0}^{\infty }| k+1\rangle \langle k| .\end{eqnarray} \tag{ A.18 }$$

Using the above relationships, we have the following Hamiltonian in the semiclassical limit $\lambda \to 0$

$$\begin{eqnarray}&&\hat{g}\to g=\displaystyle \frac{1}{2}\epsilon {r}^{2}+2\mu {J}_{n}(r)\mathrm{cos}(n\theta -\displaystyle \frac{n\pi }{2}).\end{eqnarray} \tag{ A.19 }$$

In this section, we give a detailed description about the quantum tunnelling process in phase space and the analytical behavior of 'momentum' θ in the complex plane. We also calculate the asymmetry factor δ and show its linear relationship with $1/\lambda$ for different bands. To be convenient, we define a new variable

$$\begin{eqnarray*}&&\hat{\xi }\equiv {\hat{r}}^{2}/2.\end{eqnarray*}$$

The semiclassical Hamiltonian (A.19) can be rewritten as

$$\begin{eqnarray*}&&g=\epsilon \xi +2\mu {J}_{n}(\sqrt{2\xi })\mathrm{cos}\left(n\theta -\displaystyle \frac{n\pi }{2}\right)\end{eqnarray*}$$

by new variables ξ and θ, which define the '$\xi -\theta$' phase space for our WKB calculation. For a fixed g, the general solutions of classical orbits are

$$\begin{eqnarray}&&{\theta }_{\pm }(\xi ,g)=\displaystyle \frac{\pi }{2}\pm \displaystyle \frac{1}{n}\mathrm{arccos}\left[\displaystyle \frac{g-\epsilon \xi }{2\mu {J}_{n}(\sqrt{2\xi })}\right]+k\displaystyle \frac{2\pi }{n},\end{eqnarray} \tag{ B.1 }$$

where $k=0,1,2\ ...,\ \mathrm{and}\ n-1$ represent the n branches of solutions. Here, we choose the parameters $\epsilon =0$ and $\mu =-1$. In figure B.1 , we show three classical orbits for a fixed $g\lt 0$. The two classical orbits in the first loop are indicated by red closed curves, which correspond to the following solutions

$$\begin{eqnarray}&&{\theta }_{\pm }(\xi ,g)=-\displaystyle \frac{\pi }{n}\pm \displaystyle \frac{1}{n}\mathrm{arccos}\left[\displaystyle \frac{g-\epsilon \xi }{2\mu {J}_{n}(\sqrt{2\xi })}\right],\end{eqnarray} \tag{ B.2a }$$

and

$$\begin{eqnarray}&&{\theta }_{\pm }(\xi ,g)=\displaystyle \frac{\pi }{n}\pm \displaystyle \frac{1}{n}\mathrm{arccos}\left[\displaystyle \frac{g-\epsilon \xi }{2\mu {J}_{n}(\sqrt{2\xi })}\right].\end{eqnarray} \tag{ B.2b }$$

The classical orbit in the second loop is indicated by yellow closed curve, which corresponds to the following solution

$$\begin{eqnarray}&&{\theta }_{\pm }(\xi ,g)=\pm \displaystyle \frac{1}{n}\left(\pi -\mathrm{arccos}\left[\displaystyle \frac{g-\epsilon \xi }{2\mu {J}_{n}(\sqrt{2\xi })}\right]\right).\end{eqnarray} \tag{ B.3 }$$

In the regime of $| (g-\epsilon \xi )/[2\mu {J}_{n}(\sqrt{2\xi })]| \lt 1$, two real solutions ${\theta }_{\pm }(\xi ,g)$ together represent one closed classical orbit $\theta (\xi ,g)$. In figure B.1 (left), the boundaries of classical motions are indicated by the white dashed lines, i.e., ${\xi }_{1}$, ${\xi }_{2}$ and ${\xi }_{3}$. Beyond the classical boundaries, the value of $\theta (\xi ,g)$ has imaginary part. In figure B.1 (right), we show the analytical structures of solutions ${\theta }_{\pm }(\xi ,g)$ in the complex plane. The closed curves on the real axis of θ represent classical orbits (we deviate the orbits slightly from the real axis to illustrate the shapes of orbits). There are n identical orbital branches with only a $2\pi /n$-shift of Re[θ] for each type of solution.

Zoom In Zoom Out Reset image size

In the quantum regime, the classical orbits can tunnel into each other through the classical forbidden region. In figure B.1 (left), we show the quantum tunnelling process of the two states in the first loop in phase space. The corresponding behavior of Im[θ] is depicted in figure B.1 (right). Starting from the classical boundary ${\xi }_{2}$ to the zero point of Bessel function ${\xi }^{(0)}$, the imaginary part Im[θ] increases from zero to infinite, where it jumps to another branch of solution. Then it goes back from infinite to zero as ξ changes from ${\xi }^{(0)}$ to another classical boundary ${\xi }_{3}$. After that, Im[θ] increases again from zero to infinite as ξ goes from ${\xi }_{2}$ to ${\xi }^{(0)}$, where it jumps again to another branch of solution. Finally, Im[θ] decreases from infinite to zero as ξ changes from ${\xi }^{(0)}$ to the classical boundary ${\xi }_{3}$. As we have discussed in the main text, the amplitude of quantum tunnelling rate J_l is given by the integral of the imaginary part of 'momentum' θ in the classical forbidden region ${\xi }_{2}\lt \xi \lt {\xi }_{3}$

$$\begin{eqnarray}&&| {J}_{l}| =\displaystyle \frac{\lambda }{2\pi }{{\left(\displaystyle \frac{{\rm{d}}{S}}{{\rm{d}}{g}}\right)}^{-1}| }_{g={g}_{l}}\mathrm{exp}\left(-\displaystyle \frac{2}{\lambda }{\displaystyle \int }_{{\xi }_{2}}^{{\xi }_{3}}\mathrm{Im}[\theta ]{\rm{d}}\xi \right).\end{eqnarray} \tag{ B.4 }$$

The tunnelling process can also happen through lower boundary ${\xi }_{1}$ as indicated by the white arrows in figure B.1 (left). However, the lower path is much longer than the upper path. Thus, the contribution to $| {J}_{l}|$ from the lower path is exponentially smaller than the contribution from upper path.

The jumping processes between different branches of solutions give additional phases to the quantum tunnelling rate J_l, which makes it a complex number ${J}_{l}=| {J}_{l}| {{\rm{e}}}^{-{\rm{i}}{\delta }_{l}\tau }$. As we have discussed in the main text, the connecting condition by jumping is given by the phase translation operator ${\hat{T}}_{\tau }={{\rm{e}}}^{-{\rm{i}}\tau {\hat{a}}^{\dagger }\hat{a}}$. Since ${\xi }^{(0)}\gg 1$, we can expand operator ${\hat{T}}_{\tau }$ by [1]${\hat{a}}^{\dagger }\hat{a}\approx {\xi }^{(0)}/\lambda +{\rm{i}}\partial /\partial \theta$. As a result, the connecting condition is ${\hat{T}}_{\tau }{\phi }_{l}(\theta )\approx {{\rm{e}}}^{-{\rm{i}}{\xi }^{(0)}\tau /\lambda }{\phi }_{l}(\theta +\tau )$. Thus we get the symmetry factor

$$\begin{eqnarray}&&{\delta }_{l}={\delta }_{l}^{0}+{\xi }^{(0)}/\lambda ,\end{eqnarray} \tag{ B.5 }$$

where ${\delta }_{l}^{0}$ is the residual asymmetry beyond WKB calculation. In figure B.2 (a), we compare the above linear relationships between ${\delta }_{l}$ and $1/\lambda$ for different bands to our numerical simulations. In figure B.2 (b) and figure B.2 (c), we expand the asymmetry factor to the whole field of real number ${\mathbb{R}}$ and plot it as function of $1/\lambda$ for different bands. The bands in figure B.2 (b) are all in the first loop. We see that, since the states in the first loop tunnel through the upper boundary, they all have the same slope given by ${\xi }^{(0)}$, which is the second zero point of Bessel function ${J}_{n}(\sqrt{2\xi })$. Here, we consider ${\xi }^{(0)}=0$ is the first zero point of Bessel function ${J}_{n}(\sqrt{2\xi })$ for $n\ne 0$.

Zoom In Zoom Out Reset image size

In figure B.2 (c), we show the linear relationships between ${\delta }_{l}$ and $1/\lambda$ for the bottom bands in different loops. We see their slopes are different. The reason is that the bands in different loops tunnel though different paths with different jumping points ${\xi }^{(0)}$. Like the states in the first loop, the states in other loops can tunnel through both the upper boundary and lower boundary. However, we have checked the integral $\int \mathrm{Im}[\theta ]{\rm{d}}\xi$ of the upper path is always larger than that of the lower path. Therefore, the contribution to the tunnelling rate from the upper path is exponentially smaller than the contribution from the lower path. Therefore, the slope of all the bands in the l-th ($l\gt 1$) loop is given by the l-th zero point ${\xi }_{l}^{(0)}$ of Bessel function ${J}_{n}(\sqrt{2\xi })$. In the flowing table, we compare the slopes extracted form numerical simulation to our theoretical calculation.

Band index	${\xi }^{(0)}$ (Numerical)	${\xi }^{(0)}$ (Theory)	Relative errors
I-1	642.241	651.545	−0.014
II-1	629.514	651.545	−0.034
III-1	860.600	844.308	0.019
IV-1	1021.829	1032.972	−0.011
V-1	1186.088	1225.435	−0.032
VI-1	1427.519	1424.378	0.002
VII-1	1662.219	1631.067	0.019
VIII-1	1820.811	1846.185	−0.014
IX-1	2056.534	2070.142	−0.007

In this section, we calculate the squeezing factor u of localized states near the stable points of phase space lattice. First, we determine all the extrema $({r}_{e},{\theta }_{e})$ including stable points and unstable saddle points by the derivatives of quasienergy (A.19) along both angular direction and radial direction

$$\begin{eqnarray}&&{\displaystyle \frac{\partial g}{\partial \theta }|}_{({r}_{e},{\theta }_{e})}=-2n\mu {J}_{n}({r}_{e})\mathrm{sin}\left(n{\theta }_{e}-\displaystyle \frac{n\pi }{2}\right)=0,\end{eqnarray} \tag{ C.1a }$$

$$\begin{eqnarray}&&{\displaystyle \frac{\partial g}{\partial r}|}_{({r}_{e},{\theta }_{e})}=\epsilon {r}_{e}+\mu ({J}_{n-1}({r}_{e})-{J}_{n+1}({r}_{e}))\mathrm{cos}\left(n{\theta }_{e}-\displaystyle \frac{n\pi }{2}\right)=0.\end{eqnarray} \tag{ C.1b }$$

The angular extrema can be obtained from equation (C.1a), that is, ${\theta }_{e}=l\tau /2+n\pi /2$ with $l=0,\pm 1,\pm 2,\ldots ,\pm (n-1),n$, where $\tau =2\pi /n$ is defined as lattice constant of phase space lattice. To get the radial extrema, we need to solve the equation (C.1b). The stability of these extrema $({r}_{e},{\theta }_{e})$ is determined by the second derivatives of g. If $({\partial }^{2}g/\partial {\theta }^{2})\times ({\partial }^{2}g/\partial {r}^{2}){| }_{r={r}_{e},\theta ={\theta }_{e}}\gt 0$, the extrema are stable, otherwise unstable. The second derivatives to angle θ and radius r are

$$\begin{eqnarray}&&{\displaystyle \frac{{\partial }^{2}g}{\partial {\theta }^{2}}|}_{({r}_{e},{\theta }_{e})}=-2{n}^{2}\mu {J}_{n}({r}_{e})\mathrm{cos}\left(n{\theta }_{e}-\displaystyle \frac{n\pi }{2}\right),\end{eqnarray} \tag{ C.2a }$$

$$\begin{eqnarray}&&{\displaystyle \frac{{\partial }^{2}g}{\partial {r}^{2}}|}_{({r}_{e},{\theta }_{e})}=\epsilon +\displaystyle \frac{1}{2}\mu ({J}_{n-2}({r}_{e})+{J}_{n+2}({r}_{e})-2{J}_{n}({r}_{e}))\mathrm{cos}\left(n{\theta }_{e}-\displaystyle \frac{n\pi }{2}\right).\end{eqnarray} \tag{ C.2b }$$

Below, we label the stable points (maxima and minima) and unstable saddle points by $({r}_{m},{\theta }_{m})$ and $({r}_{s},{\theta }_{s})$, respectively. We expand the quasienergy g near the stable points $({r}_{m},{\theta }_{m})$ to the second order

$$\begin{eqnarray}{g}_{{\rm{local}}} & \approx & g({r}_{m},{\theta }_{m})+\displaystyle \frac{1}{2}{\displaystyle \frac{{\partial }^{2}g}{{r}_{m}^{2}\partial {\theta }^{2}}|}_{({r}_{m},{\theta }_{m})}{({r}_{m}\theta -{r}_{m}{\theta }_{m})}^{2}+\displaystyle \frac{1}{2}{\displaystyle \frac{{\partial }^{2}g}{\partial {r}^{2}}|}_{({r}_{m},{\theta }_{m})}{(r-{r}_{m})}^{2}\\ & = & g({r}_{m},{\theta }_{m})+\displaystyle \frac{{\tilde{p}}^{2}}{2{m}_{e}}+\displaystyle \frac{1}{2}{m}_{e}{\omega }_{e}^{2}{\tilde{x}}^{2}.\end{eqnarray} \tag{ C.3 }$$

Here, we have defined coordinate $\tilde{x}=r-{r}_{m}$ and momentum $\tilde{p}={r}_{m}(\theta -{\theta }_{m})$ near the stable point. The effective mass m_e and effective frequency ${\omega }_{e}$ are given by

$$\begin{eqnarray*}&&{m}_{e}={r}_{m}^{2}{{\left(\displaystyle \frac{{\partial }^{2}g}{\partial {\theta }^{2}}\right)}^{-1}|}_{({r}_{m},{\theta }_{m})}\ \ \mathrm{and}\ \ {\omega }_{e}=\sqrt{{m}_{e}^{-1}\displaystyle \frac{{\partial }^{2}g}{\partial {r}^{2}}{| }_{({r}_{m},{\theta }_{m})}}\end{eqnarray*}$$

respectively.

Now, we define the displacement operator ${\hat{D}}_{\alpha }=\mathrm{exp}(\alpha {a}^{\dagger }-{\alpha }^{*}a)$ and the squeezing operator ${\hat{S}}_{\xi }=\mathrm{exp}\left[\displaystyle \frac{1}{2}({\xi }^{*}{a}^{2}-\xi {a}^{\dagger 2})\right]$, which have the following properties

$$\begin{eqnarray*}&&{\hat{D}}_{\alpha }^{\dagger }a{\hat{D}}_{\alpha }=a+\alpha ,\ \ {\hat{S}}_{\xi }^{\dagger }a{\hat{S}}_{\xi }={va}+{{ua}}^{\dagger }\end{eqnarray*}$$

with $\xi ={{r}{\rm{e}}}^{{\rm{i}}\theta }$. The squeezing parameters are given by $v=\mathrm{cosh}r,\ u=-{{\rm{e}}}^{{\rm{i}}\theta }\mathrm{sinh}r$ . We transform the original $\hat{g}$ to localized Hamiltonian ${\hat{g}}_{{\rm{local}}}$ at the stable point $({r}_{m},{\theta }_{m})$ by three operators, i.e.,

$$\begin{eqnarray*}&&{\hat{g}}_{{\rm{local}}}={\hat{S}}_{\xi }{\hat{D}}_{\alpha }{\hat{T}}_{{\theta }_{m}}\hat{g}{\hat{T}}_{{\theta }_{m}}^{\dagger }{\hat{D}}_{\alpha }^{\dagger }{\hat{S}}_{\xi }^{\dagger }.\end{eqnarray*}$$

Here, we first change the orientation using phase space rotation operator ${\hat{T}}_{{\theta }_{m}}={{\rm{e}}}^{-{\rm{i}}{\theta }_{m}{\hat{a}}^{\dagger }\hat{a}}$. Then we move the Hamiltonian to the position of stable point using displacement operator ${\hat{D}}_{\alpha }={{\rm{e}}}^{\alpha {a}^{\dagger }-{\alpha }^{*}a}$. Finally, we squeeze the Hamiltonian to fit the stable point using squeezing operator ${\hat{S}}_{\xi }={{\rm{e}}}^{[{\xi }^{*}{a}^{2}-\xi {({a}^{\dagger })}^{2}]/2}$. By choosing

$$\begin{eqnarray}&&\alpha =-\displaystyle \frac{{r}_{m}}{\sqrt{2\lambda }},v=\displaystyle \frac{1}{2}\left(\sqrt{{m}_{e}{\omega }_{e}}+\displaystyle \frac{1}{\sqrt{{m}_{e}{\omega }_{e}}}\right),u=\displaystyle \frac{1}{2}\left(\sqrt{{m}_{e}{\omega }_{e}}-\displaystyle \frac{1}{\sqrt{{m}_{e}{\omega }_{e}}}\right),\end{eqnarray} \tag{ C.4 }$$

we get the localized Hamiltonian as follows

$$\begin{eqnarray}{\hat{g}}_{{\rm{local}}} & = & {\hat{S}}_{\xi }{\hat{D}}_{\alpha }{\hat{T}}_{{\theta }_{m}}\hat{g}{\hat{T}}_{{\theta }_{m}}^{\dagger }{\hat{D}}_{\alpha }^{\dagger }{\hat{S}}_{\xi }^{\dagger }\\ & = & \lambda {\omega }_{e}\left({a}^{\dagger }a+\displaystyle \frac{1}{2}\right)+g({r}_{m},{\theta }_{m})+o({\lambda }^{2}).\end{eqnarray} \tag{ C.5 }$$

We investigate the quantum dynamics near the bottom of a stable state. The dissipative dynamics is modified by squeezing and can be described by an effective temperature ${\bar{n}}_{e}$. The original master equation is

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial \tau }=-\displaystyle \frac{{\rm{i}}}{\lambda }[\hat{g},\rho ]+\kappa (1+\bar{n}){ \mathcal D }[a]\rho +\kappa \bar{n}{ \mathcal D }[{a}^{\dagger }]\rho .\end{eqnarray} \tag{ D.1 }$$

The Lindblad superoperator is defined through ${ \mathcal D }[A]\rho \equiv A\rho {A}^{\dagger }-({A}^{\dagger }A\rho +\rho {A}^{\dagger }A)/2$, $\bar{n}={({{\rm{e}}}^{{\rm{\hslash }}{\omega }_{0}/{k}_{B}T}-1)}^{-1}$ is the Bose distribution and κ is the dimensionless damping scaled ω. By performing a transformation on the density operator

$$\begin{eqnarray*}&&\tilde{\rho }={\hat{S}}_{\xi }{\hat{D}}_{\alpha }{\hat{T}}_{{\theta }_{m}}\rho {\hat{T}}_{{\theta }_{m}}^{\dagger }{\hat{D}}_{\alpha }^{\dagger }{\hat{S}}_{\xi }^{\dagger },\end{eqnarray*}$$

we transform the master equation (D.1) into the following form [58]

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}\tilde{\rho }}{{\rm{d}}\tau } & = & -{\rm{i}}[{\hat{g}}_{{\rm{local}}},\tilde{\rho }]+\displaystyle \frac{\kappa }{2}\{(1+{\bar{n}}_{e}){ \mathcal D }[a]\tilde{\rho }+{\bar{n}}_{e}{ \mathcal D }[{a}^{\dagger }]\tilde{\rho }\}\\ & & +\displaystyle \frac{\kappa }{2}M(2{a}^{\dagger }\tilde{\rho }{a}^{\dagger }-{{a}^{\dagger }}^{2}\tilde{\rho }-\tilde{\rho }{{a}^{\dagger }}^{2})+\displaystyle \frac{\kappa }{2}{M}^{*}(2a\tilde{\rho }a-{a}^{2}\tilde{\rho }-\tilde{\rho }{a}^{2}).\end{eqnarray} \tag{ D.2 }$$

Here, parameter $M={{uv}}^{*}(2\bar{n}+1)$ is the squeezing number. The effective Bose distribution is given by

$$\begin{eqnarray}&&{\bar{n}}_{e}=\bar{n}{| v| }^{2}+(1+\bar{n}){| u| }^{2}={| u| }^{2}+\bar{n}(2{| u| }^{2}+1).\end{eqnarray} \tag{ D.3 }$$

Near the bottom of stable points, we can make the harmonic approximation. The squeezing number $M={{uv}}^{*}(2\bar{n}+1)$ has no contribution to the stationary distribution. The ration of probability over adjoint levels thus is given approximately by [58]${\bar{n}}_{e}/(1+{\bar{n}}_{e})$.