Quantum properties of a strongly driven Josephson junction

TL;DR

This paper investigates the quantum behavior of a strongly driven Josephson junction, revealing complex dynamics beyond classical models and highlighting the influence of quantum fluctuations on its nonlinear properties.

Contribution

It provides a detailed quantum analysis of a driven Josephson junction, extending beyond the Duffing approximation and exploring classical-quantum correspondence in nonlinear regimes.

Findings

Identification of parameter regions where Duffing approximation is valid

Analysis of quantum fluctuations affecting junction parameters

Study of quantum effects on dynamics near classical fixed points

Abstract

A Josephson junction embedded in a dissipative circuit can be externally driven to induce nonlinear dynamics of its phase. Classically, under sufficiently strong driving and weak damping, dynamic multi-stability emerges associated with dynamical bifurcations so that the often used modeling as a Duffing oscillator, which can exhibit bi-stability at the most, is insufficient. The present work analyzes in this regime corresponding quantum properties by mapping the problem onto a highly-nonlinear quasi-energy operator in a rotating frame. This allows us to identify in detail parameter regions where simplifications such as the Duffing approximation are valid, to explore classical-quantum correspondences, and to study how quantum fluctuations impact the effective junction parameters as well as the dynamics around higher amplitude classical fixed points.

Equations

$$H_{\text{JJ}}=\frac{q^2}{2C}-E_J\cos (\phi )-F\cos ({\omega }_{d}t)\phi$$

$$H_{\text{JJ}}^{\text{RF}}=\hslash \delta \omega \hat{n}-\frac{F}{2}\kappa (a+a^{†})-E_J^{\ast }\left[:J_0\left(2\kappa \sqrt{\hat{n}}\right):+{\kappa }^2\hat{n}\right],$$

$$H_{\text{JJ, 2}}^{\text{RF}}=\hslash \delta \omega \hat{n}-E_J^{\ast }\frac{{\kappa }^4}{4}({\hat{n}}^2-\hat{n})-\frac{F}{2}\kappa (a+a^{†}),$$

$$H_{\text{Duff}}^{\text{RF}}=\hslash \delta {\omega }_0\hat{n}-E_{\text{J}}\frac{{\kappa }^4}{4}\left({\hat{n}}^2+\hat{n}\right)-\frac{F}{2}\kappa (a+a^{†}),$$

$$f=\frac{F}{M\omega \delta \omega },ϵ=\frac{E_J^{\ast }}{M\omega \delta \omega }$$

$$\begin{array}{rl}g(\hat{Q},\hat{P})=& \frac{1}{2}\left({\hat{Q}}^2+{\hat{P}}^2\right)-\frac{f}{2}\hat{Q}\\ & -ϵ\left[:J_0\left(\sqrt{{\hat{Q}}^2+{\hat{P}}^2}\right):+\frac{1}{4}({\hat{Q}}^2+{\hat{P}}^2)\right].\end{array}$$

$$\begin{array}{rl}{g}_{\mathrm{Duff}}(\hat{Q},\hat{P})=& \left(\frac{\delta {\omega }_0}{2\delta \omega }-\frac{{ϵ}_0{\kappa }^2}{16}\right)\left({\hat{Q}}^2+{\hat{P}}^2\right)-\frac{f}{2}\hat{Q}\\ & -\frac{{ϵ}_0}{64}{\left({\hat{Q}}^2+{\hat{P}}^2\right)}^2,\end{array}$$

$$\frac{dQ}{dt}=-\frac{\tilde{\gamma }}{2}Q-\frac{\partial g}{\partial P},\frac{dP}{dt}=-\frac{\tilde{\gamma }}{2}P+\frac{\partial g}{\partial Q}$$

$${\partial }_t\rho =-\frac{i}{{\kappa }^2}[g,\rho ]+\frac{\tilde{\gamma }}{2}(2a\rho a^{†}-n\rho -\rho n).$$

$$g_2$$

$$\begin{gathered} \left(\begin{array}{c}\delta \dot{a}\\ \delta {\dot{a}}^{†}\end{array}\right) \\ =\left(\begin{array}{cc}-i[2-ϵ+ϵJ_0(2\kappa |\alpha |)]-\frac{\tilde{\gamma }}{2}& iϵ{\alpha }^2{J}_2(2\kappa |\alpha |)/|\alpha {|}^2\\ -iϵ({\alpha }^{\ast }{)}^2{J}_2(2\kappa |\alpha |)/|\alpha {|}^2& i[2-ϵ+ϵJ_0(2\kappa |\alpha |)]-\frac{\tilde{\gamma }}{2}\end{array}\right)\left(\begin{array}{c}\delta a\\ \delta a^{†}\end{array}\right)+\sqrt{\tilde{\gamma }}\left(\begin{array}{c}\delta a_{\mathrm{in}}\\ \delta a_{\mathrm{in}}^{†}\end{array}\right) \end{gathered}$$

