Multiple-period Floquet states and time-translation symmetry breaking in quantum oscillators

TL;DR

This paper investigates how small quantum oscillators can exhibit broken time-translation symmetry through period tripling, demonstrating robustness over long timescales and under weak decoherence.

Contribution

It introduces the concept of multiple-period Floquet states in small quantum systems and analyzes their stability and robustness against decoherence.

Findings

Period tripling persists for exponentially long times under moderate driving.

Weak decoherence further stabilizes the period tripling.

The study bridges behavior between large and small quantum systems.

Abstract

We study the breaking of the discrete time-translation symmetry in small periodically driven quantum systems. Such systems are intermediate between large closed systems and small dissipative systems, which both display the symmetry breaking, but have qualitatively different dynamics. As a nontrivial example we consider period tripling in a quantum nonlinear oscillator. We show that, for moderately strong driving, the period tripling is robust on an exponentially long time scale, which is further extended by an even weak decoherence.

Equations

$$H=H_0+H_F,H_0=\frac{1}{2}{p}^2+\frac{1}{2}{\omega }_0^2{q}^2+\frac{1}{4}\gamma q^4,$$

$$H_{\mathrm{RWA}}=(K{t}_F{)}^{-1}{\int }_0^{K{t}_F}dt{H}_K(t).$$

$$T_{t_F}\psi (t)=e^{-iE{t}_F/\hslash }U(t+t_F)\varphi (t)=e^{-iE{t}_F/\hslash }{N}_K\psi (t).$$

$$N_K=\exp (-2\pi i{a}^{†}a/K),\left[N_K,H_{\mathrm{RWA}}\right]=0.$$

$$N_K{\varphi }^{(k)}=\exp (-2\pi ik/K){\varphi }^{(k)},0\le k\le K-1.$$

$${\epsilon }^{(k)}=(E^{(k)}+\hslash {\omega }_{F}k/K)\mathrm{mod}(\hslash {\omega }_F).$$

$$U^{†}(t)[q+i(K/{\omega }_F)p]U(t)=C(Q+iP)e^{-i{\omega }_{F}t/K}.$$

$$[Q,P]=i\lambda ,\lambda =\hslash K/{\omega }_F{C}^2.$$

$$g(Q,P)=\frac{1}{4}(Q^2+P^2-1{)}^2-\frac{1}{3}f(Q^3-3PQP),$$

$${\varphi }^{(k)}(Q)=\frac{1}{\sqrt{3(1+{\delta }^{(k)})}}\sum _{m=0,1,2}{\Psi }_m(Q)e^{-2mk\pi i/3}.$$

$$\displaystyle\int_{\infty}^{Q_{*}}dQ\Bigl{[}\phi^{(k)}(Q)(\hat{g}-g_{0})\Psi_{0}(Q)\Bigr{.}$$

$$\displaystyle\Bigl{.}-\Psi_{0}(Q)(\hat{g}-g^{(k)})\phi^{(k)}(Q)\Bigr{]}=0$$

$${\Psi }_m(Q)=C_m(i{\partial }_{P}g{)}^{-1/2}{e}^{i{S}_m(Q)/\lambda }(m=0,1),$$

$${\partial }_Q{S}_m=(-1{)}^m\bar{P}(Q),g(Q,\bar{P})=g_0,$$

$$g^{(k)}-g_0=C_{\mathrm{tun}}{e}^{-S_{\mathrm{tun}}/\lambda }\cos ({\lambda }^{-1}{\Phi }_{\mathrm{tun}}-2\pi k/3),$$

$$Q_0=\frac{1}{2}\left[f+(f^2+4{)}^{1/2}\right],Q_1=Q_2=-Q_0/2,$$

$$P_0=0,P_1=-P_2=\sqrt{3}{Q}_0/2,$$

$$g_{\min }=-\frac{1}{12}f{Q}_0(Q_0^2+3),{\omega }_{\min }=[3f{Q}_0(Q_0^2+1){]}^{1/2}.$$

$$g_0=g_{\min }+\frac{1}{2}\lambda {\omega }_{\min }.$$

$${\Psi }_0(Q)=(\sqrt{\pi }{l}_q{)}^{-1/2}\exp [-(Q-Q_0{)}^2/2l_q^2],$$

$${\Psi }_0(Q)=C_0(i{\partial }_{P}g{)}^{-1/2}\exp [i{S}_0(Q)/\lambda ],$$

$$S_0(Q)={\int }_{Q_0-l_q}^{Q}d{Q}^{\prime }\bar{P}(Q^{\prime }),$$

$$\bar{P}(Q{)}^2=A(Q)+B^{1/2}(Q),A(Q)=1-Q^2-2fQ,$$

$$B(Q)=A^2(Q)-4[g(Q,0)-g_0],$$

$${\partial }_{P}g=\bar{P}(Q)B^{1/2}(Q)$$

$$C_0=({\omega }_{\min }/2\sqrt{\pi e}{)}^{1/2}.$$

$$Q_B\approx Q_0-\frac{3}{4}f.$$

$$A(Q_1)=P_1^2>0,B(Q_1)=2\lambda {\omega }_{\min }.$$

$$\mathrm{Im}\bar{P}(Q)=-{\left[-A(Q)-B^{1/2}(Q)\right]}^{1/2}(Q_B<Q<Q_0)$$

$$\mathrm{Im}\bar{P}(Q)=-{\left[(A^2+|B|{)}^{1/2}-A\right]}^{1/2}/\sqrt{2}(Q<Q_B).$$

$${\Psi }_0(Q)\approx 2C_0|{\partial }_{P}g{|}^{-1/2}\exp [-\mathrm{Im}{S}_0(Q)/\lambda ]\cos {\Phi }_0(Q),$$

$${\Phi }_0(Q)={\Phi }_0^{\prime }(Q)+{\Phi }_0^{\prime \prime }(Q).$$

$${\Phi }_0^{\prime }(Q)={\lambda }^{-1}{\int }_{Q_B}^{Q}d{Q}^{\prime }\mathrm{Re}\bar{P}(Q^{\prime }),$$

