Classical and quantum dissipative dynamics in the Josephson junction: an Arnold problem, bifurcation and capture into resonance

# Classical and quantum dissipative dynamics in the Josephson junction: an Arnold problem, bifurcation and capture into resonance

Dmitrii Pashin and Arkady M. Satanin Department of Theoretical Physics, Lobachevsky State University of Nizhny Novgorod, Nizhny Novgorod 603091, Russia    Chang Sub Kim Department of Physics, Chonnam National University, Gwangju 61186, Korea
July 14, 2019
###### Abstract

We study theoretically nonlinear dynamics in the Josephson junction, which maps onto an oscillating motion of an artificial atom under external forces and damping. With an appropriate driving condition, the Josephson particle undergoes the intriguing bistable dynamics near a saddle point in the energy landscape. The bifurcation mechanism provides a crucial component in superconducting quantum circuits with relevance to nondemolition measurements such as a high-fidelity readout of qubit states. We address the question of what the probability is of capture into either basin of attraction and answer to it both in classical and quantum dynamics. Consequently, we derive Arnold’s probability formula and numerically analyze its implementation of the controlled dynamical switching between two steady states under the various nonequilibrium conditions.

###### pacs:
05.10.-a, 85.25.Cp, 03.65.Yz, 02.50.-r

## I Introduction

Today superconducting quantum circuits with embedded Josephson junctions, in hundreds of nanometer size, can be artificially tailored, which behave like atoms characterized by well-defined discrete energy levels of confined Cooper pairs Makhlin2001 (); You2005 (); Clarke2008 (). Such a superconducting circuit provides an arena where experiments can prepare, manipulate, and readout the quantum states at macroscopic scales. When two levels, not necessarily energy eigenstates, can be isolated by an appropriate setup, the system forms quantum bits (qubits) which are the primitive building blocks of quantum computers. Accordingly, superconducting circuits have been one of the intensive research subjects from the perspectives of both fundamental interest and realization of quantum computing Vion2002 (); Greenberg2002 (); Cottet2002 (); Collin2004 (); Ithier2005 (); Siddiqi2005 (); Siddiqi2006 (); Picot2008 (); Mallet2009 (); Nakano2009 (); Meister2015 ().

In this work, we theoretically consider a current-biased Josephson junction which serves as a key component in various superconducting circuits. Our attention here is given to a control of the dissipative dynamics in the artificial system, not to a realization of a qubit or its implication in realistic circuits. The Josephson junction we consider is a nonlinear oscillator composed of an inductor, made of the gauge-invariant phase variable, connecting with an intrinsic capacitor in parallel. Also, we drive the system by applying an alternating current with a detuning from the natural frequency. This electronic component acts as a two-terminal switching device and is a promising bifurcation amplifier with substantial gain but negligible dissipation Siddiqi2004 (); Vijay:2009 (); Kakuyangi2013 (). We place the system in thermal contact with a heat reservoir so that the coupling provides an intrinsic dissipation mechanism in our model. The ensuing dynamics near a bifurcation point is very sensitive to perturbation; for instance, a weak interaction with a qubit will subserve the perturbation. By assuming the feasibility of controlling initial states of the nonlinear system, we monitor the relaxation dynamics numerically in great detail. Consequently, we are allowed to describe quantitatively how a bifurcation develops in the system, which is an essential feature for the switching function, both in the classical and quantum regimes.

The primary goal of this paper is to answer to what the probability would be of the Josephson atom being captured onto either stable steady states both in the classical and quantum regimes. This problem was studied earlier in a charged particle motion under electromagnetic fields Lifshitz1962 () and was rigorously framed into mathematical mechanics by V. I. Arnold Arnold1963 (), and is known as an Arnold problem in nonlinear dynamics Neishtadt1991 (); Haberman1995 (); Arnold:Book2006 (). Also, it was revisited recently in a simple quantum dynamics by the present authors Pashin2017 (). The switching between two coexisting stable states in the classical Josephson oscillator has drawn previously much attention among researchers. However, to our knowledge, the consideration of the probability formulation in quantum dynamics problems is rare.

We adopt the conventional one-particle Hamiltonian to describe the Josephson junction occupied by a macroscopic number of Cooper pairs, where the charging energy and power input terms play the role of kinetic and potential energies, respectively. We approximate the potential energy up through quartic terms, thus casting our problem to a Duffing oscillator in nonlinear dynamics Nayfeh1979 (). The external current adds a time-dependent term to the potential energy, which renders the effective Hamiltonian to be nonautonomous. In this scenario, we view the Josephson junction, a nonlinear oscillator or switching device, as a system of noninteracting single particles, which we name the ‘Josephson particles’. We address the Arnold problem first in the classical regime and then continue formulating it in the quantum regime. The classical description prevails when the thermal energy exceeds the energy scale defined by the characteristic frequency of the Josephson oscillator, providing an insightful physical picture of the nonlinear dynamics in the complex energy landscape. At low temperatures, however, the thermal energy is negligible compared to the oscillator energy-level spacing. Accordingly, we employ the full quantum description to account for the bifurcation probability.