$$\delta n_2$$

$${\lambda }_{\pm }=-\frac{\tilde{\gamma }}{2}\pm {\left\{{ϵ}^2{J}_2(2\kappa |\alpha |{)}^2-[2-ϵ+ϵJ_0(2\kappa |\alpha |){]}^2\right\}}^{1/2}.$$

$$\begin{array}{rl}{H}_{\text{JJ}}=& \frac{p_{\phi }^2}{2M}-M{\omega }^2\cos (\phi )-F\cos ({\omega }_{\text{d}}t)\phi \\ =& \frac{q^2}{2C}-E_J\cos (\phi )-F\cos ({\omega }_{d}t)\phi .\end{array}$$

$$\begin{array}{rl}& e^{i{\omega }_{\text{d}}t{a}^{†}a}{e}^{\pm i\sqrt{\frac{\hslash }{2M\omega }}(a+a^{†})}{e}^{-i{\omega }_{\text{d}}t{a}^{†}a}\\ & =e^{\pm i\sqrt{\frac{\hslash }{2M\omega }}{a}^{†}{e}^{i{\omega }_{\text{d}}t}}{e}^{\pm i\sqrt{+\frac{\hslash }{2M\omega }}a{e}^{-i{\omega }_{\text{d}}t}}{e}^{\frac{\hslash }{2M\omega }\frac{1}{2}[a^{†}{e}^{i{\omega }_{\text{d}}t},a{e}^{-i{\omega }_{\text{d}}t}]}\\ & =\sum _{k,u}\frac{{\left(\pm i\sqrt{\frac{\hslash }{2M\omega }}\right)}^k}{k!}{e}^{i{\omega }_{\text{d}}tk}(a^{†}{)}^k\frac{{\left(\pm i\sqrt{\frac{\hslash }{2M\omega }}\right)}^u}{u!}{e}^{-i{\omega }_{\text{d}}tu}{a}^u{e}^{-\frac{\hslash }{4M\omega }}\\ & \stackrel{\text{RWA}}{\underset{\text{k=u}}{=}}:\sum _k\frac{{\left(-\frac{\hslash }{2M\omega }\right)}^k}{k{!}^2}(a^{†}a{)}^k{e}^{\frac{-\hslash }{4M\omega }}:\\ & =e^{-\frac{{\kappa }^2}{2}}:J_0\left(2\kappa \sqrt{\hat{n}}\right):\end{array}$$

$$\begin{split}{H}^{\text{RF}}_{\text{JJ}}=&\hbar\bigl{(}\frac{\omega}{2}-\omega_{\text{d}}\bigr{)}\hat{n}+\frac{\hbar\omega}{4}+M\omega^{2}-M\omega^{2}e^{-\frac{\hbar}{4M\omega}}:J_{0}\bigl{(}2\sqrt{\frac{\hbar\hat{n}}{2M\omega}}\bigr{)}:\\ &-\frac{F}{2}\sqrt{\frac{\hbar}{2M\omega}}({a}+{a}^{\dagger})\\ =&\hbar\delta\!\omega\hat{n}-\frac{F}{2}\kappa(a+a^{\dagger})-E_{J}^{*}\left[:J_{0}\left(2\kappa\sqrt{\hat{n}}\right):+\kappa^{2}\hat{n}\right].\end{split}$$

$$:J_0\left(2\kappa \sqrt{\hat{n}}\right):=:\sum _{k=0}^{\infty }\frac{{\left(-{\kappa }^2\right)}^k}{k{!}^2}(\hat{n}{)}^k:=\sum _{k=0}^{\infty }\frac{{\left(-{\kappa }^2\right)}^k}{k{!}^2}\prod _{t=0}^{k-1}(\hat{n}-t)$$

$$F_1:=\frac{dQ}{dt}=-\frac{\tilde{\gamma }}{2}Q-\frac{\partial g}{\partial P}$$

$$F_2:=\frac{dP}{dt}=-\frac{\tilde{\gamma }}{2}P+\frac{\partial g}{\partial Q}.$$

$$J=\left(\begin{array}{cc}\frac{\partial F_1}{\partial X}& \frac{\partial F_1}{\partial P}\\ \frac{\partial F_2}{\partial X}& \frac{\partial F_2}{\partial P}\end{array}\right)$$

Full text

I Introduction