$$\mathrm{Re}\bar{P}(Q)=-{\left[(A^2+|B|{)}^{1/2}+A\right]}^{1/2}/\sqrt{2},$$

Full text

The breaking of translation symmetry in time, first proposed by Wilczek Wilczek (2012), has been attracting much attention recently. Such symmetry breaking can occur only away from thermal equilibrium Watanabe and Oshikawa (2015). It is of particular interest for periodically driven systems, which have a discrete time-translation symmetry imposed by the driving. Here, the time symmetry breaking is manifested in the onset of oscillations with a period that is a multiple of the driving period $t_{F}$. Oscillations with period $2t_{F}$ due to simultaneously initialized protected boundary states were studied in photonic quantum walks Kitagawa et al. (2012); period-two oscillations can also be expected from the coexistence of Floquet Majorana fermions with quasienergies [math] and $\hbar\pi/t_{F}$ in a cold-atom system Jiang et al. (2011). The onset of period-two phases was predicted and analyzed Khemani et al. (2016); von Keyserlingk and Sondhi (2016); Else et al. (2016); Yao et al. (2017); Khemani et al. (2016); Bairey et al. (2017) in Floquet many-body localized systems, and the first observations of oscillations at multiples of the driving period in disordered systems were reported Zhang et al. (2016); Choi et al. (2016).

In systems coupled to a thermal bath, on the other hand, the effect of period doubling has been well-known. A textbook example is a classical oscillator modulated close to twice its eigenfrequency and displaying vibrations with period $2t_{F}$ Landau and Lifshitz (2004). The oscillator has two states of such vibrations; they have opposite phases, reminiscent of a ferromagnet with two orientations of the magnetization. Several aspects of the dynamics of a parametric oscillator in the quantum regime have been studied theoretically, cf. Wolinsky and Carmichael (1988); Drummond and Kinsler (1989); Wielinga and Milburn (1993); Kryuchkyan and Kheruntsyan (1996); Marthaler and Dykman (2006); Wustmann and Shumeiko (2013); Goto (2016); Puri and Blais (2016), and in experiments, cf. Nabors et al. (1990); Wilson et al. (2010); Lin et al. (2014). For a sufficiently strong driving field, a quantum dissipative oscillator, like a classical oscillator, mostly performs vibrations with period $2t_{F}$. The interplay of quantum fluctuations and dissipation leads to transitions between the period-two vibrational states, but the rate of these transitions is exponentially small Marthaler and Dykman (2006).

The goal of this paper is to study time symmetry breaking in isolated or almost isolated driven quantum systems with a few degrees of freedom. They are intermediate between large closed systems and dissipative systems, where the nature of the symmetry breaking is very different. To this end, we analyze a driven nonlinear quantum oscillator. Time symmetry breaking in this system should not be limited to period doubling. As an illustration of a behavior qualitatively different from period doubling, we consider period tripling and find the conditions where it occurs. We also address the role of decoherence and the connection between the time symmetry breaking in the coherent and incoherent regimes.

Floquet (quasienergy) states $\psi_{\varepsilon}(t)$ are eigenstates of the operator $T_{t_{F}}$ of time translation by $t_{F}$, $T_{t_{F}}\psi_{\varepsilon}(t)\equiv\psi_{\varepsilon}(t+t_{F})=\exp(-i\varepsilon t_{F}/\hbar)\psi_{\varepsilon}(t)$. For a broken-symmetry state $\psi_{K,\varepsilon_{K}}$ with $K>1$, time translation by $t_{F}$ is not described by the factor $\exp(-i\varepsilon t_{F}/\hbar)$. Instead, $\psi_{K,\varepsilon_{K}}(t+Kt_{F})=\exp(-Ki\varepsilon_{K}t_{F}/\hbar)\psi_{K,\varepsilon_{K}}(t)$. We call $\psi_{K,\varepsilon_{K}}$ a period-$K$ Floquet state. It is an eigenstate of $T_{Kt_{F}}=(T_{t_{F}})^{K}$, but not $T_{t_{F}}$.

Multiple-period states naturally occur if the number of states of the system ${\mathbb{N}}\to\infty$, as in the case of an oscillator. For such systems the quasienergy spectrum is generally dense, cf. Hone et al. (1997). Then we can find states $\psi_{\varepsilon}$ and $\psi_{\varepsilon^{\prime}}$ with the difference of the quasienergies $|\varepsilon-\varepsilon^{\prime}|$ infinitesimally close to $\hbar\omega_{F}/K$ with integer $K{}>1$ (or to $\hbar\omega_{F}k/K$ with $k<K$); $\omega_{F}=2\pi/t_{F}$ is the driving frequency. A linear combination $\alpha\psi_{\varepsilon}(t)+\alpha^{\prime}\psi_{\varepsilon^{\prime}}(t)$ is a period-$K{}$ state. The expectation value of dynamical variables in such a state oscillates with period $Kt_{F}$. However, the oscillation amplitude will be very small as, generally, the functions $\psi_{\varepsilon}$ and $\psi_{\varepsilon^{\prime}}$ will be of a very different form.

The situation is different for an oscillator driven close to an overtone of its eigenfrequency $\omega_{0}$, i.e., for $\omega_{F}\approx K\omega_{0}$. Such an oscillator has several sets of quasienergy states where the quasienergy differences within a set are very close to $\hbar\omega_{F}/K$ in a broad parameter range, and are exactly equal to $\hbar\omega_{F}/K$ for some interrelations between the parameters, whereas off-diagonal matrix elements of the dynamical variables are large, see Fig. 1. Such states result from tunnel splitting of the states localized at the minima of the oscillator Hamiltonian in the rotating frame shown in Fig.1(c). These localized states correspond to period-$K$ vibrations in the laboratory frame, see below.

In a way, for a parametric oscillator ($K=2$) the occurrence of a period-2 state could be inferred from the results Marthaler and Dykman (2007). However, this state was not identified there and the time symmetry breaking was not addressed. In different terms, sets of states separated by $\approx\hbar\omega_{F}/K$ were found numerically for $K\gg 1$ for a special model of an oscillator in the interesting paper Guo et al. (2013); the considered states did not break time symmetry.