In the classical regime, it is well known that a nonlinear oscillator can be switched between two stable states by tuning the external force beyond some threshold magnitude Nayfeh1979 (); Bogoliubov1962 (). In the Josephson oscillator, we find it useful to parameterize the superconducting phase variable, the oscillator amplitude, in terms of ‘in-phase’ and ‘quadrature-phase’ components with respect to the driving oscillation Siddiqi2004 (). Consequently, we convert the dynamical variables into a canonically conjugate pair and solve Hamilton’s equation of motion for the effective variables. Furthermore, we smooth out the fast oscillatory time-dependence in the equations of motion from the driving current to make the dynamics approximately autonomous. Then, the phase-space trajectories are examined to learn the build-up of the bifurcation dynamics on the coarse-grained time scale.

In the quantum regime, we introduce a unitary transformation to consider the dynamics in a rotating-wave frame. Consequently, we render the Josephson particle to obey the quasi-stationary Schrödinger equation on a time scale longer than the oscillation period of the external driving. Accordingly, we solve the single-particle eigenvalue problem numerically and obtain the energy-eigenvalue spectrum and quasi-stationary eigenfunctions. The obtained single-particle states form the discrete energy levels of our macroscopic artificial atom. In practice, the superconducting circuit is embedded in a heat reservoir such as a cryogenic refrigerator. Accordingly, an initially prepared system in a pure state with specific energy by an external pumping must relax to a steady state, an incoherent superposition of mixed states. To investigate the relaxation mechanism, we model the system as embedded in a thermal reservoir. Subsequently, we set up a statistical-mechanical Liouville equation for the combined system and reservoir Blum1981 (), and derive the master equation for the probability occupancy of the Josephson levels. Then, we explore the bifurcation dynamics and show the quantum version of the Arnold formula matches with the classical-mechanical representation.

We note here that Duffing oscillators with quartic nonlinearity, driven by a time-dependent periodic force, have been studied widely among researchers in both classical and quantum dynamics. For instance, there has appeared a series of applications of the Duffing equation Dykman1979 (); Dykman1988 (); Dykman2007 (); Andre2012 (); where it is reported that fluctuations provide a mechanism to overcome the dynamical barrier, consequently to induce switching between the stable solutions, accompanying dissipation of the oscillator states. We also draw our attention to other works which are devoted to tunneling problem between two stable attractors Sazonov1976 (); Dmitriev1986 (); Wielinga1993 (); Peano2004 (); Marthaler2007 (). However, unlike our contribution to the quantum formulation, however, most of the previous works were carried out in the semiclassical approximation. In addition to the above-cited works, we also refer the review articles in which these and related subjects are described Dykman1984 (); Grifoni1998 (); Dykman2012 ().

This paper is organized as follows. First, in Sec. II, we address the Arnold question in a simple model to set the stage for the Josephson system. In Sec. III, we study the classical dynamics of the Josephson particle in phase space and establish the Arnold probability. Then, we treat the problem quantum mechanically in Sec. IV: In Sec. IV.1, we first solve the quasi-stationary eigenvalue problem of the Josephson particle. In the subsequent Sec. IV.2, by employing the density matrix, we formulate the relaxation dynamics of the pumped system in a heat bath and derive the quantum analog of the Arnold formula. Finally, in Sec. V, we provide a summary and the conclusion.

## Ii The Arnold problem

Here we recapitulate the original Arnold problem which considers a classical particle with mass in a static double-well potential Arnold1963 (). The motion of the particle is governed by the equations of motion in phase space,

 ˙x=pmand˙p=−∂V(x)∂x−γp;

where represents the potential energy and is a damping coefficient.

The separatrix is a phase portrait of the particle with energy matching with the top value () of the unstable barrier in the double well, which is an ideal trajectory in phase space in the limit of vanishing dissipation, . In Fig. 1 such a separatrix is depicted where the crossing point corresponds to the unstable equilibrium point in the potential energy. The separatrix defines the boundary in phase space separating two distinct modes in particle motion under damping: 1) All states outside the separatrix possess energies bigger than and will tend to cross a point, with time, on one of the two branches of the separatrix curve; 2) All states inside the separatrix possess energies less than and will relax to either stable equilibrium point, depending on initially which lobe they belong to.