Nonlinear quantum oscillators driven by external time-periodic forces have received substantial attention in the last decade, particularly in nanomechanical systems and mesoscopic devices, see M.Dykman (2012); Poot and van der Zant (2012) and references therein. These activities strive on the one hand to obtain an understanding of fundamental properties of quantum-classical transitions Martinis et al. (1987); Grossman and Hänggi (1993); Guo et al. (2013a) and on the other hand to make use of generic nonlinear phenomena such as bifurcations Siddiqi et al. (2004, 2005); Manucharyan et al. (2007); Vijay et al. (2009); Karabalin et al. (2011); Zorin and Makhlin (2011) in actual implementations, e.g.  for enhanced sensing schemes. In many of these situations theoretical descriptions rely on the paradigmatic case of a Duffing oscillator, a driven damped oscillator with quadratic plus quartic potential surface, as long as anharmonicities are sufficiently weak to treat them perturbatively Dykman and Krivoglaz (1980); Peano and Thorwart (2006); Serban and Wilhelm (2007); Guo et al. (2010, 2011); André et al. (2012). In this way, quantum-classical transitions have been analyzed within Wigner function approaches Kheruntsyan (1999); Katz et al. (2007); Guo et al. (2010, 2011); Gevorgyan et al. (2013) and the seminal phenomenon of quantum activation has been revealed Marthaler and Dykman (2006); Dykman et al. (2011); Ong et al. (2013). Experimentally, nonlinearities appear naturally in the key element of superconducting circuits, namely, Josephson junctions. Such circuits are used in form of Cooper pair boxes to implement artificial atoms in circuit quantum electrodynamics Nakamura et al. (1999); Vion et al. (2002); Wallraff et al. (2004), but are also exploited in Josephson bifurcation amplifiers as quantum limited detectors Siddiqi et al. (2004, 2005); Vijay et al. (2009). The nonlinear quantum properties of Josephson junctions are essential in the new field of Josephson photonics, where non-classical radiation is created by a dc-voltage biased junction placed in series with a single or multiple resonators Hofheinz et al. (2011); Chen et al. (2014); Cassidy et al. (2017); Westig et al. (2017); Rolland et al. (2018). In this latter application, the light-matter coupling constant can reach easily relatively large values compared to conventional implementations, so that nonlinear effects can become extremely pronounced Rolland et al. (2018). The full nonlinearity gives rise to a wealth of quantum dynamical features far from equilibrium which may offer new platforms to access phenomena such as critical slowing down or Kibble-Zurek-type of scenarios Kibble (1980); Zurek (1996); Zurek et al. (2005). It is also important for current research which extends concepts from the pair-creation resonance Dykman et al. (1998); Zorin and Makhlin (2011); Dykman et al. (2011) to higher-order photon resonances Guo et al. (2013b); Guo and Marthaler (2016); Denisenko et al. (2016); Zhang et al. (2017) with the ultimate goal to create multi-cat states and ensembles of entangled photons Vlastakis et al. (2013); Mirrahimi et al. (2014); H.Goto (2016); Svensson et al. (2017, 2018).

Very recently, in an extension of the dynamical regime, where the Josephson bifurcation amplifier operates, the strong driving domain of a Josephson device, realized in form of a Superconducting Interference Device (SQUID), with the phase as the relevant classical degree of freedom has been studied Jung et al. (2014). As expected, higher order anharmonicities give rise to a multitude of steady state orbits with a complex bifurcation pattern. The motivation of this present work is to explore aspects of these phenomena in the quantum realm by considering a Josephson junction with its phase being quantum and subject to an external AC-current drive such that the Duffing approximation fails. By mapping the Hamiltonian from the laboratory frame to a frame rotating with the frequency of the external source, we arrive at a quasi-energy operator that is analyzed both in absence and in presence of dissipation. This allows identifying regimes, where simplifications such as the Duffing oscillator are valid. A comparison between classical fixed points and quantum steady state distributions in Fock space reveals possible classical-quantum correspondences. Note that the set-up we consider here is very different from the voltage driven case, where the simultaneous presence of a dc- and an ac-voltage gives rise to the well-known Shapiro resonances Shapiro (1963). We also consider the role of local quantum fluctuations and squeezing close and away from classical instabilities which may experimentally appear as voltage fluctuations. This may motivate further experiments to explore the wealth of the nonlinear quantum dynamics in Josephson junctions.

The paper is organized as follows: In Sec. II the model and its mapping to the rotating frame together with simplifications is introduced. The corresponding classical system is briefly addressed in Sec. III, before in Sec. IV we analyze the energy spectra of the undamped system. The remainder, Sec. V, is devoted to the steady state behavior based on a quantum master equation. Conclusions are drawn in Sec. VI.

II Time-periodic driving and the Josephson Hamiltonian

A Josephson junction with phase difference $\varphi$ and conjugate charge operator $q$ subject to a time-periodic driving is described by

[TABLE]

with $[{q},{\varphi}]=-i2e\hbar$ so that formally $q/2e$ can be seen as the momentum operator conjugate to the position $\varphi$. The driving amplitude is $F=\phi_{0}I_{d}$, where $\phi_{0}$ is the reduced flux quantum and $I_{d}$ denotes the maximum of the bias current.