The period tripling ($K=3$) considered here for a driven oscillator is particularly interesting. It differs from the continuous Landau-type symmetry-breaking transition that occurs for period doubling, cf. Lin et al. (2015). In the presence of dissipation, the fully-symmetric (zero-amplitude) state does not loose stability. Also, in the quantum regime, there emerges a geometric phase between the broken-symmetry states localized at the minima of the effective Hamiltonian function in phase space, cf. Fig. 1(c). Thus, the period-tripling in an oscillator allows one to reveal, using a simple and physically relevant model, the generic conditions for the onset of strongly overlapping multiple-period states and to relate them to the underlying nontrivial symmetry. It also provides a platform for studying quantum tunneling between localized states in phase space. This problem is considerably different from the classical problem of tunneling in a symmetric double-well potential Landau and Lifshitz (1997) (see also Garg (2000)).

We study a most commonly used model of a nonlinear oscillator, the Duffing model, which describes parametric resonance and, as we will see, can describe period tripling; this model refers to a broad range of systems, including trapped relativistic electrons, cold atomic clouds, Josephson junction based systems, and nanomechanical systems Tan and Gabrielse (1991); Dykman (2012). Its Hamiltonian reads

[TABLE]

where $q$ and $p$ are the oscillator coordinate and momentum. The term $H_{F}\equiv H_{F}(t)$ describes the driving. In the analysis of parametric resonance, one chooses $H_{F}=-\tfrac{1}{2}q^{2}F\cos\omega_{F}t$ with $\omega_{F}\approx 2\omega_{0}$. Here we consider $H_{F}=-\tfrac{1}{3}q^{3}F\cos\omega_{F}t$ with $\omega_{F}\approx 3\omega_{0}$; the results describe also a drive $H_{F}^{\prime}=-qF^{\prime}\cos\omega_{F}t$ with $F\to 3\gamma F^{\prime}/8\omega_{0}^{2}$.

If the driving is not too strong, so that for the states of interest the expectation values of $H_{F}$ and the nonlinear term $\propto q^{4}$ are small compared to the harmonic part of $H_{0}$, the resonant oscillator dynamics can be described in the rotating wave approximation (RWA) Walls (2008). For an oscillator driven close to the $K{}$th overtone of its eigenfrequency, one makes a canonical transformation $U(t)=\exp(-ia^{\dagger}a\omega_{F}t/K)$, where $a$ and $a^{\dagger}$ are the ladder operators. The RWA Hamiltonian $H_{\rm RWA}$ is obtained by time-averaging the transformed Hamiltonian $H_{K}(t)=U^{\dagger}(t)H(t)U(t)-i\hbar U^{\dagger}(t)\dot{U}(t)$,

[TABLE]

Clearly, $H_{\rm RWA}$ is independent of time.

We now establish the relation between the eigenvalues of $H_{\rm RWA}$ and the quasienergies. If $\phi(t)$ is an eigenfunction of $H_{\rm RWA}$, i.e., $H_{\rm RWA}\phi=E\phi$, then the corresponding wave function in the lab frame is $\psi(t)=U(t)\phi(t)$, and

[TABLE]

We call $E$ the RWA energy. In Eq. (3)

[TABLE]

The above commutation relation follows from the relation $H_{K}(t+t_{F})=N_{K}^{\dagger}H_{K}(t)N_{K}$ and Eq. (2). Using the explicit form of $H_{\rm RWA}$, the commutation relation (4) was found in Ref. Guo et al., 2013 for the same operator as $N_{K}$.

Operators $N_{K}^{k}$ with $k=0,1,...,K-1$ form a cyclic group. Since eigenfunctions of $H_{\rm RWA}$ are also eigenfunctions of $N_{K}$, one can label them by a superscript $k$,

[TABLE]

Note that $H_{\rm RWA}$ has eigenfunctions with the same $k$, but different $E$. By comparing Eqs. (3) and (5) one finds that a wave function $\phi^{(k)}$ with RWA energy $E^{(k)}$ corresponds to a usual Floquet state with quasienergy

[TABLE]

As we will see, for sufficiently strong drive the eigenstates of $H_{\rm RWA}$ form multiplets with close eigenvalues $E^{(k)}$ but different $k$. The quasienergies of different states in the multiplets differ by $\approx\hbar\omega_{F}/K$.

Equation (5) allows one to write the functions $\phi^{(k)}$ in terms of the Fock states of the oscillator $|n\rangle$ defined by the condition $a^{\dagger}a|n\rangle=n|n\rangle$. Only one out of each $K$ Fock states contributes to $\phi^{(k)}$, $\phi^{(k)}=\sum_{n}C_{n}^{(k)}|Kn+k\rangle$. This relation significantly simplifies numerical diagonalization of $H_{\rm RWA}$, as the coefficients $C_{n}^{(k)}$ with different $k$ are uncoupled. More importantly, it shows that the RWA energy levels of states with different $k$ can cross when the parameters of the system vary. This crossing is seen in Fig. 1. In contrast, the RWA levels of states with the same $k$ avoid crossing.

The motion in the rotating frame is conveniently described by the coordinate $Q$ and momentum $P$, which are related to $q$ and $p$ as

[TABLE]

The parameter $C$ is the scaling factor that makes $Q$ and $P$ dimensionless,

[TABLE]

The dimensionless Planck constant $\lambda$ and the parameter $C$ in the case of a parametric oscillator, $K=2$, are given in Marthaler and Dykman (2006). For the case of period tripling, $C=(8\omega_{F}\delta\omega/9\gamma)^{1/2}$, where $\delta\omega=\frac{1}{3}\omega_{F}-\omega_{0}$ is the frequency detuning from the resonance, $|\delta\omega|\ll\omega_{F}$. In this case $H_{\rm RWA}=[8\omega_{F}^{2}(\delta\omega)^{2}/27\gamma)]\hat{g}(Q,-i\lambda\partial_{Q})$ with

[TABLE]

where $f=F/(8\omega_{F}\gamma\delta\omega)^{1/2}$ is the scaled amplitude of the driving. Of interest is the region $\gamma\delta\omega>0$, and we choose $\gamma>0$ and $\delta\omega>0$.