In the small damping limit, by applying the work-energy theorem, one can calculate the energy changed over a full cycle of the trajectory, approximately close to the separatrix, to linear order in ,

 ΔH=−γ∫p˙xdt→−γ∮pdx=−γ(SL+SR).

If the particle initiates its motion from a phase point outside the separatrix, it will cross the separatrix in time by dissipation to stochastically enter one of the lobes. It is well known that the probability of capture into each basin of attraction is proportional to the change in the Hamiltonian during one cycle of the corresponding homoclinic trajectory, see, for instance, Haberman1995 (). Accordingly, the probability of the particle to fall onto either the left (L) or right (R) equilibrium states would be proportional to the bounded area of each lobe, or . We call the resulting probability the Arnold formula that is written as

 Pα=SαSL+SR,α=L, R; (1)

which is evidently independent of the damping strength.

## Iii A single Josephson junction: Classical dynamics

The Josephson junction is macroscopically described by the persistent current of electron pairs,

 I=Icsinφ,

where is the relative phase of the condensed states between two superconducting sides with being the critical current, and the voltage across the junction,

 V=Φ0˙φ,

where is the reduced magnetic-flux quantum. Accordingly, the classical Hamiltonian for the Josephson junction in the presence of the external current may be written as

 H=12C(Φ0˙φ)2−EJcosφ−Φ0Iexφ. (2)

The first term on the right hand side (RHS) of Eq. (2) is the charging energy, in the junction; where is the intrinsic capacitance of the Josephson junction. The second term is a potential energy associated with the power input, , where , the Josephson energy. The last term describes time-dependent external control from the driving current ,

 Iex(t)=I0cos(ωt).

This Hamiltonian can be viewed as a function of two canonically conjugate variables, and , i.e. . Then, the standard Hamiltonian formulation yields

 ˙φ = ∂H∂pφ=1CΦ20pφ, ˙pφ = −∂H∂φ=−EJsinφ+Φ0Iex;

which generate the Newtonian equation of motion for the phase variable:

 CΦ20¨φ+EJsinφ=Φ0Iex.

In addition, to account for dissipation we insert the Ohmic current, with being the intrinsic resistance, in the preceding equation to generalize it as

 ¨φ+1RC˙φ+EJCΦ20sinφ=1Φ0CIex. (3)

Then, we define the various coefficients as

 γ≡1RC,ω0≡√EJCΦ20,ωp≡√I0Φ0C (4)

and further introduce the various dimensionless variables as

 τ≡ω0t,¯γ≡γω0,¯I0≡ω2pω20,and¯ω≡ωω0.

With these arrangement the equation of motion, Eq. (3) is cast into the dimensionless form,

 d2φdτ2+¯γdφdτ+sinφ=¯I0cos(¯ωτ)

which represents a nonlinear, damped harmonic oscillator with pumping. For a complete analysis we shall work in the regime where the relative phase is small across the junction so that we may approximate up to the cubic term. Then, the equation of motion to work with becomes a Duffing equation,

 ¨φ+¯γ˙φ+φ−16φ3=¯I0cos(¯ωτ), (5)

where is understood to be time-derivative with respect to dimensionless time . Note Eq. (5) takes a generic from which includes the conservative, dissipative, and external forces.

To study the Arnold problem we find it useful to parameterize the phase variable in terms of in-phase and quadrature-phase components with respect to the driving oscillation Siddiqi2004 (): We convert the dynamical variables as

 φ=qcos(¯ωτ)+psin(¯ωτ), (6) ¯ω−1˙φ=−qsin(¯ωτ)+pcos(¯ωτ) (7)

of which second equation constrains that

 ˙qcos(¯ωτ)+˙psin(¯ωτ)≡0. (8)

In order to furnish the equations of motion for the new canonical variables and we take time-derivative of Eq. (7) with holding the constraint, Eq. (8). Consequently, we get

 ¨φ=¯ω(−˙qsin(¯ωτ)+˙pcos(¯ωτ))−¯ω2φ,

that can be substituted into Eq. (5) to brings about

 −˙qsin(¯ωτ)+˙pcos(¯ωτ)=1¯ω(¯I0cos(¯ωτ)+(¯ω2−1)φ+16φ3−¯γ˙φ).