Now, in a first step we transform the above Hamiltonian to a frame rotating with the frequency of the external drive, i.e., ${H}^{\text{RF}}={U}{H}{U}^{\dagger}+i\hbar{U}^{\dagger}\partial_{t}{U}$ with the unitary operator ${U}=e^{i\omega_{\text{d}}t{a}^{\dagger}{a}}$ (similar to Marthaler and Dykman (2006)). Here $a,a^{\dagger}$ denote canonical annihilation and creation operators, respectively, according to $\varphi=\kappa(a+a^{\dagger})$ with $[a,a^{\dagger}]=1$ and ground state width in the Josephson potential $\kappa^{2}=\sqrt{E_{C}/2E_{J}}$ with charging energy $E_{C}=2e^{2}/C$. In terms of a mechanical analog, this corresponds to the ground state width of a harmonic system with frequency $\omega=\sqrt{2E_{J}E_{C}}/\hbar$ and mass $M=\hbar^{2}/2E_{C}$, i.e. $\kappa=\sqrt{\hbar/2M\omega}$. Dropping all fast oscillating terms, a straightforward calculation then leads to the Hamiltonian (1) in the rotating wave approximation, i.e.,

[TABLE]

where $\hat{n}=a^{\dagger}a$ and $J_{0}$ is the zeroth-order Bessel function with the colons indicating normal ordering. Further, one has $\delta\!\omega={\omega}^{*}-\omega_{d}$ with the renormalized plasma frequency ${\omega}^{*}=\omega\,(1+{\rm e}^{-\kappa^{2}/2})/2$ and the renormalized Josephson energy $E_{\text{J}}^{*}=E_{\text{J}}{\rm e}^{-\kappa^{2}/2}$ (cf. Joyez (2013); Rolland et al. (2018)). The rotating wave approximation limits the validity of this Hamiltonian to a domain of small detuning, $|\delta\!\omega|\ll\omega_{d}$, and also sets an upper limit to the driving amplitude, see Fig. 1. The latter stems from requiring that the energy level spacing of the autonomous system sufficiently exceeds the energy scale for driving, i.e., $\kappa F<\hbar\omega$ which is equivalent to $F/E_{J}<1$.

It is instructive to consider the situation, where $\kappa\sqrt{\langle\hat{n}\rangle}\ll 1$ so that the Bessel function in (2) can be expanded. Expressed in terms of the typical steady state value for harmonic systems at resonance $\langle\hat{n}\rangle\sim(\kappa F/\hbar\gamma)^{2}$ with damping rate $\gamma$, this implies $(F/E_{J})\ll\tilde{\gamma}$, where $\tilde{\gamma}=\gamma/\omega$ is the inverse of the Q-factor of the junction at the bare plasma frequency. This way, one finds

[TABLE]

where we kept the lowest order nonlinear term in the occupation operator.

If now $\kappa\ll 1$ meaning that $E_{J}$ by far exceeds $E_{C}$, the above result can be further simplified to

[TABLE]

where $\delta\!\omega_{0}=\omega-\omega_{d}$. This expression coincides with the rotating-wave approximated form of the Duffing oscillator Hamiltonian that has been used in various contexts recently André et al. (2012); Serban and Wilhelm (2007); Guo et al. (2010, 2011). We recall that the latter describes a periodically driven nonlinear oscillator with purely quartic anharmonicity that typically is assumed to carry a parameter (length scale) that can be tuned independently from the harmonic part of the potential. This, however, is not the case for the Josephson junction, where $E_{J}$ serves as the single parameter which controls both the plasma frequency at low amplitudes as well as the anharmonicity length scale.

It is thus the domain outside the range of validity of Eqs. (3), (4) where a deep quantum regime can be found, see Fig. 1. There, ground state fluctuations of the phase (seen as a collective degree of freedom of the circuit) influence the effective junction parameters such as plasma frequency and Josephson coupling. Formally, in the limit of large $\kappa>1$ (but still $\kappa^{2}\langle n\rangle\ll 1$ which implies $F/E_{J}\ll 1$) the expression (3) can be written as a Hamiltonian of a harmonic system, i.e., ${H}_{\text{JJ, 2}}^{\text{RF}}\to\hbar(\omega/2-\omega_{d})\,\hat{n}-F\varphi$, with resonance frequency $\omega/2$.

To analyze properties of the driven Josephson system in more detail in the following, it is convenient to introduce dimensionless quantities. For this purpose, we measure all energies in units of $M\omega\delta\!\omega$ and further introduce

[TABLE]

as dimensionless parameters. In these units and expressed in dimensionless position and momentum operators $\hat{Q}$ and $\hat{P}$, respectively, according to $a=(\hat{Q}+i\hat{P})/(2\kappa)$ and with $[\hat{Q},\hat{P}]=2i\kappa^{2}$, one obtains $g\equiv H_{\text{JJ}}^{\text{RF}}/(M\omega\delta\!\omega)$ as

[TABLE]

Likewise, the quasi-energy operator for the Duffing system corresponding to (4) and also scaled with $M\omega\delta\omega$ reads

[TABLE]