The function $g(Q,P)$ is the dimensionless Hamiltonian function in the rotating frame. It is plotted in Fig. 1. It has a three-fold rotational symmetry in the $(Q,P)$-plane. This symmetry follows from Eqs. (4) and (7), since $N_{K}$ is an operator of rotation by angle $2\pi/K$ in phase plane; the $K$-fold symmetry of $H_{\rm RWA}$ was also seen in Guo et al. (2013).

For moderately strong fields, $g(Q,P)$ has three well-separated minima positioned at the vertices of an equilateral triangle $(Q_{m},P_{m})$; we count $m=0,1,2$ counterclockwise and set $m=0$ for the vertex with $P_{0}=0$. The eigenstates of the operator $\hat{g}\equiv g(Q,-i\lambda\partial_{Q})$ with the lowest RWA energies are localized near $(Q_{m},P_{m})$. In the absence of tunneling, $\hat{g}$ has three degenerate eigenstates $\Psi_{m}$. Near their maxima, functions $\Psi_{m}$ have the form of squeezed ground states of a harmonic oscillator centered at $(Q_{m},P_{m})$ 111see Supplemental Material for the details of the calculation.

The oscillator in a state $\Psi_{m}$ has a broken time symmetry. The expectation values of dynamical variables oscillate at frequency $\omega_{F}/3$. Indeed, from Eq. (7) time translation by $t_{F}$ transforms $\Psi_{m}\to N_{3}\Psi_{m}=\Psi_{m-1}\equiv\Psi_{m+2}$. To come back to state $\Psi_{m}$, one has to increment time by $3t_{F}$. The relation $\Psi_{m+1}=N_{3}^{\dagger}\Psi_{m}$ gives the phase shift between functions $\Psi_{m+1}$ and $\Psi_{m}$. Since $N_{3}$ is a rotation operator, this phase shift is geometric in nature [34].

Tunneling between the minima lifts the degeneracy of the ground state of the operator $\hat{g}$. In contrast to the problem of tunneling in a symmetric double-well potential Landau and Lifshitz (1997), $g(Q,P)$ is not even in $Q$, it has three extrema, and two of them lie at nonzero momenta $P$.

To find the tunnel splitting, we write the wave functions in the coordinate representation, $\Psi_{m}\equiv\Psi_{m}(Q)$. The three normalized eigenstates $\phi^{(k)}$ of $\hat{g}$ with the smallest eigenvalues $g^{(k)}$ ($k=0,1,2$) have the form

[TABLE]

where $\delta^{(k)}=2{\rm Re}[\langle\Psi_{0}|\Psi_{1}\rangle\exp(-2\pi ik/3)]\ll 1$. We choose $\Psi_{0}(Q)$ to be real and normalized. Since $\Psi_{m+1}=N_{3}^{\dagger}\Psi_{m}$, we have $\Psi_{2}(Q)=\Psi_{1}^{*}(Q)$. Due to the symmetry, the functions $\phi^{(k)}$ can be shown to be orthogonal.

In the spirit of Landau and Lifshitz (1997), we calculate $g^{(k)}$ using the relation

[TABLE]

with $g_{0}$ being the eigenvalue of $\hat{g}$ in the state $\Psi_{0}$, $g_{0}\approx\min g(Q,P)$ [34]. The difference $g^{(k)}-g_{0}$ is exponentially small for a small dimensionless Planck constant $\lambda$.

To choose the upper limit $Q_{*}$ of the integral (Multiple-period Floquet states and time-translation symmetry breaking in quantum oscillators), we note that the functions $\Psi_{m}(Q)$ fall off exponentially away from the respective $Q_{m}$, with $\Psi_{0}$ and $\Psi_{1,2}$ falling off in the opposite directions in the interval $(Q_{1},Q_{0})$. We choose $Q_{*}$ within this interval and in such a way that $\Psi_{0,1,2}(Q_{*})$ are all of the same order of magnitude and thus can be kept in Eq. (10) for $\phi^{(k)}(Q)$. The result of integration (Multiple-period Floquet states and time-translation symmetry breaking in quantum oscillators) should be independent of $Q_{*}$.

The WKB wave functions $\Psi_{0,1}(Q)$ in the classically forbidden region between $Q_{1}$ and $Q_{0}$ have the form

[TABLE]

where $S_{0,1}(Q)$ is the classical action and constants $C_{0,1}$ are found from the matching to the corresponding intrawell wave functions.

It is critical for understanding the tunneling that, because the effective Hamiltonian function $g(Q,P)$ is quartic in the momentum $P$, $\bar{P}(Q)$ has a branch point $Q_{B}$ in the interval $(Q_{1},Q_{0})$. For $Q_{1}<Q<Q_{B}$, $\bar{P}(Q)$ has both imaginary and real parts. So does the action $S_{m}(Q)$. This leads to oscillations of the wave functions in the classically forbidden region. In $S_{m}(Q)$ one should keep the root with the smallest $|{\rm Im}~{}\bar{P}|$. To describe $\Psi_{0}$, Eq. (Multiple-period Floquet states and time-translation symmetry breaking in quantum oscillators) has to be modified by allowing for a complex conjugate term [34].

Calculating the integrals in Eq. (Multiple-period Floquet states and time-translation symmetry breaking in quantum oscillators) by parts, we find

[TABLE]

where $\Phi_{\rm tun}+iS_{\rm tun}=\int_{Q_{0}}^{Q_{1}}dQ^{\prime}P_{\rm cl}(Q^{\prime})+P_{1}Q_{1}/2+\lambda G$ with $P_{\rm cl}$ given by equation $g(Q,P_{\rm cl})=\min g(Q,P)$, $G$ being independent of $\lambda$ and having a contribution from the geometric phase, and $C_{\rm tun}\propto\lambda^{1/2}$ [34].

Equation (13) shows that the splitting of the eigenvalues of $H_{RWA}$ oscillates as the system parameters vary. Two eigenvalues cross each time $\lambda^{-1}\Phi_{\rm tun}=(n+n^{\prime}/3)\pi$ with integer $n,n^{\prime}$. Such crossings are seen in Fig. 1. Where the eigenvalues do not cross, they stay exponentially close to each other.