Next, after substituting Eqs. (6) and (7) for and on the RHS of the preceding equation, respectively, one can solve it, which is linearly coupled with Eq. (8), for and to get

 ¯ω˙p = + 16{q3cos4(¯ωτ)+p3cos(¯ωτ)sin3(¯ωτ)+3q2pcos3(¯ωτ)sin(¯ωτ) + 3qp2cos2(¯ωτ)sin2(¯ωτ)}−¯γ¯ω{−qsin(¯ωτ)cos(¯ωτ)+pcos2(¯ωτ)}, ¯ω˙q = −¯I0cos(¯ωτ)sin(¯ωτ)−(¯ω2−1){qcos(¯ωτ)sin(¯ωτ)+psin2(¯ωτ)} − 16{q3cos3(¯ωτ)sin(¯ωτ)+p3sin4(¯ωτ)+3q2pcos2(¯ωτ)sin2(¯ωτ) + 3qp2cos(¯ωτ)sin3(¯ωτ)}+¯γ¯ω{−qsin2(¯ωτ)+pcos(¯ωτ)sin(¯ωτ)}.

All the time-dependent terms on the RHS of the preceding equation oscillate faster than or equal to . We then take time-average of the last equations of motion over the half-period, , of the external force to smooth out the fast oscillatory time-dependence. The outcome of this slowly varying amplitude approximation (SVAA) is given as

 ˙p = −αq+βq(q2+p2)+f−¯γ2p, (9) ˙q = αp−βp(q2+p2)−¯γ2q, (10)

where the coefficients have been defined as

 α≡1−¯ω22¯ω,β≡116¯ω,andf≡¯I02¯ω.

Note that there appears a non-conventional damping term, , in Eq. (10)

Eqs. (9) and (10) constitute the coarse-grained equations of motion of the classical variables describing the Josephson dynamics. For further analysis we find it useful to construct the effective Hamiltonian function that generates the conservative parts of Eqs. (9) and (10). By inspection, we have obtained the effective Hamiltonian as

 ~H(q,p)=α2(q2+p2)−β4(q2+p2)2−fq, (11)

where is normalized with respect to the reference energy, , which gives the same energy scale with the original Hamiltonian, Eq. (2).

Unlike the usual Hamiltonian this function cannot be separated into the kinetic energy plus potential energy. However, one may still define the energy of the system as an instant value of the Hamiltonian function, given and ,

 ~E≡~H(q,p).

As a result of the SVAA, the original oscillatory time-dependence has disappeared in Eq. (11) that the Hamiltonian becomes autonomous approximately. In below, for later purposes, we shall present the numerical data for energies in rescaled units of . Also, reflecting the cubic-order expansion of the potential force , we consider the energy only in the finite range of .

In Fig. 2 we illustrate typical energy surfaces described by Eq. (11), where the numerical values for and have been estimated using the physical parameters in an experiment Siddiqi2004 (),

 α=1.30×10−1andβ=7.11×10−2.

We shall use these values throughout the calculation below. The energy landscape appears intriguing as shown in Fig. 2(a). The steady-state condition brings about three fixed points all of which happen to be at : Two of them are stable and the other is a saddle point. Between two stable extrema, one appears on the surface top, , with energy and the other at the bottom of the well, , with energy . We have performed linear stability analysis to confirm that both points are centers. The saddle point appears at with energy . These features are better viewed in Fig. 2(b) where we draw the energy contour on -plane. For reference, we give the numerical values of , , and determined from the adapted parameters for references,

 ~EL=33.1, ~Esd=2.60, and ~ER=−1.79.

In what follows we shall call the phase point at the -point and the one at the -point; where and stand for ‘left’ and ‘right’ to on axis, respectively.

Then, we view Eqs. (9) and (10) for the SVAs as a generalized Hamiltonian dynamics generated by the effective Hamiltonian, Eq. (11), with the additional dissipative terms,

 ˙p = −∂~H∂q−¯γ2p, (12) ˙q = ∂~H∂p−¯γ2q. (13)

Note that the dynamical description by Eqs. (12) and (13) contains not only a generalized force but also a generalized velocity, exhibiting an extended dynamics in the generic form Kim:2015 ().

Like the original Arnold problem, the trajectories generated by the Josephson dynamics may be classified by a separatrix which is the unperturbed phase-trajectory to that the particle with the saddle point energy tends ideally in the limit of vanishing dissipation, . In Fig. 3 we depict a separatrix resulting from the parameters and used in Fig. 2 for the same , which comprises two homoclinic orbits, a homoclinic orbit being the phase trajectory joining a saddle point to itself. According to the separatrix the mechanical states of the particle are categorized into three classes: 1) The phase points inside the inner lobe is limited in the energy window and all of which generate trajectories swirling onto the -point in time regardless of the damping strength; 2) The states confined between the outer homoclinic and the inner homoclinic (shaded area) possess energies, and relax to either the - or -point, depending on the size of damping, see figure 6.; 3) The states outside the larger homoclinic possess energies less than the saddle-point energy, . When an initial state is chosen among these phase points, the trajectory again falls onto again either - or -point, depending on the damping magnitude.