where $E_{J}^{*}$ is replaced by $E_{J}$ in $\epsilon_{0}$.

III Classical system

In order to explore quantum-classical correspondences later, we here present a brief analysis of the classical steady state dynamics of the Hamiltonian (6). A full classical limit consists in replacing operators by classical phase-space variables in the rotating frame and setting $\kappa=0$ so that parameters reduce to bare values, i.e. $\delta\!\omega=\delta\!\omega_{0}$, $E_{J}^{*}=E_{J}$. Qualitatively, one identifies a limit of small $Q^{2}+P^{2}$, where an expansion of the Bessel function in the last line of (6) leads to a dominating nonlinearity of the form $\frac{1}{2}(Q^{2}+P^{2})^{2}$, so that the classical Duffing result is regained, i.e., putting $\delta\!\omega_{0}/\delta\!\omega=1$ and $\kappa=0$ in (7). In the opposite limit of large $Q^{2}+P^{2}$, the Bessel function tends to zero so that the quasi-energy turns into a harmonic function with linear tilt, where asymptotically $g(Q,P)\rightarrow\pm\infty$ (for $\epsilon\lessgtr 2$). This is illustrated in Fig. 2 for various values of the amplitude of the driving, which breaks the $Q-P$ symmetry. In the following, $f\geq 0$ is discussed for definiteness; results for $f\leq 0$ are trivially related by $Q\rightarrow-Q$.

In case of $\epsilon\gg 1$, the relevant regime is the one of small values for $Q^{2}$ and $P^{2}$ with a weak anharmonicity. It corresponds locally to a Duffing oscillator. In contrast, for $\epsilon\gtrsim 2$ the full impact of the Bessel function comes into play. Apparently, this latter one is the most interesting domain, since then in the presence of dissipation multiple steady state orbits are accessible 111 Note that we assume a conventional Josephson junction with $E_{J}>0$ (otherwise, we shift $\phi$ by $\pi$ and procced as before), so that for the Josephson potential, which becomes softer for stronger driving, multiple solutions can only possibly occur for $\epsilon>0$. . These are obtained from the stationary solutions of the equations of motion, i.e.,

[TABLE]

similar to Dykman et al. (1998).

Results for fixed points amplitudes for different values of $\epsilon>2$ are shown in Fig. 3. As expected, one sees that for $\epsilon\gg 1$ and $f/\epsilon\ll 1$ the bifurcation pattern known from the Duffing oscillator is recovered, while this picture becomes much richer for $\epsilon\gtrsim 2$ and $f/\epsilon\sim O(1)$. There, replicas of the Duffing orbits appear at larger amplitudes. A stability analysis (see Appendix B) reveals that for values of $\epsilon\ll 2$ as well as $\epsilon\gg 2$ only a single stable orbit exists around $Q=0$. Approaching $\epsilon=2$ from either side, more and more pairs of stable and unstable orbits with larger $|Q|$ successively appear (cf. Caputo et al. (2012); Lazarides and Tsironis (2013); Jung et al. (2014); Shevchuk et al. (2014); Maizelis et al. (2014) for multistable solutions). Specifically, one finds stable orbits at the minima of the quasi-energies in Fig. 2 for $Q\geq 0$ (labeled 1 and 5 and corresponding to amplitudes colored red in Fig. 3) and at the maxima at $Q<0$ (labeled 3 and 7 and colored blue). All maxima (minima) of $g(Q,P)$ for $Q>0$ ($Q<0$) correspond to unstable orbits (green-colored, labeled 2, 4, 6, 8).

IV Bare quasi-energies and eigenfunctions

We start by analyzing the quantum properties of the system by considering the bare quasi-energy operator (6). Without driving $f=0$, the quasi-energy commutes with the number operator $\hat{n}$ and is thus diagonal in the Fock state basis $\hat{n}|n\rangle=n|n\rangle$. For small $\kappa$ and large $\epsilon$ a parabola-type spectrum for the quasi-energies $E_{n}$ results [cf. Fig. 4(a)], in contrast to the regime of $\kappa\sim 1$ and $\epsilon\gtrsim 2$ [cf. Fig. 5(a)], where we see oscillations stemming from the Bessel function similar to Fig. 2 but with substantial modifications due to the normal-ordering.

Strictly without driving, $f\equiv 0$, the corresponding eigenstates are simply Fock states. A useful connection between eigenstates and classical solutions can be found, however, by taking a minute but finite driving, $f=0.01$, which breaks the $Q-P$ symmetry. In phase space, eigenstates away from maxima or minima of the quasienergy spectrum still resemble the harmonic ones: they are centered around $Q=P=0$ as shown by their Wigner density in Fig. 4(c) and (f). For states around the top of the parabola of eigenenergies [$n\approx 25$, Fig. 4(d)], where the energy spacing between neighboring states becomes comparable to the (minute) off-diagonal driving term, constructive and destructive interference leads to a phase-space localization around the first maximum at $Q<0$ [corresponding to a stable classical orbit, cf. blue line with $A\sim 5$ in Fig. 3(a)]. The same is seen for states around higher maxima in the strongly nonlinear spectrum [Fig. 5(d) and (f)] localized at positions, where higher order maxima in the classical quasi-energy appear. Despite their association with classical solutions the eigenstates can be strongly squeezed (and displaced) and carry substantial negative parts in the Wigner functions, a signature of their pronounced non-classicality.