If the oscillator is in a superposition of two states $\phi^{(k)}$ and $\phi^{(k^{\prime})}$, the expectation values of its variables have period $3t_{F}$ provided the observation time is smaller than the exponentially long time $|\Omega_{kk^{\prime}}|^{-1}$, where the frequency $\Omega_{kk^{\prime}}=\lambda^{-1}[g^{(k)}-g^{(k^{\prime})}]\delta\omega$ is determined by the tunnel splitting. The Fourier spectra of the expectation values generally have components at frequency $\omega_{F}/3\pm\Omega_{kk^{\prime}}$; in particular, the coordinate and momentum have just one of these components. This behavior is characteristic also of the oscillator in intrawell states $\Psi_{m}$, which are superpositions of $\phi^{(1,2,3)}$. The oscillator fluorescence spectrum will display peaks at $\omega_{F}/3\pm\Omega_{kk^{\prime}}$ as well.

It is instructive to compare these results with the period-doubling associated with the topologically protected Floquet boundary states in extended systems Kitagawa et al. (2012); Jiang et al. (2011). To some extent, such states are analogous to the symmetry-protected states $\phi^{(k)}$. If tunneling between the Floquet boundary states can be disregarded, similar to disregarding oscillator tunneling, their combination becomes a multiple-period state. However, their overlap is exponentially small, in contrast to the functions $\phi^{(k)}$.

The intrawell states $\Psi_{m}$ are particularly important in the presence of dissipation. Even if the dissipation rate $\Gamma$ is extremely small, but exceeds the exponentially small frequencies $\Omega_{kk^{\prime}}$, instead of coherent tunneling between the wells of $g(Q,P)$, the oscillator performs incoherent interwell hopping with typical rate $W<|\Omega_{kk^{\prime}}|$ [34]. This hopping corresponds to flips of the vibration phase. On times small compared to $W^{-1}$ the oscillator stays in the multiple-period state inside a well. This is the exact analog of the classical behavior of a dissipative oscillator, including a parametric oscillator, where the multiple-period state is seen on times short compared to the reciprocal rate of interstate switching.

A promising type of oscillator for observing period tripling are modes of microwave cavities coupled to Josephson junctions. Recently there have been studied systems where inelastic Cooper pair tunneling leads to an effective driving of a cavity mode that nonlinearly depends on the mode coordinate and has a tunable frequency $2eV/\hbar$ determined by the voltage $V$ across the Josephson junction Hofheinz et al. (2011); Armour et al. (2013); Gramich et al. (2013). There are also other possibilities to resonantly excite multiple-period modes in microwave cavities 222P. Delsing, D. Esteve, and F. Portier, private communications.

In conclusion, we studied a quantum oscillator driven close to an overtone of its eigenfrequency and showed that a small quantum system can display coherent multiple-period dynamics. We explicitly described this dynamics for the previously unexplored nontrivial case of period tripling and established the relation to protected boundary Floquet states in extended systems and to multiple-period states in dissipative systems.

We are grateful to G. Refael, M. Rudner, and S. Sondhi for the discussions and correspondence. YZ and SMG were supported by the U.S. Army Research Office (W911NF1410011) and by the National Science Foundation (DMR-1609326).; JG and JA were supported in part by the German Science Foundation through SFB/TRR 21 and the Center for Integrated Quantum Science and Technology (IQST); MID was supported in part by the National Science Foundation (Grant No. DMR-1514591).

I The Intrawell Wave Functions of the RWA Hamiltonian

We consider the dynamics of the oscillator driven close to three times its eigenfrequency in the rotating wave approximation (RWA). The scaled RWA Hamiltonian function $g(Q,P)$, which is given by Eq. (10) of the main text and is plotted there in Fig. 1, has three symmetrically located minima at points $(Q_{m},P_{m})$ with $m=0,1,2$,

[TABLE]

The minimal value $g_{\min}$ of $g(Q,P)$ and the dimensionless frequency of classical vibrations about a minimum $\omega_{\min}=\left(\det[\partial^{2}_{x_{i}x_{j}}g(x_{1},x_{2})]\right)^{1/2}$ (the derivatives are calculated at a minimum of $g$) are

[TABLE]

The frequency $\omega_{\min}$ is the same for all minima. So is also the lowest eigenvalue $g_{0}$ of the Hamiltonian $\hat{g}(Q,-i\lambda\partial_{Q})$ in the neglect of tunneling. To the lowest order in the dimensionless Planck constant $\lambda$ it corresponds to the lowest eigenvalue of a harmonic oscillator,

[TABLE]

The calculation of the tunnel splitting is done below by first finding the intrawell wave functions $\Psi_{m}(Q)$ near their maxima inside the well, then finding the geometric phase shift between different $\Psi_{m}$, and then explicitly writing down the WKB tails of functions $\Psi_{m}$ in the classically forbidden regions, which are given by Eq. (13) of the main text. Since $\Psi_{2}(Q)=\Psi_{1}^{*}(Q)$, we only need to find $\Psi_{0}(Q)$ and $\Psi_{1}(Q)$.

I.1 The wave function $\Psi_{0}(Q)$

Near the minimum $(Q_{0},P_{0})$ we have $g(Q,P)\approx g_{\min}+\frac{1}{2}(Q_{0}^{2}+1)(Q-Q_{0})^{2}+\frac{3}{2}fQ_{0}P^{2}$. The wave function $\Psi_{0}(Q)$ is Gaussian for $|Q-Q_{0}|\ll|Q_{1}-Q_{0}|$ and can be chosen to be real,

[TABLE]

with $l_{q}=[\lambda\omega_{\min}/(Q_{0}^{2}+1)]^{1/2}$ being the localization length.

We are interested in the tail of $\Psi_{0}$ for $Q$ between the minima of $g(Q,P)$, i.e., for $Q_{1}<Q<Q_{0}-l_{q}$. The WKB form of $\Psi_{0}(Q)$ is given by Eq. (13) of the main text, which we here write explicitly,