We have performed the numerical experiments to monitor the ensuing trajectories by solving Eqs. (12) and (13). To run the dynamics, an arbitrary initial state has been chosen in each of three distinguished regions in phase space, divided by the separatrix, for two different damping strengths. The outcome is as follows: 1) When an initial state is chosen inside the inner lobe, the trajectory falls onto the -focus regardless of the damping magnitude. 2) When initial state is located inside the shaded area in Fig. 3, the trajectory falls onto either the -focus or -focus, depending on the damping strength. 3) When we choose an initial state outside the large separatrix, again like the case 2), there is a finite probability for the resulting trajectory to fall onto either focus for a finite damping. In Fig. 4 we depict the trajectories for the case 3).

Having learned that the dynamics, generated from initially chosen states placed outside the seraratrix in phase space, manifests a bistable behavior in the vanishing damping limit, it is of interest to determine what the respective probability of the resulting trajectories to relax stochastically onto either the - or the - point is. It constitutes the Arnold problem. To study the Arnold question quantitatively, we consider the change in the energy along a segment of the separatrix for motion with a small damping (), that can be calculated as

 Δ~H = ∫(∂~H∂p˙p+∂~H∂q˙q)dτ (14) = ∫{(˙q+¯γ2q)˙p+(−˙p−¯γ2p)˙q}dτ = ¯γ2∫(q˙p−p˙q)dτ.

Although the phase-trajectories are not exactly periodic for the considered dissipative dynamics, we may assume that in the limit of vanishing damping the particle returns to an initial phase point after a cycle within an error of . Accordingly, the energy change of the particle along the outer closed loop, an homoclinic orbit, on the separatrix is evaluated as

 Δ~Hout = ¯γ2∮out(qdp−pdq)=¯γ2{(SL+SR)−(−(SL+SR))} (15) = ¯γ(SL+SR),

where is the phase-space area enclosed in the inner lobe, and is the area enclosed in the outer lobe subtracted by , thus the shaded area. Similarly, the energy change along the inner closed loop, again a homoclinic orbit, is calculated to be

 Δ~Hin=¯γ2(−SR−SR)=−¯γSR. (16)

Thus, the energy increases when the particle evolves around the outer loop in the counterclockwise manner; while it loses energy along the clockwise inner loop on the separatrix.

Now we choose an initial state outside the external contour of the separatrix, see Fig. 3, whose corresponding energy is below , and simulate its motion for a small damping . Fig. 5 illustrates the trajectory at early times; where it is seen that the trajectory swirls about the outer region of the separatrix counterclockwise. The particle keeps evolving while gradually reducing its radius of curvature with time. In this process the energy of the particle progressively augments and must reach the value at an instant. At this moment we say that the particle ‘enters’ onto the outer separatrix, and have indicated such an entering point by the mark on the separatrix loop in Fig. 3. We watch, upon entering, that the ensuing trajectory follows the outer separatrix for some time and turns its direction clockwise to swirl around the inner separatrix. As time elapses further, the trajectory falls either onto the - or onto the -basin of attraction, depending on the sign of the successive energy change after the entering. In Fig. 6 we present this finding from two selected initial states. We observe in Fig. 6(a) that the trajectory started out from the initial point has fallen into the -basin of attraction and continues relaxing to the -point. In contrast, in Fig. 6(b) we see that the initial state evolves into the -basin of attraction.

Next we analyze dynamics upon ‘entering’ the outer separatrix in great detail. The Arnold problem addresses the question in the vanishing damping limit, . In principle, the limiting trajectory would tend to the separatrix on an extremely long time scale, which cannot be achieved numerically. For the sake of argument, we take a portion of the outer homoclinic orbit as an approximate trajectory of motion just after entering takes places. To be more specific, let us denote as the energy at point in Fig. 3 whose value equals . After entering, the particle continually evolves along the aforementioned separatrix segment with gaining energy by the continuing counterclockwise motion around the outer loop and with losing energy along the inner branch clockwise. As the particle approaches to the saddle point, marked by in Fig. 3, its accumulated energy becomes

 ~E=~EA+δ~Eout+δ~Ein,