For strong driving, energy eigenstates can no longer be indexed by $n$ (as they are not eigenstates of the number operator), but are ordered by decreasing energy (and for clarity indexed by $m$ in the following). It is illuminating to monitor how spectra evolve with growing $f$ from the ones at $f=0$. For the two $\kappa$-values shown in Fig. 4(b) and 5(b), spreading of the eigenenergies with increasing driving is observed. In particular, for $f=0$ there are clusters of energetically-close states around the minima/maxima of Fig. 4(a) and 5(a). Even a small driving leads to a strong repulsion among these near-degenerate states. For finite driving, no good quantum number beside the energy remains, so that generically anticrossings appear in the spectra.

Wigner functions at finite driving are shown in Fig. 6 for those eigenstates that can be associated with classically stable solutions for that driving strength [Fig. 6(a), cf. Fig. 2 and 3(b)]. For instance, the highest-energy eigenstate ($m=0$) corresponds to the stable solution (local minimum) with highest quasi-energy located around $(Q,\,P)=(-4,\,0)$; $m=14$ (c) corresponds to another local minimum at negative $Q\approx-16$ with strong squeezing; and Fig. 6(d) ($m=48$) displays the near-Gaussian eigenstate, which is found around the central well $Q=P=0$.

V Steady state quantum dynamics

Let us now turn to the non-equilibrium quantum dynamics. For this purpose, we employ a master equation formulation at $T=0$ for the dynamics of the reduced density operator of the Josephson system, i.e.,

[TABLE]

Steady state solutions of (9) are approached in the long time limit and can then be compared either to corresponding findings for the Duffing case or to the properties of the non-dissipative system. Thereby, $\kappa$ plays the role of a tuning parameter which allows exploring quantum-classical transitions and semiclassical fixed points. Practically, it is convenient to represent the density operator in the Fock state basis of the ground state oscillator, i.e. $\rho_{n,m}=\langle n|\rho|m\rangle$.

In the range of small values for $\kappa$ (and large $\epsilon$), one finds steady state eigenstate populations that are very close to the results of the corresponding Duffing system. Figure 7 shows results in a regime, where the equivalent classical system relaxes to one of the two stable orbits depending on the initial conditions (cf. Fig. 3), while quantum mechanically one always finds a bimodal distribution. This can be understood as mixing of the two classically stable states Guo et al. (2010, 2011); Katz et al. (2007). In fact, when one relates in a semiclassical sense the occupation number $n$ to the amplitude $A$ of the classical orbits according to $n\sim A^{2}$, one finds that the locations of the maxima of the distribution match approximately the stable steady state amplitudes in Fig. 3(a) at the respective driving, while there is no weight in the distribution at the unstable steady state amplitude. Both peaks are, however, rather broad reflecting quantum fluctuations. Correspondingly, when inspecting the full density, e.g., Re$\{\rho_{n,m}\}$ in Fig. 8, one observes two isolated domains in Fock space with off-diagonal contributions (coherences) concentrated close to the diagonal. The appearance of a checkerboard sign pattern for the large amplitude solution can be traced back to the fact, that it is located at negative Q-values in phase space. Consequently, there are no sign changes for the small-amplitude situated at positive Q. There are minor deviations in the location and height of the maxima of the population distributions when comparing the Duffing approximation to results for the full Josephson Hamiltonian. Directly setting the respective quasi-energies for $f=0$ side by side reveals increasing differences for $n$ values beyond the energy maximum, see Fig. 4(a). Overall, however, in this regime the Josephson steady state is well described by making the Duffing approximation.

The situation is very different for $\kappa\approx 1$ and smaller $\epsilon$, see Fig. 9. The full Josephson population distribution in (a) contains in the range $0\leq n\leq 90$ now three pronounced maxima as well as a very weak one around $n=4$ almost hidden behind the strong background of the peak at $n=0$ (see inset). Interestingly, maxima in the population distribution for large $n$ may exceed those of smaller $n$. The Duffing case, in contrast, exhibits only a single peak around $n=0$.