[TABLE]

where $\bar{P}(Q)$ is given by equation $g(Q,\bar{P})=g_{0}$ and $\partial_{P}g$ is calculated for $P=\bar{P}(Q)$. For the branch of $\bar{P}$ that we are interested in

[TABLE]

with Im $\bar{P}<0$ for $Q<Q_{0}$; we keep the correction $\propto\lambda$ to secure matching to Eq. (16).

For $Q$ close to $Q_{0}$ and $Q<Q_{0}-l_{q}$, we have $A(Q)<0,B(Q)>0$, and $A(Q)+B^{1/2}(Q)<0$. Therefore $\bar{P}(Q)$ is purely imaginary and the same is true for the function

[TABLE]

with $i\partial_{P}g>0$. Accordingly, $\Psi_{0}(Q)$ exponentially decays with increasing $Q_{0}-Q$. The prefactor $C_{0}$ is determined by matching Eqs. (16) and (I.1) for $Q$ close to $Q_{0}$ but $Q_{0}-Q\gg l_{q}$,

[TABLE]

As $Q$ decreases, first $B(Q)$ becomes equal to zero at point $Q_{B}$. To the leading order in $\lambda\ll 1$

[TABLE]

For still smaller $Q$, $A(Q)$ changes sign to positive. This happens for $Q_{B}>Q>Q_{1}\equiv-Q_{0}/2$. Importantly,

[TABLE]

In the explicit form, the imaginary part of the momentum in the classically forbidden region is

[TABLE]

As discussed in the main text, the level splitting crucially depends on the oscillations of the wave function under the barrier. These oscillations start with the decreasing $Q$ at $Q=Q_{B}$. Near $Q_{B}$ we have $B(Q)\approx\partial_{Q}B(Q_{B})(Q-Q_{B})$, whereas $A(Q_{B})<0$. Therefore $\bar{P}\approx-i|A(Q_{B})|^{1/2}+(i/2)|\partial_{Q}B(Q_{B})/A(Q_{B})|^{1/2}(Q-Q_{B})^{1/2}$ for small $Q-Q_{B}>0$, i.e., $Q_{B}$ is a branching point of $\bar{P}(Q)$. We have to go around above and below this point in the complex plane to obtain the wave function for $Q<Q_{B}$, following the standard procedure Landau and Lifshitz (1997). As a result, we find for $Q<Q_{B}$

[TABLE]

Here, the phase $\Phi^{\prime}_{0}(Q)$ comes from the real part of the action,

[TABLE]

whereas $\Phi_{0}^{\prime\prime}(Q)$ comes from the prefactor, with account taken of going around $Q_{B}$ in the complex plane,

[TABLE]

The choice of ${\rm Re}\bar{P}$ and ${\rm Im}\bar{P}$ in Eqs. (I.1) and (I.1) corresponds to writing $B^{1/2}=i|B|^{1/2}$ in Eq. (I.1) for $\bar{P}^{2}$ in the region where $B(Q)<0$.

The WKB approximation (I.1) breaks down near $Q_{1}$, as $B(Q)$ becomes $\sim\lambda$ and $|\partial_{P}g|$ becomes small. However, we do not need to calculate the wave function $\Psi_{0}(Q)$ in this region, as seen from Eq. (12) of the main text.

I.2 The wave function $\Psi_{1}(Q)$

The minimum of $g(Q,P)$ at $(Q_{1},P_{1})$ corresponds to a nonzero momentum $P_{1}>0$. Therefore the wave function $\Psi_{1}$ centered at $Q_{1}$ is complex valued even near its maximum. Calculating $\Psi_{1}$ involves three steps: finding it inside the well of $g(Q,P)$ near $Q_{1},P_{1}$; finding the geometric phase, that relates $\Psi_{1}$ and $\Psi_{0}$ given that $\Psi_{0}$ is chosen in the form (16), and then finding the tail of $\Psi_{1}$ in the classically forbidden range.

I.2.1 The intra-well wave function and the geometric phase

Using the explicit form (I) of $Q_{1},P_{1}$, to the second order in $\delta Q=Q-Q_{1},\delta P=P-P_{1}$ we write the Hamiltonian near $(Q_{1},P_{1})$ as

[TABLE]

The expression for $\Psi_{1}$ for $|\delta Q|\ll Q_{0}-Q_{1}$ then reads

[TABLE]

The Gaussian-width parameter $\beta$ is now complex-valued. So is also the prefactor $C_{1,{\rm intra}}$, which has a phase factor $\exp(i\theta_{1})$.

The phase $\theta_{1}$ has a geometric nature. It is determined by the fact that, as indicated in the main text, $\Psi_{1}$ and $\Psi_{0}$ are related by the transformation of rotation in phase plane, $\Psi_{1}=N_{3}^{\dagger}\,\Psi_{0}$. Here, $N_{3}=\exp(-2\pi ia^{\dagger}a/3)$ with $a=(2\lambda)^{-1/2}(Q+iP)\equiv(2\lambda)^{-1/2}(Q+\lambda\partial_{Q})$. To calculate $\theta_{1}$, we consider a coherent state in the coordinate representation

[TABLE]

and set $\alpha=Q_{0}/\sqrt{2\lambda}$, so that the wave function $|\alpha\rangle$ is centered at $Q_{0}$ and thus strongly overlaps with $\Psi_{0}$. Since the function $\Psi_{1}$ is obtained from $\Psi_{0}$ by applying to $\Psi_{0}$ the operator $N_{3}^{\dagger}$, we can write the overlap integral as $\langle\alpha|\Psi_{0}\rangle=\langle\alpha|N_{3}\Psi_{1}\rangle=\langle\alpha\exp(2\pi i/3)|\Psi_{1}\rangle$. The “rotated” state $|\alpha\exp(2\pi i/3)\rangle$ strongly overlaps with $\Psi_{1}$. Therefore the above overlap integrals can be calculated using the explicit Gaussian form of $\Psi_{0}$ and $\Psi_{1}$ near their maxima. With account taken of the normalization of $\Psi_{1}$, this gives

[TABLE]