where represents the energy gain along the elapsed segment of the outer loop, and represents the energy loss along the inner loop, which equals from Eq. (16). If , the arrival energy of the particle at exceeds the saddle point value . Accordingly, the particle must evolve into the -basin of attraction to eventually reach the L-point. This situation is depicted in Fig. 6(a). On the other hand, if , the net change in energy becomes negative so that the approaching energy to the saddle point becomes less than the actual saddle point energy . Then, the particle trajectory must fall into the -basin of attraction as depicted in Fig. 6(b). Thus, entering states onto the separatrix are categorized into two classes by the very entering position at which the net-energy change balances, , so : Those entering states onto the lower segment () will be attracted to the -basin with an accumulated loss of energy . And, those states entering onto the upper segment () will gain more energy than and consequently accumulate a net energy to fall into the -basin of attraction. Recall that the energy change along the complete loop, plus , adds up to give the full area enclosed by the outer homoclinic orbit, [Eq. (15)]. Accordingly, because , the energy change must be . The desired bifurcation probability is known to be proportional to the energy changes Haberman1995 (). Thus, the chances of an initial state, which enters onto either loop segment, subsequently relaxing to either steady state, are proportional to the area of the corresponding basin of attraction, i.e. or :

 Pα=¯γSα¯γ(SL+SR)→SαSL+SR,α=L, R;

which takes the same form as Eq. (1).

In short, our analysis has revealed that the Arnold formula also holds in the classical Josephson dynamics: The probability of the particle, with initial energy below the bifurcation value whose classical phase point resides outside the larger homoclinic orbit in the separatrix, being captured onto one of the two stable steady-state points is proportional to the phase-space area swept by the corresponding homoclinic loop.

## Iv Quantum dynamics in the Josephson junction

Now we address the Arnold problem in the quantum Josephson particle and will derive a formula that describes the bistable dynamics in quantum domain.

### iv.1 Single electron-pair spectrum

To consider the Arnold problem quantum mechanically we first translate the dynamical variables in the classical Hamiltonian, Eq. (2) into the Hermitian operators,

 φ→^φandN=12eCΦ0˙φ→^N,

where is the difference in the number of the Cooper pairs, being the net charge stored, across the Joshepson junction. The operators and are canonically conjugate with each other to satisfy the commutator,

 [φ,^N]=i.

Then, the quantum behavior of the Josephoson junction is described by the Hamiltonian operator,

 ^H=Ec^N2−EJcos^φ−Φ0Iex^φ, (17)

where is the charging energy per electron pair. The last term on the RHS of Eq. (17) is associated with the driving current, see Sec. III.

To proceed we suppose that the relative phase is small as we have done in the classical calculation to expand the potential operator as

 cos^φ≈1−^φ2/2!+^φ4/4!.

Also we introduce the creation and annihilation operators to linearly transform as

 ^N = i(EJ8EC)1/4(^a†−^a), (18) ^φ = (EC2EJ)1/4(^a†+^a). (19)

One can confirm that the new operators satisfy the commutation relation for harmonic oscillators,

 [^a,^a†]=1.

Subsequently, by plugging Eqs. (18) and (19) into Eq. (17) we convert the Hamiltonian, Eq. (17) into

 ^H = ℏω0(^a†^a+12)−ε(^a+^a†)4−f0cos(ωt)(^a+^a†)−EJ, (20)

where with being the natural frequency. And, other definitions and are related with the parameters and in the classical Hamiltonian, Eq. (11), as

 ε≡148EC=13ℏω√EC2EJβ, f0≡I0Φ0(EC2EJ)1/4=ℏω(2EJEC)1/4f.

The obtained Hamiltonian is time-dependent with periodicity of the driving current , .

Here we find it useful to define the unitary transformation,

 ∣Ψ⟩→∣ΨRWA⟩=^U∣Ψ⟩, (21)

where and are the kets governed by the original and the transformed , respectively, and the unitary operator is defined to be

 ^U=eiω^a†^at (22)

which is a quantum version of the classical parametrization, Eqs. (6) and (7). Subsequently, we show that the two Hamiltonians are related to each other in accordance with

 ^HRWA=^U^H^U†−ℏω^a†^a, (23)

and the quantum dynamics is described by the time-dependent Schrödinger equation in the rotating-wave frame,

 iℏ∂∂t∣ΨRWA⟩=^HRWA∣ΨRWA⟩. (24)

The explicit form of the transformed Hamiltonian in the RWA is given as

 ^HRWA = (ℏ(ω0−ω)−6ε)(^a†^a+12)−6ε(^a†^a)2 − ε(e−4iωt^a4+e−2iωt^a3^a†+e−2iωt^a2^a†^a+e−2iωt^a^a†^a2+e2iωt^a^a†3 + e−2iωt^a†^a3+e2iωt^a†^a^a†2+e2iωt^a†2^a^a†+e2iωt^a†3^a+e4iωt^a†4) − f0cos(ωt)(eiωt^a†+e−iωt^a)−EJ,