The maxima at $n=4,\,26,\,66$ are clearly related to the maxima in the bare quasi-energy spectrum, see Fig. 5(a). Comparing the location of the Josephson peaks with the classical steady state amplitudes $A$ in Fig. 3 and the topology of the quasi-energy in Fig. 2 the following picture appears: only the central stable orbit located at $Q\approx 0$ (labeled 1 with $A^{2}\approx 0.04$) and the orbits with $Q\leq 0$ (labeled 3, 7 with $A^{2}\approx 4,\,26$, respectively) can be associated with maxima in the populations $\rho_{n,n}$, while the one located at $Q>0$ (labeled 5 with $A^{2}\approx 17$) seems to be absent. Classically, some intuition about the stability of orbits against fluctuations can be gained from quasi-energy plots such as Fig. 2, where orbit 5 resides in a very shallow minimum, and bifurcation diagrams, such as Fig. 3, where a driving of $f=0.8$ is just below the critical driving strength at which the orbit disappears. While, generically, the quantum mechanical picture is very complex Marthaler and Dykman (2006); Dykman et al. (2011), such intuition about classical (in)stability nonetheless hints, whether an orbit will carry considerable weight in the steady state distributions or is washed out by fluctuations and thus explains the absence of a peak for orbit $5$ at $n\approx A^{2}\approx 17$.

The full Fock space density shown in Fig. 10 displays basically three distinct domains according to the three pronounced peaks with no discernible coherences between the distinct peaks. Only for the weak peak at $n=4$ and the central $n=0$ peak substantial off-diagonal contributions appear, indicating a type of ‘hybridization’ of two classical steady state orbits. Again the checkerboard pattern appears for solutions located at negative $Q$.

To study fluctuations around fixed points for larger amplitudes, a semiclassical analysis is illuminating (cf. Armour et al. (2013)). For this purpose, we introduce $a\rightarrow\alpha+\delta a$, with $\alpha=Ae^{-i\phi}$ being a classical fixed point with amplitude $A=\sqrt{Q^{2}+P^{2}}/2\kappa$ and phase $\phi={\rm arctan}(\frac{P}{Q})$ while $\delta a$ ($\delta a^{\dagger}$) denote local annihilation (creation) operators $[\delta a,\delta a^{\dagger}]=1$. Accordingly, one considers a displaced reduced density operator $\sigma=D^{\dagger}(\alpha)\rho D(\alpha)$ with the displacement operator $D(\alpha)=\exp(\alpha a^{\dagger}-\alpha*a)$ and neglects higher than second order terms in the quasi-energy operator (6), i.e. $g\to g_{2}$, with

[TABLE]

with the Bessel functions $J_{0}$ and $J_{2}$. In particular, one has $\langle n\rangle=|\alpha|^{2}=A^{2}$ and $\langle\delta a\rangle=0$ and $\sigma$ follows from (9) with $g$ replaced by $g_{2}$. Accordingly, the Heisenberg equations for the displaced operators read

[TABLE]

where the noise operators $a_{\rm in},a_{\rm in}^{\dagger}$ obey: $\langle a_{\rm in}\rangle=\langle a_{\rm in}^{\dagger}\rangle=0$, $\langle a_{\rm in}(t)a_{\rm in}(t^{\prime})\rangle$=$\langle a_{\rm in}^{\dagger}(t)a_{\rm in}^{\dagger}(t^{\prime})\rangle$=$\langle a_{\rm in}^{\dagger}(t)a_{\rm in}(t^{\prime})\rangle=0$ and $\langle a_{\rm in}(t)a_{\rm in}^{\dagger}(t^{\prime})\rangle=\delta(t-t^{\prime})$. Note that here $\alpha$ is a function of the parameters $f$ and $\epsilon$. This way, one finds for the number fluctuations the expression

[TABLE]

The eigenvalues of the matrix in (11) determine the stability of the fixed point solutions and are obtained as

[TABLE]

In the harmonic limit, this reduces to $\lambda_{\pm}=-\gamma/2\pm 2i$ and fixed points are always stable. Outside, a change of stability occurs when $\lambda_{+}=0$ by tuning $f$ (and thus $|\alpha|$) and/or $\epsilon$ which implies that number fluctuations (or energy fluctuations) according to (V) diverge. Figure 11 illustrates this behavior for fixed points corresponding to maxima at $n=0$ and $n=26$ in Fig. 9; in Fig. 3(b) these may roughly be associated with the classical orbits of least amplitude (label 1 in Fig. 3, red line in Fig. 11) and with $A\approx 5$ (label 5 in Fig. 3, blue line in Fig. 11). For the fixed point emerging from the harmonic one with increasing driving, number fluctuations start from the result for a purely coherent state $\delta n_{2}=1$ and remain close to this value for a broad range of driving strengths. Close to the instability point $\delta n_{2}$ sharply grows indicating strong quadrature squeezing. The classical fixed point with higher amplitude exhibits two instability points and possesses only a narrow range of driving, where $\delta n_{2}$ is close to the coherent state value; outside strong quadrature squeezing dominates. Qualitatively, this originates from the fact that around a bifurcation at $f_{c}$ with amplitude $A_{c}$ number fluctuations scale as $\delta n_{2}\sim A_{c}/|f-f_{c}|^{\beta}$ with $|A-A_{c}|\propto|f-f_{c}|^{\beta}$ with $\beta<1$. In the present case, one has $\beta\approx 1/2$ (cf. Fig. 3) which implies that the range of $f$-values with $\delta n_{2}>1$ becomes very broad for larger amplitudes $A_{c}$. However, between these regions narrow domains are found, where the system reduces again to a coherent state with $\delta n_{2}=1$. This is a consequence of $J_{2}(2\kappa|\alpha|)=0$ in (V).

