Preparing quasienergy states on demand: A parametric oscillator

TL;DR

This paper demonstrates a method to prepare specific quasienergy states in a nonlinear parametric oscillator using adiabatic techniques and Landau-Zener transitions, analyzing the resulting quantum states and emitted radiation.

Contribution

It introduces a controlled protocol for on-demand quasienergy state preparation in a nonlinear oscillator, including state characterization and emission spectrum analysis.

Findings

Successful preparation of even quasienergy states from the ground state.

Ability to create superpositions via Landau-Zener transitions.

Characterization of the emitted radiation spectrum.

Abstract

We study a nonlinear oscillator, which is parametrically driven at a frequency close to twice its eigenfrequency. By judiciously choosing the frequency detuning and linearly increasing the driving amplitude, one can prepare any even quasienergy state starting from the oscillator ground state. Such state preparation is effectively adiabatic. We find the Wigner distribution of the prepared states. For a different choice of the frequency detuning, the adiabaticity breaks down, which allows one to prepare on demand a superposition of quasienergy states using Landau-Zener-type transitions. We find the characteristic spectrum of the transient radiation emitted by the oscillator after it has been prepared in a given quasienergy state.

Figures

Equations

$$H(t)=\frac{p^2}{2}+\frac{1}{2}{q}^2[{\omega }_0^2+F\cos ({\omega }_{F}t)]+\frac{\gamma }{4}{q}^4.$$

$${\psi }_{\epsilon }(t)=e^{-i\epsilon t/\hslash }{u}_{\epsilon }(t),u_{\epsilon }(t+t_F)=u_{\epsilon }(t),$$

$$H_{\mathrm{RWA}}=-\hslash \delta {\omega }_F\hat{n}+\frac{\hslash V}{2}({\hat{n}}^2+\hat{n})+\frac{\hslash \tilde{F}}{2}(a^2+a^{+2})$$

$$\epsilon =[E+(1-P_E)\hslash {\omega }_F/4]\mathrm{mod}(\hslash {\omega }_F),$$

$$u_{\epsilon }(t)=\exp [i(1-P_E){\omega }_{F}t/4]U(t){\varphi }_E,$$