in obtaining which we have used the operator identities,

 eiω^a†^at^a†=^a†eiω^a†^ateiωtandeiω^a†^at^a=^aeiω^a†^ate−iωt,

which give rises to, for instance, . Next, we take an average of the preceding over the half-period of the driving current as we performed it in the classical SVAA. Consequently, we obtain the coarse-grained Hamiltonian for quantum analysis of the Arnold problem, up to a constant, as

 ^HRWA→^Hrwa≡ℏ¯ω0^a†^a−6ε(^a†^a)2−f02(^a†+^a), (25)

where we have set . Note all the oscillatory time-dependences have been smoothed out in Eq. (25) and the Hamiltonian becomes approximately time-independent, more precisely autonomous.

Thus, the system is in the stationary state in the slowly varying RWA and is described, in general, as

 ∣Ψrwa(t)⟩=∑jaje−iℏ~Ejt∣ϕj⟩, (26)

where and are solutions to the energy-eigenvalue equation, We note that the energy, defined here on the coarse-grained time scale, is different from the ‘quasienergy’ that holds when a periodically driven Hamiltonian is considered strictly Shirley1965 (); Zeldovich1966 (); Ritus1966 ().

 ^Hrwa∣ϕj⟩=~Ej∣ϕj⟩. (27)

We have solved the time-independent Schrödnger equation, equatioin (27) in the number representation (Fock basis),

 ∣ϕj⟩=∑c(j)n∣n⟩, (28)

where are the eigenkets of the number operator and is the expansion coefficient. In doing so, we have used the numerical values, adopted from Siddiqi2004 (),

 ω=0.878ω0,ε=3.28⋅10−5ℏω0,andf0=0.89ℏω0;

which are equivalent to the parameters used in Sec. III. The natural frequency for tunnel junction is estimated as which gives a temperature scale . Also, the charging energy per Cooper pair and Josephson energy are estimated to be erg and erg.

The result is given in Fig. 7 where we present the energy spectrum as a function of the average phase, , of the Josephson junction at quantum level , which corresponds to the slowly varying amplitude in classical dynamics. In classical dynamics, quantum average must be replaced with time average. The time-average of the classical representation, Eq. (6), for over the symmetric half-period, , brings about , the average amplitude. It is seen that, as we follow the energy levels from below, the corresponding stationary position of the particle moves from the origin toward the bottom of the -well and continues to the positive direction but not quite reaches the saddle point. Then, it turns its direction and keeps climbing up to reach the -point at the highest level (). The states marked by the open circles inside the well, , cover the RWA levels . And, intriguingly, they are nearly doubly-degenerated, see the inset. The closest levels to inside the -well are and belonging to the inner lobe and the outside the large separatrix in Fig. 3, respectively. We have numerically checked that increasing the number of basis vectors simply produces more negative energies without affecting the upper branch in the spectrum. Note the energy spectrum of the quasi-Hamiltonian is bounded from above, not from below.

Having solved the problem in the RWA, we now consider the dynamics governed by the original Hamiltonian by performing the inverse transformation of Eq. (21), which is done approximately in the SVAA. The temporal development of the state is then given by

 ∣Ψ(t)⟩=e−iω^a†^ate−iℏ^Hrwat∣Ψ(0)⟩, (29)

where is an arbitrary initial state that may be expanded in the rotating-wave basis ,

 Ψ(0)=∑aj∣ϕj⟩.

For later purposes we shall recast the above Eq. (29) into

 ∣Ψ(t)⟩=∑jaje−iℏ~Ejt∣ψj(t)⟩, (30)

where we have used the expansion, Eq. (28) to define

 ∣ψj(t)⟩≡e−iω^a†^at∣ϕj⟩=∑nc(j)ne−inωt∣n⟩, (31)

which satisfies the periodicity of ,

 ∣ψj(t+2π/ω)⟩=∣ψj(t)⟩.

In particular, if the system is prepared initially in a RWA eigenstate, say , it evolves simply as

 ∣Ψ(t)⟩→e−iℏ~Ejt∣ψj(t)⟩≡∣ΨEj(t)⟩, (32)

which we interpret as a quasi-stationary state. Although the original Hamiltonian is non-autonomous, the energy of the system in the quasi-stationary state may be still defined as the expectation value, . which can be calculated using Eq. (23):

 ⟨ΨEj∣(^U†^Hrwa^U+ℏω^a†^a)∣ΨEj⟩ = ⟨Ψrwa∣^Hrwa∣Ψrwa⟩j+ℏω⟨ΨEj∣^a†^a∣ΨEj⟩ = ~Ej+ℏω∑nn|c(j)n|2.