VI Conclusion

The quantum dynamics of a Josephson junction driven by an external time-periodic force coupling to its phase is investigated in the regime, where driving is so strong that the conventionally employed reduction to a Duffing oscillator fails. Since there is no static force (dc-current) applied, known resonances in form of Shapiro steps do not exist Shapiro (1963). Starting from the Hamiltonian in the laboratory frame, we investigate the system in a rotating frame and reveal under which conditions simplifications such as, for example, the Duffing approximation are valid. Ground state fluctuations of the phase induce a renormalization of the bare Josephson coupling, a phenomenon that has been recently detected in a somewhat different setting Joyez (2013); Rolland et al. (2018). The junction’s response is thus resonant only if the driving frequency matches the renormalized plasma frequency of the junction. Further, the quantum dynamics in the rotating frame and in presence of dissipation is studied within a master equation approach. Results are discussed in comparison to classical steady state orbits and quantum results for the Duffing system. The main difference to the latter is the occurrence of a complex pattern of multiple fixed points in the classical limit leading quantum mechanically to multiple maxima in the population distribution. Coherences between classically well separated orbits are absent while they are found to be strong for closely adjacent orbits. While classically the system locks to specific fixed points with narrow domains, where changes of stabilities due to bifurcations occur, quantum mechanically, substantially enhanced number fluctuations (or energy fluctuations) corresponding to strong quadrature squeezing are seen over broad ranges of driving strengths. In very narrow domains, the system collapses to coherent states again with strongly suppressed fluctuations though. In the laboratory frame large energy fluctuations induced by time-dependent currents manifest themselves as large fluctuations of ac-voltages and thus may be accessible experimentally. Our study sheds light on Josephson physics in a regime that has not been explored yet in depths, in some contrast to the fact that Josephson junctions serve as paradigmatic nonlinear elements for mesoscopic superconducting circuits. It may allow to elucidate transitions between the well-known Duffing system and higher order nonlinear oscillators and may thus contribute to the ongoing experimental activities to explore nonlinear quantum dynamics far from equilibrium.

acknowledgments

We are grateful to L.Guo for discussions and for financial support by the German Science Foundation through SFB/TRR 21 and the Center for Integrated Quantum Science and Technology (IQST).

Appendix A Transformation of the Josephson Hamiltonian into the rotating frame

We start with the Josephson Hamiltonian

[TABLE]

We transform the above Hamiltonian to a frame rotating with the frequency of the external drive, i.e., ${H}^{\text{RF}}={U}{H}{U}^{\dagger}+i\hbar{U}^{\dagger}(\partial_{t}{U})$ with the unitary operator ${U}=e^{i\omega_{\text{d}}t{a}^{\dagger}{a}}$. Here $a,a^{\dagger}$ denote canonical annihilation and creation operators, respectively, according to $\hat{Q}=\kappa(a+a^{\dagger})$ with $[a,a^{\dagger}]=1$ and ground state width of the Josephson potential $\kappa^{2}=\sqrt{E_{C}/2E_{J}}$ with charging energy $E_{C}=2e^{2}/C$.

By applying the unitary transformation operators to the nonlinear parts of the Josephson Hamiltonian

[TABLE]

we used the Baker-Campbell-Hausdorff-formula $e^{X+Y}=e^{X}e^{Y}e^{-\frac{[X,Y]}{2}}$ and the Bessel function $J_{0}(z)=\sum_{k=0}^{\infty}(-1)^{k}(\frac{z}{2})^{2k}\frac{1}{(k!)^{2}}$. After the transformation we find the Hamiltonian in the rotating frame

[TABLE]

To investigate higher orders it is more convenient to write the Bessel function

[TABLE]

including a product of $\hat{n}$, where a normal ordering is not necessary any more.

Appendix B Classical stability analysis

To perform a stability analysis of the steady state solutions of

[TABLE]

and

[TABLE]

we compute the extrema for the Jacobi Matrix

[TABLE]

using $\frac{\partial}{\partial X}(\frac{J_{1}(\Delta)}{\Delta})=\frac{\partial}{\partial\Delta}(\frac{J_{1}(\Delta)}{\Delta})\frac{\partial\Delta}{\partial X}=-\frac{J_{2}(\Delta)}{\Delta}\frac{\partial\Delta}{\partial X}$. Stable solutions are related to a negative real part of the eigenvalues, unstable ones to a positive real part of the eigenvalues.