$$\displaystyle\text{where }P_{E}=\left\{\begin{array}[]{ll}-1,&\phi_{E}\text{ is odd}\\ 1,&\phi_{E}\text{ is even}.\end{array}\right.$$

$$H_{\mathrm{RWA}}(Q,P)=(F^2/6\gamma )g(Q,P),$$

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

$$E_n={\bar{E}}_n-{\bar{E}}_0,{\bar{E}}_n=\frac{1}{2}\hslash V{\left(n+\frac{1}{2}-\frac{\delta {\omega }_F}{V}\right)}^2.$$

$$\delta E_n=-\hslash V\frac{{\tilde{F}}^2}{4V^2}\frac{2{\bar{E}}_n/\hslash V-(\delta {\omega }_F/V{)}^2-3/4}{2{\bar{E}}_n/\hslash V-1}$$

$$H_{\mathrm{RWA}}=\frac{\hslash V}{2}\left(a^{+2}+\frac{\tilde{F}}{V}\right)\left(a^2+\frac{\tilde{F}}{V}\right)-\frac{\hslash {\tilde{F}}^2}{2V}.$$

$$i\hslash {\partial }_t\varphi (t)=H_{\mathrm{RWA}}(t)\varphi (t).$$

$$H_{\mathrm{RWA}}(t)=\hslash \left(\begin{array}{cc}\Delta & \nu (t)\\ \nu (t)& -\Delta \end{array}\right),$$

$$U_{\sigma }^{†}{H}_{\mathrm{RWA}}{U}_{\sigma }=H_{\mathrm{LZ}}=\hslash \left(\begin{array}{cc}\nu (t)& \Delta \\ \Delta & -\nu (t)\end{array}\right).$$

$$C_{\downarrow }(t)$$

$$C_{\uparrow }(t)$$

$$\theta (t)$$

$${\alpha }_{\uparrow }\approx 1-\frac{i}{12}\frac{s}{{\Delta }^2},{\alpha }_{\downarrow }\approx -\frac{i}{4}\frac{s}{{\Delta }^2}$$

$${\partial }_t\rho$$

$$\hat{\Gamma }\rho$$

$$W(\Omega ,t)=Q[\Omega ,t-t_0,{\rho }_0(t_0)]\Omega {\xi }^2(\Omega ),$$

$$Q[\Omega ,t-t_0,{\rho }_0(t_0)]=2\mathrm{Re}{\int }_{t_0}^{t}d{t}^{\prime }{e}^{i(\Omega -{\omega }_F/2)(t-t^{\prime })}$$

$$\times \mathrm{Tr}[a^{†}(t^{\prime }-t_0)a(t-t_0){\rho }_0(t_0)],$$

$$E_{\mathrm{rad}}(\Omega )={\int }_{t_0}^{\infty }dtQ[\Omega ,t-t_0,{\rho }_0(t_0)-{\rho }_{\mathrm{st}}].$$

$$\epsilon u_{\epsilon }(t)=(H(t)-i\hslash {\partial }_t)u_{\epsilon }(t).$$

$$\epsilon u_{k,n}=$$

$$M_{k,n}^{k^{\prime },n^{\prime }}=$$

$$+\frac{1}{4}F{q}_{n{n}^{\prime }}^2({\delta }_{k^{\prime },k+1}+{\delta }_{k^{\prime },k-1})$$

$$G_{k,2n^{\prime }}\equiv G_{n^{\prime },2n^{\prime }}$$

$$=\{u_{n^{\prime }+k^{\prime },2n^{\prime }+2k^{\prime }},k^{\prime }=-n^{\prime },-n^{\prime }+1,...\}=G_{0,0}.$$

$$G_{k,2n^{\prime }+1}\equiv G_{0,1}.$$

$$(\epsilon -{\mathcal{E}}_{2k}+k\hslash {\omega }_F)u_{k,2k}$$

$$=\tilde{F}\left[k(2k-1)u_{k-1,2k-2}+(k+1)(2k+1)u_{k+1,2k+2}\right],$$

$$(\epsilon -{\mathcal{E}}_{2k+1}+k\hslash {\omega }_F)u_{k,2k+1}$$

$$=\tilde{F}\left[k(2k+1)u_{k-1,2k-1}+(k+1)(2k+3)u_{k+1,2k+3}\right]$$

$$g\approx (\mu +1)(Q-Q_0{)}^2+P^2+g_{\min },$$

$$Q-Q_0=\sqrt{\frac{\lambda }{2}}(\mu +1{)}^{-1/4}(b^{†}+b),$$

Full text

I Introduction

Periodically driven quantum systems are described by quasienergy (Floquet) states, which are a time-domain analog of Bloch states in spatially periodic systems Shirley (1965); Zel’dovich (1967); Ritus (1967); Sambe (1973). The new physics associated with quasienergy states has been attracting much interest recently. Examples include topological Floquet states, artificial gauge fields, and new many-body phases Kitagawa et al. (2010); Lindner et al. (2011); Goldman et al. (2014); Bukov et al. (2015); Peano et al. (2016); Khemani et al. (2016); von Keyserlingk and Sondhi (2016); Khemani et al. (2016); Zhang et al. (2017); Choi et al. (2016); Bairey et al. (2017).

Preparation of Floquet states is often discussed in the adiabatic framework assuming that the periodic field is slowly turned on, cf. D’Alessio and Rigol (2015); Heinisch and Holthaus (2016); Weinberg et al. (2016); Ho and Abanin (2016) and references therein. The analysis for many-body systems is complicated by the effect of heating, and much progress has been made by studying systems that display many-body localization, as it may alleviate the heating. Recently, adiabatic state preparation was considered also for a parametrically driven nonlinear oscillator Goto (2016); Puri and Blais (2016). In contrast to many-body systems, the energy spectrum here is discrete, which simplifies the problem. However, a potential complication, and also potentially new and interesting features stem from the fact that the quasienergy levels for weak driving can display degeneracy, or a specific type of degeneracy, which we call the reduced-band (RB) degeneracy.

The goal of this paper is to study preparation of quasienergy states in a small quantum system in the case where the quasienergy states can display degeneracy or the RB degeneracy for weak driving. In optics terms, this case corresponds to either a multiphoton resonance or a subharmonic resonance, where the distance between the energy levels of the system is close to either a multiple or a fraction of the radiation frequency multiplied by $\hbar$. Multiphoton resonance leads to Rabi oscillations described in Ref. Larsen and Bloembergen, 1976 for a nonlinear oscillator using perturbation theory. In terms of the Floquet states, when the driving frequency is close to the oscillator eigenfrequency, such oscillator can display simultaneous multiple anticrossing of the quasienergy levels Dykman and Fistul (2005).

We will use as a model a driven quantum oscillator. Such model is interesting as it describes a broad range of physical systems, from molecular vibrations Larsen and Bloembergen (1976) to the modes of nonlinear optical and microwave cavities to Josephson junctions Dykman (2012). Here we study the features of the Floquet dynamics that emerge when an oscillator is driven parametrically and the drive frequency $\omega_{F}$ is close to twice the oscillator eigenfrequency.

To explain how the multiphoton and subharmonic resonances are seen in the quasienergy spectrum, we note that quasienergies of a system $\varepsilon_{n}$ and the quasienergy level spacing $\varepsilon_{n}-\varepsilon_{m}$ are defined modulo $\hbar\omega_{F}$. In the limit of zero driving $\varepsilon_{n}-\varepsilon_{m}$ is simply related to the spacing ${\cal E}_{n}-{\cal E}_{m}$ of the corresponding energy levels of the system, $\varepsilon_{n}-\varepsilon_{m}=({\cal E}_{n}-{\cal E}_{m})\!\!\!\mod(\hbar\omega_{F})$. The standard multiphoton resonance for weak driving occurs if ${\cal E}_{n}-{\cal E}_{m}$ is a multiple of $\hbar\omega_{F}$, and then $\varepsilon_{n}-\varepsilon_{m}=0$ for a given pair of states $(n,m)$, i.e., the quasienergies are degenerate. In contrast, in the case of a subharmonic resonance, ${\cal E}_{n}-{\cal E}_{m}$ can be a fraction of $\hbar\omega_{F}$. In particular, for the parametric resonance in an oscillator one can have $|{\cal E}_{n}-{\cal E}_{n+1}|=\hbar\omega_{F}/2$ (a more general resonant condition is discussed below, cf. Fig. 2). In this case $|\varepsilon_{n}-\varepsilon_{n+1}|=\hbar\omega_{F}/2$. This is the RB degeneracy, as the quasienergies would coincide if they were defined modulo $\hbar\omega_{F}/2$. Such degeneracy is nontrivial, since if the system is prepared in a superposition of the RB-degenerate states, it displays period doubling: the state is reproduced (up to a trivial phase factor) after twice the driving period, rather than after one period.

In what follows we show that, by slowly turning on resonant parametric drive, it is possible to prepare on demand various quasienergy states starting from the ground state of the oscillator ($n=0$). Importantly, this can be done in a finite time and with high accuracy without using special pulse-shaping techniques, but just by increasing the amplitude of the drive linearly in time. Such scenario is easy to implement in the experiment. We also study preparation of a superposition of quasienergy states starting from the ground state. Such preparation can be accomplished using non-adiabatic transitions for the driving frequency $\omega_{F}$ tuned close to multiphoton resonance, so that $\varepsilon_{m}-\varepsilon_{0}$ is small for the targeted $m$. Again, this relies on a simple linear increase of the driving amplitude. However, the nonadiabatic dynamics in this case turns out to be different from the conventional Landau-Zener dynamics.

The paper is organized as follows. In Sec. II, we present the model of a parametric nonlinear oscillator and discuss its quasienergy spectrum. We show the evolution of the spectrum with the varying driving frequency in the limit of zero drive amplitude and the occurrence of the degeneracy and the RB degeneracy of the quasienergy levels as the system goes through multiphoton or subharmonic resonance. In Sec. III, we present the Wigner distribution for the quasienergy states prepared from the oscillator ground state by slowly ramping up the amplitude of the driving in the absence of degeneracy. We demonstrate the possibility to prepare a Floquet state “on demand” and the rich structure of its Wigner distribution. The only constraint is that the resulting Floquet states are “even” with respect to inversion in phase space. In Sec. IV, we consider preparation of a superposition of two quasienergy states via a non-adiabatic transition when the system is close to degeneracy for weak field. In Sec. V, we briefly discuss the adiabaticity in the presence of dissipation. In Sec. VI we study fluorescence of the oscillator driven into a Floquet state, and in particular the characteristic transient spectrum of the fluorescence. Sec. VI contains concluding remarks.

II RWA Hamiltonian and quasienergy spectrum

The Hamiltonian of a weakly nonlinear parametric oscillator with coordinate $q$ and momentum $p$ has the form

[TABLE]

We assume that the driving amplitude $F$ and the nonlinearity are comparatively small, $F,\gamma\langle q^{2}\rangle\ll\omega_{0}^{2}$, and the driving frequency $\omega_{F}$ is close to resonance, $|\omega_{F}-2\omega_{0}|\ll\omega_{0}$; without loss of generality, we consider $F,\gamma>0$. A quantum parametric oscillator described by Eq. (1) has been realized in various platforms, from optical and microwave cavities to nanomechanical systems, cf. Nabors et al. (1990); Wilson et al. (2010); Dykman (2012); Lin et al. (2014).

For a periodically modulated quantum system, there exists a complete set of solutions to the Schrödinger equation called Floquet states, which are eigenfunctions of the operator $T_{t_{F}}$ of time translation by the modulation period ${}t_{F}$,

[TABLE]

Parameter $\varepsilon$ is called quasienergy or Floquet eigenvalue. For the parametric oscillator with Hamiltonian (1), ${}t_{F}=2\pi/\omega_{F}$.

A standard procedure to find quasienergy states and quasienergies is to plug the solution Eq. (2) into the Schrödinger equation, and then solve the resulting equation for $u_{\varepsilon}(t)$ using Fourier series expansion; see Appendix. For a driven oscillator, a much simpler way to find quasienergies is to go to the rotating frame at frequency $\omega_{F}/2$ by applying the standard unitary transformation $U(t)=\exp[-i\omega_{F}a^{\dagger}at/2]$, where $a$ and $a^{\dagger}$ are the oscillator ladder operators. In the rotating wave approximation (RWA) we disregard fast oscillating terms in the transformed Hamiltonian $U^{\dagger}HU-i\hbar U^{\dagger}\dot{U}$, which gives the RWA Hamiltonian

[TABLE]

where $\hat{n}=a^{\dagger}a$, $\delta\omega_{F}=\omega_{F}/2-\omega_{0}$ is the detuning frequency, $\tilde{F}=F/4\omega_{0}$, and $V=3\gamma\hbar/4\omega_{0}^{2}$.

The Hamiltonians $H$ and $H_{\rm RWA}$ commute with the parity operator $\hat{P}=\exp(-ia^{\dagger}a\pi)$ Haroche and Raimond (2006) that transforms $q\to-q,p\to-p$. Therefore, an eigenstate $\phi_{E}$ of $H_{\rm RWA}$ has definite parity $P_{E}=\pm 1$; here $E$ is an eigenvalue of $H_{\rm RWA}$, which can be called the RWA energy, $H_{\rm RWA}\phi_{E}=E\phi_{E}$. As a consequence, the corresponding time dependent state in the lab frame $\Phi_{E}(t)\equiv\exp(-iEt/\hbar)U(t)\phi_{E}$ is a Floquet state of Eq. (2). The quasienergy $\varepsilon$ and the periodic factor in the Floquet wave functions $u_{\varepsilon}$ are immediately expressed in terms of the RWA energy $E$ and the eigenfunction $\phi_{E}$,

[TABLE]

The aforementioned RB degeneracy where the quasienergies differ by $\hbar\omega_{F}/2$ occurs if $H_{\rm RWA}$ has degenerate states $\phi_{E}$. Such degeneracy is possible for a parametric oscillator for a finite driving amplitude Marthaler and Dykman (2007). A driven oscillator also provides a platform for investigating more complicated cases of RB degeneracy Zhang et al. (2017).

The understanding of the spectrum of $H_{\rm RWA}$ can be gained by looking at the Hamiltonian function $H_{\rm RWA}$ in the phase space of the oscillator in the rotating frame, i.e., by writing $H_{\rm RWA}$ in terms of the scaled quadratures $P$ and $Q$ defined as $Q=i(a-a^{\dagger})\sqrt{\lambda/2},P=(a^{\dagger}+a)\sqrt{\lambda/2}$. Here, $\lambda=V/2\tilde{F}$ is the dimensionless Planck constant. In these variables

[TABLE]

where $\mu=2\omega_{F}(\delta\omega_{F})/F$ Marthaler and Dykman (2007). The eigenstates of the Hamiltonian $H_{\rm RWA}$ can be written in the $Q$-basis, $\phi_{E}\equiv\phi_{E}(Q)$. The parity operator $\hat{P}$ is then the inversion operator, $\hat{P}\phi_{E}(Q)=\phi_{E}(-Q)$.

For $\mu+1>0$, function $H_{\rm RWA}(Q,P)$ has two minima located at $P=0,Q=\pm\sqrt{\mu+1}$. Function $H_{\rm RWA}(Q,P=0)$ is shown in Fig. 1. For sufficiently strong driving, where the two wells become deep and well-separated, the low-lying eigenstates of $H_{\rm RWA}$ are symmetric or anti-symmetric superpositions of intra-well states.

In the opposite limit of weak driving, $F\rightarrow 0$, the Hamiltonian $H_{\rm RWA}$, Eq. (3), is trivially diagonalized in the basis of the oscillator Fock states. What is interesting, however, is that the order of the RWA eigenstates in the rotating frame can be changed compared to the order of the Fock states in the laboratory frame. From Eq. (3), for $F=0$ the eigenvalues $E_{n}$ of $H_{\rm RWA}$ can be written in a suggestive form

[TABLE]

From Eq. (8), $E_{n}$ considered as a continuous function of $n$ is a simple parabola with a minimum at $n=\delta\omega_{F}/V-1/2$; see Fig. 2. For $\delta\omega_{F}/V<1/2$, $E_{n}$ quadratically increases with the increasing $n$; see the top line in Fig. 2. However, as the ratio $\delta\omega_{F}/V$ increases, $E_{n}$ bends over and has a minimum at some positive $n$. Of course, the actual RWA energies are determined by $E_{n}$ with integer $n$. When $\delta\omega_{F}/V>1$, the state with the lowest $E_{n}$ is no longer the Fock state $|0\rangle$. For instance, for $\delta\omega_{F}/V=1.8$ (blue dots, which lie on the second from top line in Fig. 2), this state is $|1\rangle$.

The reordering of the quasienergy states described by Eq. (8) is essential for preparing quasienergy states on demand. Indeed, if the oscillator is initially in the ground state, then by tuning the driving frequency and increasing the driving strength, we make this state an arbitrary even in $Q$ quasienergy state, i.e., an arbitrary superposition of Fock states $|m\rangle$ with even $m$. We also note that, for certain values of $\delta\omega_{F}/V$, there can be degenerate RWA levels (the green and brown dots, which lie on the two lowest curves in Fig. 2 and are connected by dashed lines). We will discuss such degeneracy later in details.

The driving mixes Fock states with the same parity. The evolution of the RWA spectrum with the increasing $F$ is shown in Fig. 3 for different values of the detuning $\delta\omega_{F}/V$. A common trend is that RWA energy levels of the same parity repel each other, whereas neighboring levels of opposite parity attract each other and form pairs for large $\tilde{F}/V$. As mentioned above, such pairs for large $\tilde{F}/V$ are even and odd superposition of “intra-well” states of $H_{\rm RWA}$. The distance between the states within the pairs is determined by interwell tunneling Marthaler and Dykman (2007).

II.1 Special features of the RWA spectrum

We find that, somewhat counterintuitively, the RWA levels do not cross each other as $F$ changes. Therefore any gaps that are present at $F\rightarrow 0$ will remain open for any finite $F$. For instance, Figs. 3a and b refer to the cases where the Fock state $|0\rangle$ is the first and the third lowest RWA eigenstate at $F\rightarrow 0$, respectively. As $F$ increases, it remains the first and the third lowest RWA eigenstate. Such non-crossing feature will be important for the preparation of quasienergy states by slowly turning on the driving.

A remarkable feature of the RWA spectrum is that, when the ratio $\delta\omega_{F}/V$ is a positive integer, there is a set of simultaneously doubly-degenerate levels of opposite parity regardless of the value of $F$. For $F\rightarrow 0$, this can be readily seen from Eq. (8) (cf. Dykman and Fistul (2005) where a similar feature was found in the case of the driving at frequency close to $\omega_{0}$). When $\delta\omega_{F}/V=k,k=1,2,3..$, the minimum of $E_{n}$ as a continuous function of $n$ is reached at half odd integer $n=k-1/2$. Since $E_{n}$ is a symmetric function of $n$ with respect to the minimum, the levels separated by $\Delta n=2m+1$ are degenerate, that is, $E_{k+m}=E_{k-(m+1)}$, for $m=0,1..,k-1$. The green curve in Fig. 2 (the third from the top) refers to the case $k=2$, where the degeneracy condition is met.

The degeneracy of the RWA energy levels persists for nonzero $F$, as can been seen in Fig. 3c. At weak driving, this follows from the perturbation theory. To the second order in $F$, the correction $\delta E_{n}$ to $E_{n}$ is

[TABLE]

The dependence of $\delta E_{n}$ on the level number $n$ is exactly the same as that of $E_{n}$. cf. Eq. (8). Therefore, if $E_{n}=E_{n^{\prime}}$, then $\delta E_{n}=\delta E_{n^{\prime}}$. Note that the perturbation theory still applies even if there are degenerate levels of opposite parity since there is no coupling between them. At strong driving, such degeneracy corresponds to the vanishing of tunnel splitting found in Ref. Marthaler and Dykman (2007).

For the special case $\delta\omega_{F}/V=1$, $H_{\rm RWA}$ can be factored Puri and Blais (2016),

[TABLE]

In this case the coherent states $|\pm\alpha\rangle$, $\alpha=\sqrt{-\tilde{F}/V}$, are exact degenerate eigenstates of $H_{\rm RWA}$ for arbitrary driving strength. However, no such eigenstates are known for other values of $\delta\omega_{F}/V$.

If the ratio $\delta\omega_{F}/V$ is a half-integer, $\delta\omega_{F}/V=(2k+1)/2,k=1,2,3,...$, the minimum of function $E_{n}$ for $F\rightarrow 0$ is reached at integer $n=k$. Again, due to the parabolic dependence of $E_{n}$ on $n$, levels $E_{k\pm m}$ are degenerate for $m=1,2,...,k$. For instance, the lowest (brown) curve in Fig. 2 refers to the case $k=2$. The degeneracy of the levels of the same parity occurs when the driving frequency equals to one of the transition frequencies of the undriven oscillator. This can be seen by rewriting $E_{n}$ as $E_{n}=-n\hbar\omega_{F}/2+\mathcal{E}_{n}$, where $\mathcal{E}_{n}=n\hbar\omega_{0}+\hbar Vn(n+1)/2$ is the $n$th energy level of the oscillator in the absence of driving. Clearly, the degeneracy condition $E_{k+m}=E_{k-m}$ is equivalent to $\mathcal{E}_{k+m}-\mathcal{E}_{k-m}=m\omega_{F}$, which is the m-photon resonance condition for transition from $\mathcal{E}_{k-m}$ to $\mathcal{E}_{k+m}$. The degeneracy is lifted at finite $F$ due to level repulsion, as shown in Fig. 3d.

III Effectively-adiabatic

preparation of quasienergy states and the Wigner distribution

The observation that the quasienergy levels of the same parity do not approach each other with the increasing field $F$ is critical for state preparation. It allows one to prepare a quasienergy state by slowly turning on the field, provided the states are non-degenerate for $F\to 0$.

We consider ramping up the driving amplitude $\tilde{F}$ linearly with speed $s$ starting at $t=0$, $\tilde{F}(t)=s_{0}t$. If $\sqrt{s_{0}}$ is small compared to $\omega_{0}$, the time evolution of the oscillator wave function $\phi(t)$ can be described in the RWA,

[TABLE]

We will solve this equation assuming that initially, for zero driving, the system is in the ground state of the oscillator, $\phi(Q,t=0)=|0\rangle$.

The results of the numerical solution of Eq. (10) are illustrated in Fig. 4. The values of $\delta\omega_{F}/V$ were chosen in such a way that, in one case ($\delta\omega_{F}=0$), the state of the system remains close to the eigenstate of $H_{\rm RWA}$ with the lowest eigenvalue $E_{n}$, whereas in the other case ($\delta\omega_{F}/V=1.8$) it is close to the third lowest-$E_{n}$ state, cf. Fig. 3(b). The quality of the adiabatic approximation for the chosen parameters can be characterized by the inner product of the state $\phi(Q)$ at the end of ramp-up and the corresponding stationary RWA eigenstate $\phi_{E}(Q)$ calculated for $\tilde{F}=\tilde{F}_{\rm final}$. This inner product is 0.997 and 0.98 for the cases shown in Fig. 4a and Fig. 4b, respectively, which shows that the adiabatic approximation is very good.

The final value of the field amplitude $\tilde{F}_{\rm final}$ in Fig. 4 refers to the case where the Hamiltonian function $H_{\rm RWA}(Q,P)$, Eq. (II), has a pronounced double-well structure, cf. Fig. 1. For $\delta\omega_{F}=0$, the state $\phi(Q)$ is well described by a symmetric superposition of the lowest intra-well states in Fig. 1, $\phi(Q)=(\phi_{L}+\phi_{R})/\sqrt{2}$ where $\phi_{L}$ and $\phi_{R}$ refer to the left and right well, respectively. Near their maxima, functions $\phi_{L,R}$ are given by squeezed coherent states with equal amplitude and opposite phases, $\phi_{L,R}\propto\exp[-(Q\pm Q_{0})^{2}/2\lambda\eta]$ where $Q_{0}=\sqrt{\mu+1}$ is the position of the right well and $\eta=1/\sqrt{\mu+1}$ characterizes the state squeezing, see Appendix B. The adiabatic preparation of such “cat” state has been discussed in Refs. (Goto (2016); Puri and Blais (2016)).

In contrast, for the case in Fig. 4b, the driving brings the system to an excited state of $H_{\rm RWA}$. The state $\phi(Q)$ for $t=\tilde{F}_{\rm{final}}/s_{0}$ is no longer a superposition of the lowest intra-well states but, for the chosen $\delta\omega_{F}/V$, the superposition of the second lowest intra-well states, $\phi(Q)=(\phi_{L}^{\prime}+\phi_{R}^{\prime})/\sqrt{2}$. Near their maxima, functions $\phi_{L,R}^{\prime}$ are well described by a displaced and squeezed Fock state $|1\rangle$: $\phi_{L,R}^{\prime}\propto(Q\pm Q_{0})\exp[-(Q\pm Q_{0})^{2}/2\lambda\eta].$ Since the RWA energy levels for small $F$ in this case are closer than for $\delta\omega_{F}=0$, in particular the Fock states $|0\rangle$ and $|2\rangle$ have close RWA energies, we had to use a much slower increase of the driving amplitude to attain high fidelity of the prepared large-$F$ state.

IV Preparing a superposition of quasienergy states nonadiabatically

As the driving amplitude $F$ is ramped up, the non-adiabaticity can mix quasienergy states of the same parity. The mixing is particularly strong if the quasienergy gap that separates the states is small. As shown in Sec. II, this gap is controlled by the driving frequency. In this section, we consider a situation where two nearest quasienergy states of the same parity have close quasienergies for $F\to 0$, whereas the quasienergies of other states are significantly different, so that mixing with these other states can be disregarded for slowly varying $F(t)$. We show that, by ramping up the driving amplitude linearly in time, we can prepare, with high accuracy, a desired coherent superposition of the chosen two quasienergy states.

We assume that the states with close quasienergies for $F\to 0$ are $|n-1\rangle$ and $|n+1\rangle$, which means that $\delta\omega_{F}/V\approx n+1/2$. As the drive is ramped up, these states are mixed with each other. Concurrently, they are mixed with other states of the same parity. However, this mixing is nonresonant and therefore is weaker.

The picture of the state evolution is as follows. The resonant mixing leads to a redistribution of the initial population between the resonating states and to a separation of their quasienergies already for a comparatively weak field, see Fig 5. The increase of the field afterwards does not change the state populations, even though it modifies the states by increasingly strongly admixing them to other states of the same parity.

To describe the initial stage of the evolution we project the Hamiltonian $H_{\rm RWA}$ onto the subspace formed by the states $|n-1\rangle$ and $|n+1\rangle$, subtract the mean RWA energy $(E_{n+1}+E_{n-1})/2$, and disregard the coupling to other states. Then the Hamiltonian becomes

[TABLE]

where $\Delta=(E_{n-1}-E_{n+1})/2\hbar$, $\nu(t)=\sqrt{n(n+1)}\tilde{F}(t)$. For a field that linearly increases in time $\nu(t)=s{}t$.

It is convenient to re-write the Hamiltonian (11) in the conventional form used in the analysis of the Landau-Zener tunneling. Making a unitary transformation $U_{\sigma}=(1/\sqrt{2})(\sigma_{z}+\sigma_{x})$ ($\sigma_{x,z}$ are Pauli matrices), we obtain

[TABLE]

Note that the vectors $\left(\begin{matrix}1\\ 0\end{matrix}\right)$ and $\left(\begin{matrix}0\\ 1\end{matrix}\right)$ for the Hamiltonian (12) are, respectively, the wave functions $(|n-1\rangle+|n+1\rangle)/\sqrt{2}$ and $(|n-1\rangle-|n+1\rangle)\sqrt{2}$.

The only difference of the evolution of the states we consider here from the standard Landau-Zener scenario is that the initial condition for the Schrödinger equation $i\hbar\dot{\phi}(t)=H_{\rm LZ}\phi(t)$ is set for $t=0$ and the problem is considered on the semi-axis $t\geq 0$. It is convenient to seek the wave function as $\phi(t)=(1/\sqrt{2})\sum_{\alpha=\pm}C_{\alpha}(t)[|n-1\rangle+\alpha|n+1\rangle]$. We will be interested in the solution that corresponds to the initial condition where the smaller-$n$ state is occupied while the larger-$n$ state is empty, $C_{+}(0)=C_{-}(0)=1/\sqrt{2}$. As in the Landau-Zener problem, the solution to the Schrödinger equation can be expressed in terms of the parabolic cylinder functions; see Appendix C.

In Fig. 6, we show the result for the coefficient $C_{\uparrow}(t)$, which is equal to the projection $\langle\phi^{\rm ad}_{\uparrow}(t)|\phi(t)\rangle$ of the wave function $\phi(t)$ on the upper branch (the higher energy branch in Fig. 5) of the adiabatic solutions $\phi_{\uparrow,\downarrow}(t)$ of the Schrödinger equation, $H_{\rm LZ}\phi^{\rm ad}_{\uparrow,\downarrow}(t)=\pm[\nu^{2}(t)+\Delta^{2}]^{1/2}\phi^{\rm ad}_{\uparrow,\downarrow}(t)$. The result is in full agreement with the numerical solution of the Schrödinger equation.

Of primary interest is the long time behavior of $C_{\uparrow,\downarrow}(t)$. It can be obtained from the asymptotic expansion of the parabolic cylinder functions (see Appendix C), or directly by solving the Schrödinger equation in the WKB approximation,

[TABLE]

Here, $\theta(t)$ is the dynamical phase $\int_{0}^{t}dt^{\prime}\sqrt{\nu^{2}(t^{\prime})+\Delta^{2}}$ associated with the adiabatic solutions in Fig. 5,

[TABLE]

The expressions for the parameters $\alpha_{\uparrow,\downarrow},\beta_{\uparrow,\downarrow}$ in Eq. (13) follow from the general solution of the Schrödinger equation; the explicit form of $\alpha_{\uparrow,\downarrow}$ is given in Appendix C.

The coefficients $C_{\uparrow,\downarrow}(t)$ approach their asymptotic values $\propto\alpha_{\uparrow,\downarrow}$ as $1/t$ and oscillate as $\exp[\pm i\theta(t)]$. We note that, for $t\rightarrow\infty$, we have $C_{\uparrow}\to C_{+}(t)$ and $C_{\downarrow}\to C_{-}(t)$, i.e., Eq. (13) directly gives the coefficients $C_{\pm}$ of the expansion of the wave function in the symmetric and antisymmetric combination of functions $|n\pm 1\rangle$.

Figure 7 shows the asymptotic value $|C_{\uparrow}(\infty)|^{2}=|\alpha_{\uparrow}|^{2}$ as a function of the Landau-Zener parameter $\Delta^{2}/s{}$. In the adiabatic limit $\Delta^{2}/s{}\gg 1$ and for the case $\Delta>0$, where the system starts from the upper branch, $\phi(0)=\phi^{\rm ad}_{\uparrow}(0)$, we have

[TABLE]

($|\alpha_{\uparrow}|^{2}+|\alpha_{\downarrow}|^{2}=1$). In distinction from the Landau-Zener problem, where the non-adiabatic transition probability approaches zero exponentially as $\exp(-\pi\Delta^{2}/s{})$, here it approaches zero as $(\Delta^{2}/s)^{-2}$. This special feature is due to the initial condition in the considered problem being set at $t=0$ rather than $t\to-\infty$.

In the strongly non-adiabatic case, $\Delta^{2}/s\ll 1$, if $\phi(0)=|n-1\rangle$, in the long-time limit the state of the system ultimately approaches an equal superposition of the eigenstates $\phi^{\rm ad}_{\uparrow,\downarrow}$ of $H_{\rm LZ}$: $|\alpha_{\uparrow}|\approx|\alpha_{\downarrow}|\approx 1/\sqrt{2}.$ This can be seen from Eq. (12); see also Appendix C. In the case $\Delta=0$, the states $(|n-1\rangle\pm|n+1\rangle)/\sqrt{2}$ are exact eigenstates for any time $t$. Therefore, $\phi(t)$ will remain in an equal superposition of these two states for any time; note, however, that the states depend on time differently.

An instructive case is when the oscillator is in the ground state before the driving is applied and the detuning of the driving frequency $\delta\omega_{F}$ is close to $3V/2$. Here, if the field is ramped up fast, the oscillator will end up in equally populated adiabatic states, which corresponds to two equally populated even interwell states in Fig. 1.

V Adiabaticity in the presence of dissipation

Coupling to the environment leads to decoherence of the quasienergy states. It reduces the fidelity of the state preparation. Here we consider the constraint on the dissipation in the case of state preparation by slowly ramping up the driving field. To achieve high fidelity, one needs to increase the field at a rate larger than the relaxation rate, but smaller than the reciprocal spacing of the relevant RWA energies divided by $\hbar$. For a state $\phi_{E}$, this means that the decay rate of this state $\Gamma_{E}$ should be small compared to $\Delta_{E}$, where $\hbar\Delta_{E}$ is the instantaneous difference between the quasienergy of the state $\phi_{E}$ and the nearest state of the same parity. The parity constraint here is the consequence of the fact that the field mixes only the same-parity states.

The RWA level spacing $\hbar\Delta_{E}$ can be estimated where the driving is weak, $\tilde{F}\ll V$, or strong, $\tilde{F}\gg V$. For weak driving, the RWA eigenstates are close to the Fock states. From the results of Sec. II, $\Delta_{E}\sim V$ and depends on the ratio $\delta\omega_{F}/V$, cf. Fig. 2. At strong driving, $\hbar\Delta_{E}$ is given by the spacing of the intrawell energy levels of the Hamiltonian $H_{\rm RWA}(Q,P)$; see Fig. 1. It is determined by the frequency $\omega_{\rm min}$ of oscillations about the minima of $H_{\rm RWA}(Q,P)$, which gives $\Delta_{E}\approx 2[(\delta\omega_{F}+\tilde{F})\tilde{F}]^{1/2}$; see Appendix B.

We illustrate the effect of dissipation using the well-known model Mandel and Wolf (1995) where the kinetics in the rotating frame is described by the Markov master equation for the density matrix $\rho$ of the form

[TABLE]

Here, $\Gamma$ is the oscillator relaxation rate and we assume that the temperature of the environment is sufficiently low, $k_{B}T\ll\hbar\omega_{0}$.

The decay rate $\Gamma_{E}$ of an RWA eigenstate $\phi_{E}$ can be estimated as the decay rate of the diagonal matrix element of the density matrix $\langle\phi_{E}|\rho|\phi_{E}\rangle$. Assuming that the system is in state $\phi_{E}$, i.e., $\rho=|\phi_{E}\rangle\langle\phi_{E}|$, and taking into account that the matrix elements of the ladder operators on the states of the same parity are zero, we find from Eq. (16) $\Gamma_{E}=2\Gamma\langle\phi_{E}|a^{\dagger}a|\phi_{E}\rangle$. At weak driving, $\Gamma_{E}\sim\Gamma$. At strong driving $\Gamma_{E}$ is determined by the rate of transitions between the intrawell states of $H_{\rm RWA}$ Marthaler and Dykman (2006), $\Gamma_{E}\sim\Gamma\tilde{F}/V$.

From the above estimates, the adiabaticity condition $\Gamma_{E}\ll\Delta_{E}$ requires that $\Gamma\ll V,|\delta\omega_{F}|$ for weak driving and $\Gamma\ll V$ for strong driving. Fig. 8 illustrates the evolution of $\Delta_{E}$ and $\Gamma_{E}$ of an RWA eigenstate $\phi_{E}$ with the varying driving amplitude $F$. For the case shown in the figure, the state $\phi_{E}$ has the lowest RWA eigenenergy. At large $\tilde{F}/V$, both $\Delta_{E}$ and $\Gamma_{E}$ increase linearly with $F$ as we expect from the analysis above. The slope of $\Gamma_{E}$ as a function of $\tilde{F}$ increases as $\Gamma/V$ increases. It coincides with the slope of $\Delta_{E}$ for $\Gamma/V\approx 2$ as shown by the green curve. For the condition $\Gamma_{E}\ll\Delta_{E}$ to be satisfied for any $F$, one needs to have $\Gamma/V\ll 2$. For $\Gamma/V\gtrsim 2$, in the considered case $\Gamma_{E}$ and $\Delta_{E}$ as a function of $\tilde{F}$ can cross each other.

VI Transient radiation from quasienergy states

Decay of a parametrically driven oscillator is accompanied by emission of excitations into the surrounding medium. The most familiar picture is decay of optical/microwave cavity modes into propagating electromagnetic waves. Detection of the radiation from the cavity provides a way of characterizing the cavity modes. Radiation from the modes in a non-steady state, such as a quasienergy state, is transient. After a time of the order of the mode relaxation time, the system relaxes to a steady state, the radiation becomes steady and does not depend on the quasienergy state the system had been staying in. To identify a quasienergy state from the radiation, one needs to collect the transient radiation.

We model the radiation field by a set of oscillators enumerated by subscript $k$, with quasi-continuous frequencies $\omega_{k}$ and with Hamiltonian $H_{\rm rad}=\sum_{k}\hbar\omega_{k}b^{\dagger}_{k}b_{k}$. We assume that the coupling of the considered oscillator to this field is bilinear in the ladder operators of the oscillator and the radiation, $H_{i}=\sum_{k}\xi_{k}(b_{k}+b_{k}^{\dagger})(a+a^{\dagger})$, where $\xi_{k}$ are the coupling parameters. The total Hamiltonian is $H_{\rm total}=H_{0}+H_{\rm rad}+H_{i}$. Operator $H_{0}$ is the Hamiltonian of the oscillator and the non-radiative thermal reservoir to which the oscillator is coupled. We assume that this reservoir and the radiation field are at the same temperature, which we assume to be sufficiently low, $k_{B}T\ll\hbar\omega_{0}$. The coupling to the reservoir leads to relaxation of the oscillator with typical relaxation rate $\Gamma$, cf. Eq. (16).111The oscillator is characterized also by a much longer rate, which is related to the dissipation-induced transitions between the wells in Fig. 1.

If the coupling to the radiation field is weak, it can be considered as a perturbation to the non-radiative dynamics. The power of the radiation emitted into a spectral range $d\Omega$ around frequency $\Omega$ is given by the change of the energy of the radiation field in this interval per unit time $W(\Omega,t)d\Omega=\frac{d}{dt}\sum_{k}\delta(\omega_{k}-\Omega)d\Omega\langle\hbar\omega_{k}b_{k}^{\dagger}b_{k}\rangle$. To the lowest order in the coupling strength $\xi_{k}$, we have in the resonant region where $\Omega$ is close to $\omega_{F}/2$ Dykman (1975)

[TABLE]

where $\xi^{2}(\Omega)=\hbar^{-1}\sum_{k}|\xi_{k}|^{2}\delta(\Omega-\omega_{k})$.

In deriving Eq. (17) we assumed that the coupling to the radiation is switched on at time $t_{0}$; $\rho_{0}(t_{0})$ is the density matrix of the oscillator and the non-radiative environment. Equation (17) is written in the rotating frame used above to find the quasienergy states of the oscillator, with the time counted off from $t_{0}$.

The two-time correlation function in Eq. (17) can be found by solving the quantum kinetic equation. As an initial condition to this equation we choose the density matrix $\rho_{0}(t_{0})$ in the form of a product of the oscillator density matrix $\rho(t_{0})$ and the density matrix of the non-radiative environment in thermal equilibrium. Such choice is justified, since a weak coupling to the dissipative (non-radiative) environment allows preparing the oscillator in a certain state at time $t_{0}$ given that the preparation time is short compared to the relaxation time. The following evolution on the time scale, which largely exceeds both $\omega_{F}^{-1}$ and the time it took to prepare the state, can be described by assuming that at $t_{0}$ there is switched on not only the coupling to the radiation field, but also the stronger (but still weak) coupling to the non-radiative environment. Corrections to the dynamics due to the switching are well-understood, they are small in the considered case Dykman (1975).

The time evolution of the oscillator density matrix in the rotating frame is then often described by Eq. (16). To study transient radiation, we set $\rho(t_{0})=|\phi_{E}\rangle\langle\phi_{E}\rangle$, where $\phi_{E}$ is a RWA eigenstate in which the oscillator is prepared at $t_{0}$.

For not very strong driving, $\tilde{F}\lesssim V$, function $\phi_{E}$ has a contribution of only a few Fock states. Respectively, the oscillator will radiate only a few photons as it comes to the stationary state. Then, rather than measuring the radiation power $W(\Omega,t)$ it is more feasible to measure the total energy emitted over the transient time. The observation time should exceed the relaxation time to enable sufficient spectral resolution.

The energy of the transient radiation has to be separated from the energy that the oscillator emits in the stationary state. This can be done by noting that the latter energy is proportional to the observation time. The spectral power density (power per unit frequency) in the stationary regime is given by Eq. (17) written for $t\to\infty$ Dykman (2012). Therefore one can define the transient radiation spectral density as the integral over time of the difference of the emitted power (17) and the power emitted in the stationary regime. Writing this spectral density as $\Omega\xi^{2}(\Omega)E_{\rm rad}(\Omega)$, we obtain

[TABLE]

Here, $\rho_{\rm st}$ is the stationary density matrix of the driven oscillator and the non-radiative environment.

The spectral density $E_{\rm rad}(\Omega)$ is given by the difference between the irradiated energy and the energy that would be irradiated into the same spectral interval if the system were stationary. This difference is accumulated over a sufficiently long time that largely exceeds the relaxation time. By construction, it can be positive or negative.

As the oscillator decays from the initial state $\phi_{E}$, it emits radiation at frequencies $\omega_{F}/2+(E-E^{\prime})/\hbar$, where $E^{\prime}$ is the RWA energy of a state $\phi_{E^{\prime}}$ into which the oscillator can make a dipolar transition from $\phi_{E}$. In contrast, in the stationary state, the oscillator generally can be found in the both states $\phi_{E},\phi_{E^{\prime}}$, with different probabilities. Depending on these probabilities, it radiates at the both frequencies $\omega_{F}/2\pm(E-E^{\prime})/\hbar$ generally with different intensities. As a result, in the spectrum $E_{\rm rad}(\Omega)$ one may expect a peak or a dip at $\omega_{F}/2+(E-E^{\prime})/\hbar$, but only a dip at $\omega_{F}/2-(E-E^{\prime})/\hbar$.

Figures 9 (a) and (b) show the spectrum $E_{\rm rad}(\Omega)$ when the oscillator is initially in the RWA eigenstate $\phi_{E}$ prepared from the vacuum $|0\rangle$ by adiabatically ramping up the driving field. The driving frequency is chosen so that $\phi_{E}$ has the second lowest RWA energy among even states; see the insets. The transient radiation is dominated by transitions from the state $\phi_{E}$ to the lowest odd state $\phi_{E^{\prime}}$. In this case, for a strong driving field $E-E^{\prime}\approx 2\hbar[\tilde{F}(\delta\omega_{F}+\tilde{F})]^{1/2}$ corresponds to the spacing between the two lowest intrawell states of $H_{\rm RWA}$ in Fig. 1. For weak driving, $E-E^{\prime}\approx\hbar(\delta\omega_{F}-V)$; the frequency $\omega_{F}/2-(E-E^{\prime})/\hbar=\omega_{0}+V$ is the frequency of the transition from the first excited state to the ground state of the undriven oscillator. Figures 9 (a) and (b) refer not to these limiting cases but to the intermediate field strengths.

As expected, the spectrum $E_{\rm rad}(\Omega)$ displays a peak at $\omega_{F}/2+(E-E^{\prime})/\hbar$ for relatively strong driving and a small dip at this frequency for weak driving. It also displays a characteristic pronounced dip at $\omega_{F}/2-(E-E^{\prime})/\hbar$ in the both cases. In addition, for strong driving, the spectrum has a negative narrow peak at $\omega_{F}/2$ due to the interwell transitions Dykman (2012). For a comparison, Figs. 9 (c) and (d) show the steady-state radiation power spectrum $Q_{\rm st}(\Omega)=Q(\Omega,\infty,\rho_{\rm st})$ for the same parameters as in Figure 9(a) and (b), respectively.

VII Conclusions

We have studied preparation of quasienergy states of a nonlinear oscillator. We found that various states can be prepared with high accuracy in a finite time by simply linearly increasing in time the amplitude of the parametric driving. The driving frequency $\omega_{F}$ was chosen to be close to twice the oscillator eigenfrequency $\omega_{0}$, so that strong excitation of the oscillator could be achieved for a comparatively weak driving field. The prepared state sensitively depends on the interrelation between the detuning of the driving frequency $\delta\omega_{F}=\omega_{F}/2-\omega_{0}$ and the nonequidistance $V$ of the oscillator energy levels due to the nonlinearity (in frequency units).

An important factor for the state preparation is that the quasienergy states are either even or odd with respect to inversion in the phase space. The states of different parity are not coupled by the driving. This allows one to prepare on demand an arbitrary even quasienergy state just by slowly ramping up the driving, if the oscillator is initially in the ground state. The resulting states have very different structures in phase space, as evidenced by the Wigner tomography. A similar analysis shows that an arbitrary odd state can be prepared, if initially the oscillator is in the first excited state.

A remarkable feature of the system related to its symmetry is that the oscillator energy levels calculated in the rotating wave approximation do not cross or anti-cross with the increasing driving amplitude. Rather the neighboring RWA energy levels of even and odd states approach each other pairwise. At the same time, the levels of the opposite-parity states can cross with varying $\delta\omega_{F}$. This crossing does not lead to crossing of the quasienergy levels. Where the RWA energy levels cross, the quasienergy levels are separated by $\hbar\omega_{F}/2$.

It is also important for the state preparation that, in the limit of zero driving, the RWA energy spectrum can simultaneously have several double-degenerate levels. Such degeneracy corresponds to either a multi-photon or a subharmonic resonance. By tuning the driving frequency, one can bring the RWA energy levels closer or further away from the pairwise degeneracy.

The degeneracy of same-parity states provides an effective way of preparing superpositions of quasienergy states. It is based on non-adiabatic transitions induced by the increasing driving amplitude. The field leads to the state mixing that depends on how fast it is increased. The problem differs from the standard Landau-Zener problem, since the initial state is close to degeneracy and the field is ramped up in a finite time. As a result, for a linearly increasing field, the probability of the non-adiabatic transition falls off as a power law, rather than exponentially, with the Landau-Zener parameter $\Delta^{2}/s$, where $\Delta$ is the level spacing and $s$ is the ramp-up speed.

Dissipation due to the coupling to a thermal reservoir reduces the fidelity of the state preparation. However, away from the level degeneracy, the effect of the dissipation is small, if the oscillator nonlinearity parameter $V$ exceeds the decay rate $\Gamma$. Then one can ramp up the driving at a rate that is much smaller than the quasienergy level spacing, yet much larger than the decay rate. Fluctuations of the system parameters and of the driving power can also reduce the fidelity. Their effect is small if their bandwidth is small compared to $V$ or if they are sufficiently weak, so that their effect does not accumulate over the duration of the state preparation.

Because of dissipation, the parametric oscillator prepared in a given quasienergy state will ultimately come to a stationary state. Our results show that the prepared state can be characterized by studying the transient radiation of the oscillator. This method is complimentary to the commonly used Wigner tomography. It can be particularly useful for investigating quasienergy states of cavity modes in microwave cavities, the area of much current interest. The above analysis suggests a simple way of preparing various quasienergy states in such cavities as well as in other systems, for example, Josephson junctions, that can be modeled by nonlinear quantum parametric oscillators.

VIII Acknowledgements

This work was supported in part by the National Science Foundation (Grant No. DMR-1514591); YZ was also partly supported by the U.S. Army Research Office (W911NF1410011) and by the National Science Foundation (DMR-1609326).

Appendix A Fourier series for quasienergy states

The eigenvalue problem for the periodic part $u_{\varepsilon}(t)$ of the Floquet wave function defined in Eq. (2) reads

[TABLE]

Since $u_{\varepsilon}(t)$ and $H(t)$ are both periodic in time, it is convenient to expand them in Fourier series. It is also convenient to write $u_{\varepsilon}(t)$ in the basis of the Fock states $|n\rangle$ of the harmonic oscillator with frequency $\omega_{0}$. Then $u_{\varepsilon}(t)=\sum_{k,n}u_{k,n}\exp(-ik\omega_{F}t)|n\rangle$ and Eq. (19) takes the form of the standard eigenvalue problem

[TABLE]

where $q^{2}_{nm}=\langle n|q^{2}|m\rangle$ and $\mathcal{E}_{n}$ is the $n$th energy level of the Duffing oscillator in the absence of driving; to the leading order in the nonlinearity $\mathcal{E}_{n}=\hbar[\omega_{0}n+V(n^{2}+n)/2]$. The sum runs over $k=0,\pm 1,\pm 2,...$ and $n=0,1,2,...$.

The matrix elements $q^{2}_{nm}$ are nonzero for $n=m$ and $n=m\pm 2$. Therefore the driving term $\propto F$ couples $u_{k,n}$ to $u_{k\pm 1,n\pm 2},u_{k\pm 1,n}$. However, only the coupling to $u_{k+1,n+2}$ and $u_{k-1,n-2}$ is resonant, since the diagonal elements of matrix $\hat{M}$ for such $u$ are close; for example, $(\mathcal{E}_{n}-k\hbar\omega_{F})-[\mathcal{E}_{n+2}-(k+1)\hbar\omega_{F}]=2\hbar\delta\omega_{F}-\hbar V(2n+3)$ is small compared to $\hbar\omega_{F}$. Therefore, one can limit the analysis to a set $G_{k,n}$ of the variables $u_{k^{\prime},n^{\prime}}$ resonantly coupled to $u_{k,n}$. It has the form $G_{k,n}=\{u_{k+k^{\prime},n+2k^{\prime}},k^{\prime}\in\mathbb{Z}\text{ and }k^{\prime}\geq-n/2\}$. This is the rotating wave approximation in the Floquet formulation (19).

The sets $G_{k,n}$ with different $k$ but the same $n$ are equivalent: indeed, changing $k\to k_{1}$ corresponds to changing $\varepsilon\to\varepsilon+(k-k_{1})\hbar\omega_{F}$ in Eq. (20). Since $\varepsilon$ is defined modulo $\hbar\omega_{F}$, such change makes no difference. We can then simplify $G_{k,n}$ as follows. Consider first even $n$, i.e., $n=2n^{\prime}$, and set $k=n^{\prime}$,

[TABLE]

In the last equation, we simply redefined $k^{\prime}$ to absorb $n^{\prime}$ in the new definition.

Similarly, for odd $n$, where $n=2n^{\prime}+1$,

[TABLE]

The simplification described by Eqs. (21) and (22) allows one to reduce Eq. (20) to two sets of equations,

[TABLE]

where $\tilde{F}=F/4\omega_{0}\approx F/2\omega_{F}$ and we used the explicit form of the matrix elements $\langle n|q^{2}|n\pm 2\rangle$.

Equation (23) coincides with the RWA Schrödinger equation $E\phi_{E}=H_{\rm RWA}\phi_{E}$ if one writes $\phi_{E}$ in the basis of the Fock states and replaces $\varepsilon$ with $E$ using Eq. (6).

Appendix B Semiclassical analysis of RWA Hamiltonian

For completeness, here we present, following Marthaler and Dykman (2006), the description of the scaled RWA Hamiltonian function $g(Q,P)$ for large driving. For $\mu<-1$, function $g$ has one minimum at $(Q,P)=(0,0)$. For $-1<\mu<1$, the minimum at (0,0) becomes a saddle point and there appears two minima located at $(Q,P)=(\pm Q_{0},0),Q_{0}=\sqrt{\mu+1}$. For $\mu>1$, the saddle point at $(0,0)$ becomes a minimum again and there appear two saddle points at $(Q,P)=(0,\pm\sqrt{\mu-1})$.

Of primary interest in this paper is the regime $0\leq\mu<1$ where the quasinergy spectrum can display degeneracy and RB degeneracy. We expand $g$ about the minimum at $(Q_{0},0)$ to second order in $Q-Q_{0}$ and $P$,

[TABLE]

where $g_{\rm min}=-(\mu+1)^{2}/4$.

Introducing ladder operators $b,b^{\dagger}$ defined as

[TABLE]

[TABLE]

($[b,b^{\dagger}]=1$), we write the Hamiltonian $g(Q,-i\lambda\partial_{Q})$ for low-lying intrawell eigenstates in the form

[TABLE]

The eigenstates of operator $b^{\dagger}b$ give the intra-well states used in the main text.

Appendix C Non-adiabatic transition amplitude

The equations for $C_{\pm}(t)$ in Sec. IV can be rescaled to the form of Weber differential equation,

[TABLE]

The general solution to this equation is a linear combination of two parabolic cylinder functions Whittaker and Watson (1990),

[TABLE]

Coefficients $A_{\pm},B_{\pm}$ can be found from the initial values of $C_{\pm}(0)$ with account taken of the relation $i\hbar\dot{C}_{\pm}(0)=\Delta C_{\mp}(0)$.

Using the asymptotic expansion $D_{q}(z)\approx\exp(-z^{2}/4)z^{q}$ for $|z|\rightarrow\infty,|\rm{arg}\,z|<\frac{3}{4}\pi$, we find to the first order in $1/t$

[TABLE]

where $\theta(t)$ is given by Eq. (14).

For $|\Delta|\ll\nu(t)$ we have $C_{\uparrow}\approx C_{+}+(\Delta/2\nu)C_{-}$ and $C_{\downarrow}\approx C_{-}-(\Delta/2\nu)C_{+}$. One can then immediately find the coefficients $\alpha_{\uparrow,\downarrow},\beta_{\uparrow,\downarrow}$ in Eq. (13). In particular, $\alpha_{\uparrow}=B_{+}\alpha_{+},\alpha_{\downarrow}=B_{-}\alpha_{-}.$

Of primary interest to us is the limiting value $C_{\uparrow,\downarrow}(\infty)\propto\alpha_{\uparrow,\downarrow}$. For the considered initial condition $C_{+}(0)=C_{-}(0)=1/\sqrt{2}$, we find that

[TABLE]

where the upper sign refers to $\alpha_{\uparrow}$ and the lower sign refers to $\alpha_{\downarrow}$; $\Gamma(x)$ is the gamma function.

The expressions for $\alpha_{\uparrow,\downarrow}$ in the adiabatic limit $p\rightarrow\infty$ can be obtained from Eqs. (29) using the asymptotic form of the gamma function $\Gamma(z)$ for $|z|\to\infty$, cf. Abramowitz and Stegun (1972). They were used in Eq. (15).