Thus, we get the quasi-stationary energy associated with the state as

 Ej=~Ej+ℏω∑nn|c(j)n|2. (33)

Note we use the notation to denote the energy expectation value of , distinguished from for RWA energies. The obtained approximate energy eigenvalues of are time-independent and the states of the system evolve in the quasi-stationary manner on the coarse-grained time scale longer than .

We present the energy spectrum of , Eq. (33) in Fig. 8 which is bounded from below in contrast to that of . Because of the contribution from the second term on the RHS in Eq. (33), the quasi-eigenvalue corresponding to each level has been changed in comparison with Fig. 7. The splitting structure in the energy band arises from the nearly degenerated levels which occupy the RWA-energy well alternatively in Fig. 7: The second term on the RHS in Eq. (33) produces distinguished contributions to two nearly degenerated states. Consequently the quasi-degeneracy in the spectrum is lifted to bring about two distinctive branches in the original spectrum. The highest level in Fig. 7 is not the highest level anymore in Fig. 8. And, the lowest level on the red-colored branch in the energy spectrum gives the ground state of the physical Hamiltonian . Its paired partner near the bottom of the -well in the classical energy contour in Fig. 7 appears on the upper (black-colored) branch with a larger energy. Also the sparse levels near the mark B are the ones that approach close to the saddle point in Fig. 7. The two pronounced levels at and near the mark B correspond to the highest degenerated states in the RWA well in Fig. 7. For reference, we have checked that the level spacing between adjacent levels on the red branch is in order of . Also, we have drawn the energy levels only up to by reflecting the quadruple-order expansion of the potential operator.

Next, we examine time-evolution of the quasi-stationary state , see Eq. (32). The corresponding wavefunction is given explicitly in -representation as

 ΨEj(q,t)=e−iℏ~Ejtψj(q,t), (34)

where

 ψj(q,t)=∑nc(j)ne−iωntHn(q),

where are the eigenfunctions of the harmonic oscillator.

In Fig. 9 we draw the squared amplitude of the wavefunctions, , for a few levels. The pronounced characteristic of the quasi-stationary states is that their wave functions oscillate resonantly with the external frequency, , while the shape of the wave packets being modulated over the period. The two particularly well-localized wavepackets in Fig. 9 are from the levels and in Fig. 8, which correspond to the top of -hill and the bottom of -well in the classical RWA energy contour, respectively, see Fig. 7. The wavepacket marked by II is the one corresponding to , marked by B, which is mostly close to the saddle point in the classical energy contour. The wavepackets marked by I and III are from the two split states, and in Fig. 8, which, in turn, are the degenerated levels in Fig. 7. In fact, all the paired levels, one on black and its partner on red branches in Fig. 8, which are split from the corresponding nearly-degenerated levels inside the RWA energy well, manifest the same feature, one level is delocalized and the other localized in the pair. In the classical picture, the localized levels, with nearly same RWA energies in the window , belong to the inner lobe and the delocalized levels belong to the region outside the large lobe in the separatrix. Exceptionally we found that the paired, split levels of (488,489), (490,491), and (492, 493), which are inside the well and very close to , exhibit all extended states.

### iv.2 Density matrix formulation of the Josephson dynamics in a boson bath

So far, we have considered the single-particle spectrum of the Josephson particle and quasi-stationary dynamics in a pure state without dissipation. We now suppose that the Josephson junction described by , Eq. (20) is placed in a heat reservoir at an absolute temperature . For the composite system of the Josephson junction and the reservoir, we write the total Hamiltonian as

 ^Htot=^H+^HR+^V, (35)

where is the Hamiltonian of the reservoir composed of a number of bosonic modes,

 ^HR=∑iℏΩi^b+i^bi, (36)

and is the interaction between the system and the reservoir, simply modeled as

 ^V=(^a++^a)∑iκi(^b+i+^bi), (37)

where is the coupling constant between the Josephson particle with th bath mode. So, our model introduces dissipation via the system-reservoir interaction in the independent particle picture. Because of the interaction the Josephson particle must be, in general, in a mixed state. To describe the system in the density-matrix formulation, we first set up the quantum Liouville equation for the total system.

 iℏ∂^ρtot(t)∂t=[^Htot,^ρtot(t)], (38)

where denotes the density operator of the total system. The system alone is described by which is the reduced density operator from ,

 ^ρ(t)=TrR^ρtot(t),

where indicates the trace over the reservoir states.

Our goal here is to construct a closed equation for which describes the desired dissipative quantum dynamics of the system. To this end, we find it useful to define the transformation via

 ~ρtot(