I.2.2 The wave function $\Psi_{1}$ in the classically forbidden region

It is clear from Eq. (12) of the main text that we need to find the tail of the wave function $\Psi_{1}$ in the classically forbidden region only for $Q>Q_{1}$. It is given by Eq. (13) of the main text. In a more explicit form

[TABLE]

where $\bar{P}(Q)$ is given by Eqs. (I.1) and (I.1), $l_{q}^{\prime}=[\lambda/{\rm Re}\,\beta]^{1/2}$ . Equation (29) corresponds to choosing $B^{1/2}(Q)=i|B(Q)|^{1/2}$ for $B(Q)<0$ and to $\partial_{P}g$ calculated for $P(Q)=\bar{P}(Q)$, i.e., $\partial_{P}g=\bar{P}B^{1/2}$. For $Q_{B}-Q\gg Q-Q_{1}\gg l^{\prime}_{q}$ we have $-\bar{P}(Q)\approx P_{1}+i\beta(Q-Q_{1})$, as expected from Eq. (I.2.1). By matching Eqs. (I.2.1) and (29), we find

[TABLE]

Because we count the action $S_{1}$ off from $Q_{1}+l_{q}^{\prime}$, there emerges an extra phase factor in $C_{1}$ due to the oscillations of the wave function inside the “potential well” centered at $(Q_{1},P_{1})$.

II Tunnel splitting of the scaled RWA energy levels

The scaled RWA energies $g^{(k)}$ give the values of the quasienergies $\varepsilon^{(k)}$ of the driven oscillator, $\varepsilon^{(k)}=(\Xi g^{(k)}+\hbar\omega_{F}k/3)\!\!\!\mod(\hbar\omega_{F})$, where $\Xi=|8\omega_{F}^{2}(\delta\omega)^{2}/27\gamma|$, see Eqs. (6) and the text above Eq. (9) of the main text. The explicit expressions for the wave functions (I.1) and (29) allow us to calculate the scaled energies $g^{(k)}$ using Eq. (12) of the main text, which we reproduce here for convenience,

[TABLE]

Functions $\phi^{(k)}$ are sums of functions $\Psi_{m}(Q)$ weighted with factors $\exp(-2\pi imk/3)/\sqrt{3}$. For $Q_{*}$ well inside the interval $(Q_{1},Q_{0})$ , we have $\int_{\infty}^{Q_{*}}\Psi_{0}^{2}(Q)dQ=-1$. Taking into account that overlapping of the functions $\Psi_{1,2}(Q)$ with $\Psi_{0}(Q)$ is exponentially small, we rewrite Eq. (II) as

[TABLE]

It is important that the product $\Psi_{0}(Q)\Psi_{1}(Q)$ has two terms. One of them is $\propto\exp\{i[S_{0}(Q)+S_{1}(Q)]\}$. It smoothly depends on $Q$, because $S_{0}(Q)+S_{1}(Q)={\rm const}$ for $Q_{1}<Q<Q_{0}$. The other term is $\propto\exp\{-i[S_{0}^{*}(Q)-S_{1}(Q)]\}$, it is a fast oscillating function of $Q$. The contribution of this term to the integrals (II) is exponentially small and exponentially sensitive to the change of $Q_{*}$ on the scale $\propto\lambda$. Therefore this term should be disregarded.

Using the explicit form of the operator $g(Q,-i\lambda\partial_{Q})$ and integrating by parts, from Eq. (II) we obtain

[TABLE]

This expression is somewhat inconvenient, as $\bar{P}$ is calculated with account taken of the term $\propto\lambda$. It is easy to see that $\bar{P}(Q)\approx{}P_{\rm cl}(Q)+\frac{1}{2}\lambda\omega_{\min}/\partial_{P}g$, where $P_{\rm cl}$ is given by the value of $\bar{P}$ calculated for $\lambda=0$. This approximation breaks down near $Q_{0},Q_{B}$ and $Q_{1}$ where $\partial_{P}g$ goes to zero. Similar to Ref. Garg, 2000, for $Q_{0}>Q>Q_{B}$ one can write

[TABLE]

Here, ${}B_{\rm cl}(Q)=(16f/3)(Q-Q_{1})^{2}(Q-Q_{B})$ is the value of $B(Q)$ calculated for $\lambda=0$. A similar transformation can be made for $\int_{Q_{1}+l_{q}^{\prime}}^{Q}dQ^{\prime}\bar{P}(Q^{\prime})$ in the region $Q_{1}<Q<Q_{B}$.

We now have to consider the vicinity of $Q_{B}$. Formally, the quantum correction to $\bar{P}$ diverges at $Q_{B}$. However, the divergence is integrable. Therefore Eq. (II) applies all the way till $Q=Q_{B}$, and one can use the value of $Q_{B}$ given by Eq. (20).

The final result for the difference of the scaled RWA energies is Eq. (14) of the main text,

[TABLE]

with real $S_{\rm tun}$ and $\Phi_{\rm tun}$,

[TABLE]

Here,

[TABLE]

and $\theta_{1}$ is given in (I.2.1).

The explicit expression (35) is in an extremely good agreement with the numerical calculations. This can be seen from Fig. 1 in the main text. A more detailed comparison is shown in Fig. 2. Equation (35) simplifies in the limit of comparatively strong drive, $f\gg 1$. The leading order terms in $S_{\rm tun}$ and in $\Phi_{\rm tun}$ are quadratic in $f$. Numerically, the asymptotic regime is reached for comparatively large $f$, where the tunneling amplitude becomes very small.

III Quantum Diffusion over the broken-symmetry states

The dynamics of the driven oscillator system can be strongly changed by an already very weak dissipation. Two types of dissipative processes can be conditionally separated. One of them causes transitions between the states that belong to the same multiplet formed by the tunnel splitting of a quantized state of motion inside a well of $g(Q,P)$. In particular, in this paper we considered such multiplet $\phi^{(k)}$ formed by the tunnel splitting of the lowest quantized intrawell state. The other dissipative process leads to transitions between the intrawell states.

In terms of the dissipation mechanisms, an important type of physical dissipative processes are transitions between the Fock states of the oscillator with emission or absorption of excitations of the thermal reservoir to which the oscillator is coupled. Another mechanism is fluctuations of the oscillator eigenfrequency due to the coupling to a reservoir or due to an external noise. It leads to dephasing of the vibrations, but not to an appreciable energy exchange with the reservoir. There may be also dissipation channels that are induced by the driving field; however, for the considered comparatively weak resonant field they are not important.

We note first that the dephasing does not mix the states within the tunnel-split multiplets. Indeed, as indicated in the main text, the wave functions $\phi^{(k)}$ can be written in terms of the Fock states of the oscillator $|n\rangle$ as $\phi^{(k)}=\sum_{n}C_{n}^{(k)}|3n+k\rangle$. The coupling to a thermal bath, which leads to dephasing, has the form $H^{\rm(ph)}=a^{\dagger}aH_{b}^{\rm(ph)}$, where $a,a^{\dagger}$ are the oscillator ladder operators and $H_{b}^{\rm(ph)}$ is an operator that depends on the dynamical variables of the bath only. Clearly, such coupling is diagonal in the $\phi^{(k)}$ basis.

The simplest coupling that leads to the oscillator energy relaxation is linear in $a,a^{\dagger}$. To the lowest order of the perturbation theory, for the well-understood conditions, it is described by the term $\dot{\rho}_{d}$ in the equation for the oscillator density matrix $\rho$ in slow time compared to $\omega_{F}^{-1}$,

[TABLE]

where $2\Gamma$ is the energy decay rate of the oscillator and $\bar{n}=[\exp(\hbar\omega_{0}/k_{B}T)-1]^{-1}$ is the oscillator Planck number. In what follows we assume that $\bar{n}=0$; an extension to a nonzero Planck number is straightforward and does not affect the result.

The goal of this section is to show that, even where $\Gamma$ is extremely small, but exceeds the tunnelling frequency $\Omega_{kk^{\prime}}=\lambda^{-1}\delta\omega(g^{(k)}-g^{(k^{\prime})})$, the oscillator dynamics changes qualitatively compared to the coherent dynamics. Instead of coherent tunneling between the intrawell states $\Psi_{m}$, which have broken time translation symmetry, the oscillator performs random hopping between the wells.

We first discuss the effect of the dissipation (III) by disregarding the dissipation-induced transitions between the intrawell states. In this approximation, one can describe the evolution of the oscillator in terms of the kinetic equation for the matrix elements $\rho_{mm^{\prime}}\equiv\langle\Psi_{m}|\rho|\Psi_{m^{\prime}}\rangle$. The interwell tunneling can be mapped onto the tight-binding model with Hamiltonian

[TABLE]

where we use the convention $|\Psi_{3}\rangle\equiv|\Psi_{0}\rangle$. The hopping integral is $t_{\rm tun}=(\hbar\delta\omega/2\lambda)C_{\rm tun}\exp[(-S_{\rm tun}+i\Phi_{\rm tun})/\lambda]$ with $C_{\rm tun},S_{\rm tun}$, and $\Phi_{\rm tun}$ given by Eq. (II).

To the leading order in $\lambda$, we have $\langle\Psi_{m}|a|\Psi_{m^{\prime}}\rangle=(2\lambda)^{-1/2}(Q_{m}+iP_{m})\delta_{mm^{\prime}}$. Therefore, from Eq. (III), off-diagonal matrix elements $\rho_{mm^{\prime}}$ decay with rate $\propto\Gamma/\lambda$. If this rate exceeds $|\Omega_{kk^{\prime}}|\sim|t_{\rm tun}|/\hbar$, then over time $\sim\lambda/\Gamma$ the off-diagonal matrix elements decay to their quasi-stationary values, which are determined by the diagonal matrix elements $\rho_{mm}$. The latter vary much slower,

[TABLE]

Parameter $W$ is the rate of hopping between the wells of $g(Q,P)$, it is much smaller than the tunneling frequency $|\Omega_{kk^{\prime}}|$. The hopping is a Poisson process in the slow time, it is incoherent and is a discrete analog of diffusion. The above analysis is in the spirit of the theory of quantum diffusion in solids Kagan (1992) and its analog in systems with a small number of potential wells Dykman and Tarasov (1978).

The role of the dissipation-induced intrawell transitions is more subtle. Even for $T=0$, these transitions lead to an occupation of excited intrawell states, cf. Marthaler and Dykman (2006). On the time scale determined by $1/\Gamma$, near the minimum of a well there is progressively formed a Boltzmann-type distribution over the states. The stationary ratio of the populations of the neighboring states can be shown to be $(1+2fQ_{0}-\omega_{\min})/(1+2fQ_{0}+\omega_{\min})$.

The tunnel splitting increases for higher-lying intrawell states. However, near the minimum of $g(Q,P)$ this increase is slow. The tunneling action $S_{\rm tun}(n)$ varies with the intrawell level number $n$ as $|\partial S_{\rm tun}(n)/\partial n|=\lambda\omega_{\min}\tau_{n}$. Here, $\tau_{n}$ is the dimensionless imaginary time of interwell tunneling given by Im $\int dQ/\partial_{P}g$, where the classical momentum is calculated for $g(Q,P)=g_{\min}+\lambda\omega_{\min}(n+1/2)$. This time is logarithmically large for small $n$. Therefore, for small $\Gamma$ but still $\Gamma\gg|t_{\rm tun}|/\hbar$, tunneling via excited intrawell states weakly renormalizes the rate $W$ in Eq. (III).

Even if for highly excited intrawell states, with $n$ exceeding some critical $n_{\rm cr}\propto 1/\lambda$, the hopping integral exceeds $\Gamma/\hbar$, interwell switching via these states will be very slow, as the occupation of these states will be small. We note that, if $\hbar\Gamma$ exceeds the hopping integral for almost all intrawell states, interwell switching may occur via dissipation-induced transitions over the interwell barrier of $g(Q,P)$, i.e., over the saddle point of $g(Q,P)$ seen in Fig. 1 of the main text. This is the dominating switching mechanism for a parametric oscillator Marthaler and Dykman (2006).