Discrete time symmetry breaking in quantum circuits: exact solutions and tunneling

TL;DR

This paper provides exact solutions and analysis of discrete time symmetry breaking in quantum circuits, revealing how nonlinearities and anharmonicities influence stability, tunneling, and multistability in superconducting quantum systems.

Contribution

It offers the first exact analytic solutions for steady states and tunneling times in discrete time quantum symmetry breaking, incorporating effects of anharmonicities in superconducting circuits.

Findings

Anharmonic terms increase tunneling rates.

Exact solutions for steady states and tunneling times are derived.

Superconducting circuits exhibit regimes of tristability and bistability.

Abstract

We discuss general properties of discrete time quantum symmetry breaking in degenerate parametric oscillators. Recent experiments in superconducting quantum circuit with Josephson junction nonlinearities give rise to new properties of strong parametric coupling and nonlinearities. Exact analytic solutions are obtained for the steady-state of this single-mode case of subharmonic generation. We also obtain analytic solutions for the tunneling time over which the time symmetry-breaking is lost above threshold. We find that additional anharmonic terms found in the superconducting case increase the tunneling rate, and can also lead to new regimes of tristability as well as bistability. Our analytic results are confirmed by number state calculations.

Equations

$$H=\hslash \sum _{n=1}^6\sum _{k=1}^2{H}_k^{(n)},$$

$$H_k^{(1)}$$

$$H_k^{(2)}$$

$$H_k^{(3)}$$

$$H_k^{(4)}$$

$$H_1^{(5)}$$

$$H_1^{(6)}$$

$$H_{rev}=\hslash \sum _{n=3}^6{H}^{(n)}.$$

$$a_k^{†}\left(t\right)=a_k^{†}(0)e^{-ik\omega t},$$

$$H^{SR}=H_{rev}-\sum _{k}k\hslash \omega a_k^{†}{a}_k.$$

$$\dot{\rho }=-\frac{i}{\hslash }\left[H^{SR},\rho \right]+\sum _{k,j>0}\frac{1}{j}{\gamma }_k^{(j)}{\mathcal{L}}_k^{(j)}\left[\rho \right].$$

$${\mathcal{L}}_k^{(j)}\left[\rho \right]=2\hat{O}\rho {\hat{O}}^{†}-\rho {\hat{O}}^{†}\hat{O}-{\hat{O}}^{†}\hat{O}\rho +2n_k^{(j)}\left[\left[\hat{O},\rho \right],{\hat{O}}^{†}\right].$$

$$\hat{\rho }=\int \int d\mu \left(\mathbf{\alpha },{\mathbf{\alpha }}^{+}\right)P\left(\mathbf{\alpha },{\mathbf{\alpha }}^{+}\right)\hat{\Lambda }\left(\mathbf{\alpha },{\mathbf{\alpha }}^{+}\right).$$

$$⟨a_j^{†}\dots a_k⟩=\int \int d\mu \left(\mathbf{\alpha },{\mathbf{\alpha }}^{+}\right)P\left(\mathbf{\alpha },{\mathbf{\alpha }}^{+}\right)\left[{\alpha }_j^{+}\dots {\alpha }_k\right].$$

$$\frac{\partial P}{\partial t}=\sum _{n=1}^6\sum _{k=1}^2{\mathcal{D}}_k^{(n)}P,$$

$${\mathcal{D}}_k^{(1)}$$

$${\mathcal{D}}_1^{(2)}$$

$${\mathcal{D}}_k^{(3)}$$

$${\mathcal{D}}_k^{(4)}$$

$${\mathcal{D}}_1^{(5)}$$

$${\mathcal{D}}_1^{(6)}$$

$$\frac{\partial P}{\partial t}=\left[-{\partial }_{\mu }{A}_{\mu }+\frac{1}{2}{\partial }_{\mu }{\partial }_{\nu }{D}_{\mu \nu }\right]P.$$

$${\gamma }_k={\gamma }_k^{(1)}+i{\Delta }_k,g_1={\gamma }_1^{(2)}+i\chi .$$

$$\vec{A}\left(\vec{\alpha}\right)=\left[\begin{array}[]{c}\mathcal{E}_{1}-\gamma_{1}\alpha_{1}-g_{1}\alpha_{1}^{+}\alpha_{1}^{2}+\kappa\alpha_{1}^{+}\alpha_{2}\\ \mathcal{E}_{1}^{*}-\gamma_{1}^{*}\alpha_{1}^{+}-g_{1}^{*}\alpha_{1}\alpha_{1}^{+2}+\kappa\alpha_{1}\alpha_{2}^{+}\\ \mathcal{E}_{2}-\gamma_{2}\alpha_{2}-\kappa\alpha_{1}^{2}/2\\ \mathcal{E}_{2}-\gamma_{2}^{*}\alpha_{2}^{+}-\kappa\alpha_{1}^{+2}/2\end{array}\right].$$

$$\underline{\underline{D}}\left(\vec{\alpha}\right)=\left[\begin{array}[]{cc}\underline{D}_{1}\\ &\underline{0}\end{array}\right],$$

$${\underline{D}}_1$$

$$\frac{d\vec{\alpha }}{dt}=\vec{A}\left(\vec{\alpha }\right)+\underline{\underline{B}}\left(\vec{\alpha }\right)\vec{\zeta }\left(t\right),$$

$$⟨{\zeta }_i\left(t\right){\zeta }_j\left(t^{\prime }\right)⟩={\delta }_{ij}\delta \left(t-t^{\prime }\right).$$

$$\underline{\underline{B}}\left(\vec{\alpha}\right)=\left[\begin{array}[]{cc}\underline{B}_{1}\\ &\underline{0}\end{array}\right],$$

$${\underline{B}}_1$$

$$\hat{\rho }={\int }_{\mathcal{C}}d\overrightarrow{{\alpha }_1}{P}_1\left({\vec{\alpha }}_1\right)\int d\overrightarrow{{\alpha }_2}{P}_2\left({\vec{\alpha }}_2|{\vec{\alpha }}_1\right)\hat{\Lambda }\left(\mathbf{\alpha },{\mathbf{\alpha }}^{+}\right).$$

$${\dot{\alpha }}_2={\mathcal{E}}_2-{\gamma }_2{\alpha }_2-\kappa {\alpha }_1^2/2.$$

$${\alpha }_2^{(a)}=\frac{{\mathcal{E}}_2-\kappa {\alpha }_1^2/2}{{\gamma }_2}.$$

$$P_2\left({\vec{\alpha }}_2|{\vec{\alpha }}_1\right)=\delta \left({\vec{\alpha }}_2-{\vec{\alpha }}_2^{(a)}\right).$$

$${\hat{\rho }}_1=\int \int d\mu \left(\vec{\alpha }\right)P_1\left(\vec{\alpha }\right){\hat{\Lambda }}_1\left(\vec{\alpha }\right).$$

Full text

††: New J. Phys.

Keywords: time symmetry breaking, quantum optics, quantum tunneling, phase space methods, degenerate parametric oscillation.

1 Introduction

Quantum time symmetry breaking is a widespread phenomenon in non-equilibrium quantum optics and superconducting quantum circuits. For quantum clocks and lasers, time and phase are inter-related. Therefore, quantum time symmetry breaking occurs in the electromagnetic phase. This is implicit in the use of coherent states, which have a well-defined phase, to describe lasers [1, 2]. Number state and coherent state descriptions are both complete descriptions in quantum mechanics. However, coherent state expansions allows one to recognize more readily that time translational symmetry is broken through an observation of the phase.

Discrete time quantum symmetry breaking takes place in intra-cavity subharmonic generation, otherwise known as degenerate parametric oscillators [3]. For quantum optical systems, an exact solution for the steady-state quantum density matrix [4] is known for the special case of a single-mode system with no detunings or anharmonicity. From this, one can calculate quantum tunneling between time phases [5]. Schrödinger cats might also seem possible in the steady-state [6]. However, this is not the case [7], although transient Schrödinger cat formation is possible [7, 8, 9] for strong enough coupling. Tunneling in these systems demonstrates the existence of long range time order, which has been confirmed in optical experiments [10] that measured extremely narrow subharmonic line-widths.

Classical spontaneous symmetry breaking is widespread in quantum optics and related fields, including pattern formation with translational and cylindrical symmetry breaking [11], noncritical squeezing with rotational symmetry breaking [12], and super-solid formation with continuous translational symmetry breaking [13]. Similar to normal crystals which have a repeating pattern in space, systems with quantum time symmetry breaking repeat themselves in time. In the case of spontaneous discrete time-translation symmetry breaking oscillation takes place at a fraction of the frequency of some periodic driving force [14]. Here we consider the extent to which intra-cavity subharmonic generation leads to discrete time symmetry breaking. This requires an investigation of the time-scale for restoration of the original symmetry from a symmetry-broken state.

As well as transient macroscopic superpositions and tunneling, subharmonic generators can also generate space-time ordering [15]. These phenomena are often termed as time crystals [16, 14], since they combine discrete time symmetry breaking with spatial ordering. Time crystals have also been observed in spin systems [17, 18] and Bose-Einstein condensates [19] experimentally.

As an example of subharmonic generators, degenerate parametric oscillation has been investigated for a long time in optics [3, 4, 5, 7, 8, 15]. These devices convert a high-frequency input optical mode into two equal frequency modes, whose frequency is half of the input mode’s, by means of the parametric nonlinearity. Below threshold, they are quantum squeezed state generators. These are used to reduce quantum noise in gravity-wave detectors [20, 21]. Networks of above threshold parametric devices are now being used as analog quantum computers for NP-hard optimization [22, 23].

The most general case of single mode degenerate parametric oscillation involves arbitrary detunings, nonlinear losses and anharmonic nonlinearities in the fundamental mode. This is characteristic, for example, of superconducting cavity experiments [24]. In these investigations, two superconducting microwave cavities are coupled through a Josephson junction in a bridge transmon configuration. Mode one, also termed “the storage”, holds the mode with time symmetry breaking, and is designed to have minimal single-photon dissipation. Mode two, “the readout”, is over-coupled to a transmission line, and it removes entropy from the storage.

These experiments used Josephson junctions to generate a coupling that exchanges pairs of photons in the storage with single photons in the readout. Single-photon dissipation for both modes $a_{1}$ and $a_{2}$ are included. In these systems, as well as other parametric devices with small mode volumes, there is a large anharmonic nonlinearity, and possible detunings. These additional effects change the physics, and modify both the steady-state quantum behavior and the tunneling rates compared to previous studies.

Quantum tunneling has also been studied in single-mode nonlinear photonic resonators, both theoretically [25, 5, 26] and experimentally [27]. In our more general case there are several parameters to control the tunneling rates, which were not included in previous work. We study the mean-field solutions in detail, demonstrating the existence of a universal phase diagram with monostable, bistable and tristable phases. Including quantum fluctuations, we obtain the full quantum distribution in all cases, and focus on the limiting case of pure quantum tunneling at the midpoint of the hysteresis curve in the bistable region. This is the most long-lived tunneling region. We show that in this regime there is an exponentially long tunneling time at large driving. It is this case where quantum time symmetry breaking has the longest lifetime before the discrete symmetry is restored.

We first obtain the exact steady-state solution expressed in form of a complex potential function. The regime where damping exceeds coupling at the single-photon level is studied in detail. We find that quantum tunneling occurs in this parameter region, slowly equilibrating a system with discrete time-symmetry breaking to the steady-state. The tunneling time is obtained analytically within a novel, complex potential-barrier approximation in the complex P-representation valid at large particle number. The results agree with those numerically obtained in the number-state basis for relatively small photon number. This opens the way to studying quantum tunneling in multi-mode non-equilibrium systems, as in the coherent Ising machine [22, 23].

The paper is organized as follows. In section 2 we introduce the Hamiltonian of the two-mode quantum circuit system. Taking the driving, damping and nonlinearities into account, the master equation for the time evolution of the system is derived. We use the generalized P-representation to obtain a conditional Fokker-Planck equation, and then employ adiabatic elimination of the strongly damped pump mode to obtain a simpler one-mode Fokker-Planck equation. In section 3 the mean-field limit is analysed, showing the existence of stable, bistable and tristable phases.

The quantum steady-state and quantum tunneling are studied in section 4, demonstrating that a novel complex manifold method can be used to calculate tunneling rates. In section 5, agreement between the analytic results and numerical calculations using number states is obtained. Finally, section 6 gives an outlook and our conclusions, with proofs of the complex tunneling rate method given in an A.

2 Time-evolution of the system

2.1 General model Hamiltonian

Quantum optical and quantum circuit physics are closely related. The principal difference is that quantum circuits operate at much lower temperatures, and at microwave rather than optical frequencies. To treat both cases, we consider a general model for nonlinear interactions and damping of two coupled bosonic modes of an open system. We include driving, damping and both coherent and dissipative nonlinearities in the degenerate parametric oscillator.

We suppose that $a_{k},a_{k}^{\dagger}$ are the $k$-th mode annihilation and creation operators of modes at two different frequencies $\omega_{k}$ in coupled resonant cavities. There is an overall quantum Hamiltonian given by

[TABLE]

where $H_{k}=\sum_{n}H_{k}^{(n)}$ is the Hamiltonian for the $k$-th mode, with driving, damping, and nonlinear terms. We define $H^{(n)}=\sum_{k}H_{k}^{(n)}$ as the sum over modes for the $n-th$ type of interaction. The detailed structure of these terms is as follows:

[TABLE]

Here $\hat{\Gamma}_{k}$ are the external reservoir coupling operators with reservoir Hamiltonians $H_{k}^{R}$, and $\mathcal{E}_{k}$ are the envelope amplitudes of external coherent driving fields at angular frequency $k\omega$ for the $k$-th mode. The nonlinear parameters are $\kappa$ for subharmonic generation and $\chi$ for anharmonicity. We suppose $\omega_{2}\simeq 2\omega_{1}\simeq 2\omega$, so the system can be externally driven simultaneously at fundamental and subharmonic frequencies, although we include detunings as well.

In many typical experiments [24], the inputs are in the higher frequency mode, $a_{2}$. Thus, the subharmonic driving field is zero for such experiments. We will keep this term in the general Hamiltonian and exact solutions, for completeness, but set it to zero in the tunneling calculations below.

2.2 Master Equation

The system Hamiltonian $H_{rev}$ is defined as the reversible part of the Hamiltonian without reservoir couplings, so that:

[TABLE]

The rotating frame system Hamiltonian $H^{S}$ is obtained by subtracting the driving frequency terms. This is used to give a picture such that only slowly varying behavior is retained in the state equation, while the operators evolve at their respective driving frequencies. Therefore, we define an interaction picture such that:

[TABLE]

and evolve the density matrix using the subtracted Hamiltonian, $H^{SR}$. From now on we let $a_{k}^{\dagger}=a_{k}^{\dagger}(0)$, and define:

[TABLE]

This transformation has the effect of changing the mode frequencies in the Hamiltonian to relative detunings, so that $\omega_{k}\rightarrow\Delta_{k}=\omega_{k}-k\omega$. The resulting master equation [28, 29] for the quantum density matrix $\rho$, on tracing over all the reservoirs in the Markovian approximation, is:

[TABLE]

Here $\gamma_{k}^{(j)}$ are the linear and nonlinear amplitude relaxation rates for $j$-boson relaxation in the $k$-th mode, and $n_{k}^{(j)}$ are the corresponding thermal occupations of the reservoirs. These are $n_{k}^{(j)}=1/\left[\exp\left(j\hbar\omega_{k}/k_{B}T\right)-1\right]$, where $j\omega_{k}$ is the resonant frequency of the $k$-th mode reservoir for the $j$-photon damping process, and $k,j=1,2$. Note that either resonant mode could have single-photon and/or two-photon losses [30], but for simplicity we only include nonlinear losses for $k=1$. Higher order decoherence is also possible [31], but not treated here.

We set $\hat{O}=\hat{a}_{k}^{j}$ to describe general $j$-photon damping in the $k$-th mode, giving the the master equation super-operator for decoherence as

[TABLE]

We will put $n_{2}^{(1)}=n_{k}^{(2)}=0$, for simplicity, since these reservoirs have at least twice the frequency of the fundamental reservoirs, and will have lower thermal occupations, which we neglect. The energy relaxation rate in each mode for single-particle decay is $2\gamma_{k}^{(1)}$. In the exact solutions presented in later sections, we also assume that $n_{1}^{(1)}=0$. This allows us to investigate the important exact quantum tunneling solutions in the low-temperature limit, although finite temperature effects in the fundamental reservoirs are included in the next section for completeness. We include thermal occupation for single-photon processes at this stage, to show how this alters the resulting equations in high temperature cases where such effects are important.

2.3 Fokker-Planck equation

We now introduce the generalized P-representation developed in the reference [32]. This expands the quantum density matrix in terms of a complete operator basis $\hat{\Lambda}\left(\bm{\alpha},\bm{\alpha}^{+}\right)$, and a P-representation $P\left(\bm{\alpha},\bm{\alpha}^{+}\right)$ so that:

[TABLE]

The operator basis uses off-diagonal coherent state projection operators, and has the form, $\hat{\Lambda}\left(\bm{\alpha},\bm{\alpha}^{+}\right)=\left|\bm{\alpha}\right\rangle\left\langle\bm{\alpha}^{+*}\right|/\left\langle\bm{\alpha}^{+*}\right|\left.\bm{\alpha}\right\rangle$, where $\left|\bm{\alpha}\right\rangle$ is a two-mode coherent state. This includes a normalizing factor so that the P-distribution integrates to unity. The integration measures $d\mu$ can be chosen as either a volume measure on a full eight-dimensional complex space over the combined vector $\vec{\alpha}=\left(\bm{\alpha},\bm{\alpha}^{+}\right)$, or as a contour integral on a complex manifold, which is explained in detail below. The advantage of this approach is that the resulting Fokker-Planck equation is exact without truncation. This is not the case, for example, with the Wigner representation. In the quantum-dominated regime with a strongly damped high-frequency mode, we show that there is an exact solution for the equilibrium steady-state of the resulting single-mode Fokker-Planck equation.

All normally ordered quantum correlation functions are moments of the distribution, since:

[TABLE]

Using standard operator identities, the resulting Fokker-Planck equation has the form that:

[TABLE]

where we introduce the notation that $\mathcal{D}_{k}^{(n)}$ is a differential operator acting on the P-function. The derivative operators are: $\partial_{k}\equiv\partial/\partial\alpha_{k}$, $\partial_{k}^{+}\equiv\partial/\partial\alpha_{k}^{+}$, and the individual terms involved that correspond to each Hamiltonian coupling are:

[TABLE]

Here h.c. denotes a term generated from the hermitian conjugate operator identities, in which coefficients are conjugated, and $\alpha_{k}\rightarrow\alpha_{k}^{+}$. Due to the freedom to make phase rotations in defining the mode operators, we can define $\mathcal{E}_{2},\kappa$ to be real parameters without loss of generality.

By combining terms, we obtain the following Fokker-Planck equation, where $\mu=1,\ldots 4$, $\vec{\partial}=\left(\partial_{\bm{\alpha}},\partial_{\bm{\alpha}^{+}}\right)$, and an Einstein summation convention is used to sum over repeated indices:

[TABLE]

In this form it is convenient to introduce complex single particle decay rates $\gamma_{k}$ and two-particle decay rates $g_{1}$, so that these parameters can be combined into complex rate terms:

[TABLE]

With these definitions, the combined deterministic or drift term becomes:

[TABLE]

The corresponding combined diffusion coefficient is then:

[TABLE]

where the individual mode diffusion sub-matrices are:

[TABLE]

while $\Gamma_{1}\equiv 2\gamma_{1}^{(1)}n_{1}^{(1)}$ is the thermal noise coefficient.

2.4 Stochastic Equations

If all modes have strong enough damping, so that all boundary terms vanish in the Fokker-Planck equation, there is a corresponding stochastic equation for the positive P-representation [32], which can be written in a combined vector form as:

[TABLE]

where the Gaussian noise term $\zeta$ has a vanishing mean, and the only nonzero correlations are:

[TABLE]

The corresponding combined stochastic coefficient is then:

[TABLE]

where the individual mode noise sub-matrices are:

[TABLE]

In the next section, we will focus on the steady-state solutions in the zero temperature limit, in order to understand the steady-state properties of maximal quantum coherence. Although these stochastic equations are useful when the damping rates of both modes are large compared to the nonlinearities, they have no known analytic solutions. In addition, these stochastic equations can have boundary terms at strong coupling, leading to instabilities. For this reason, we turn next to a hybrid method. This allows us to derive a solution for the complex P-representation of the subharmonic mode, which allows us to analytically calculate the tunneling rate.

2.5 Hybrid representation

We will consider the case where the second harmonic mode is strongly damped, and the first harmonic mode is not, as in many experiments. This will be treated in a hybrid measure, where the second harmonic mode is treated stochastically, while the first harmonic is expanded on a complex manifold. In this case, we extend methods used previously through defining a conditional $P_{2}$ distribution for the second harmonic mode that depends on the amplitude of the first mode, so that:

[TABLE]

Since this $\vec{\alpha}_{2}$ mode is strongly damped, it can be readily solved on the relevant short time-scales. This is equivalent to a standard characteristic function solution of a first order partial differential equation:

[TABLE]

In the limit of $\gamma_{2}^{(1)}\gg\gamma_{1}^{(1)}$, the second harmonic mode is rapidly damped to a deterministic solution,

[TABLE]

There is a similar equation for $\alpha_{2}^{+}$, which means that in the adiabatic limit

[TABLE]

For simplicity, we assume that in this strong damping limit the corresponding detuning $\Delta_{2}$ is negligible, and therefore $\gamma_{2}$ is treated as a real parameter. The details of adiabatic elimination in the full quantum theory are treated in the next section.

2.6 Quantum adiabatic elimination

We will treat the quantum effects in the adiabatic limit, but with quantum noise included. We now introduce a reduced P-representation obtained by tracing over the second-harmonic mode, so that $\hat{\rho}_{1}=\Tr_{2}\left(\hat{\rho}\right)$. If we expand the quantum density matrix in terms of a single-mode operator basis $\hat{\Lambda}_{1}\left(\vec{\alpha}\right)$, and a single-mode P-representation $P\left(\vec{\alpha}\right)$, where $\vec{\alpha}=\vec{\alpha}_{1}$, one then obtains:

[TABLE]

The operator basis uses coherent state projection operators, as before, but with a simpler form, $\hat{\Lambda}_{1}\left(\vec{\alpha}\right)\equiv\left|\alpha\right\rangle\left\langle\alpha^{+*}\right|/\left\langle\alpha^{+*}\right|\left.\alpha\right\rangle$. All normally ordered single-mode quantum correlation functions are moments of the distribution, since:

[TABLE]

With this approach, the elimination of the $\alpha_{2}$ amplitude results in a single-mode Fokker-Planck equation,

[TABLE]

Here $\gamma\equiv\gamma_{1}$ and the combined effective nonlinear loss is:

[TABLE]

This leads to a combined nonlinear coefficient $g$, where:

[TABLE]

We have also defined

[TABLE]

with $\mathcal{E}\equiv\kappa\mathcal{E}_{2}/\gamma_{2}$. Physically, $\bar{\mathcal{E}}\left(\alpha\right)$, together with its conjugate, $\bar{\mathcal{E}}^{*}\left(\alpha^{+}\right)$, are input fields that include depletion. The notation $hc$ indicates hermitian conjugate terms obtained by the replacement of $\alpha\rightarrow\alpha^{+}$, and the conjugation of all complex parameters. In the present case that $\Delta_{2}=0$, we have a real $\gamma_{2}$ and a real $\mathcal{E}$, but $g$ is still generally complex because of the nonzero anharmonic nonlinearity.

It is simpler for the detailed analysis of this problem to use dimensionless parameters, which we define as follows:

[TABLE]

We also introduce a relative phase, $e^{i\theta}=g/|g|=n/\epsilon$. This gives a Fokker-Planck equation in a more universal form, as:

[TABLE]

Physically, time is scaled relative to the two-photon driving rate, while $c$ is a complex dimensionless single-photon loss and detuning. The important scaling parameter $n$ is the photon number at which saturation of the mode occupation occurs due to the nonlinear losses.

2.7 Equivalent Hamiltonian

After adiabatically eliminating the second harmonic mode [4, 24], we obtained a new Fokker-Planck equation, as described above. This corresponds exactly to the dynamics for the fundamental quantum system with $a\equiv a_{1}$, governed by the adiabatic Hamiltonian,

[TABLE]

together with a single-photon loss $\gamma^{(1)}=\gamma_{1}^{(1)}$ and an effective two-photon loss $\gamma^{(2)}$ which is defined above.

The master equation of the density matrix $\rho_{1}$ that is equivalent to the new Fokker-Planck equation given above is then obtained as

[TABLE]

where we have used a zero-temperature limit. Similar systems have been studied before where quantum squeezing and bifurcation have been found [33, 34, 35, 36], but here we focus on the issue of quantum tunneling.

Equation (2.38) can be solved numerically using a number-state expansion provided the photon number is not too large. However, to gain more insight from analytic results, it is also possible to use the P-representation directly, given the single-mode Fokker-Planck equation (2.31). In the following sections, we will employ both methods.

3 Mean field theory predictions

Before carrying out the full quantum theory calculations, we obtain the predictions of mean field theory for the degenerate parametric oscillator, where the system is described by the equation (2.19) but without fluctuations. The resulting deterministic equations are a good approximation when the system has negligible quantum and thermal noise. Since the noise terms are generated by nonlinearities, this implies that the nonlinearity at the single photon level is much smaller than the damping. From the dimensionless Fokker-Planck equation, equation (2.36), the mean-field limit corresponds to $n\rightarrow\infty$.

Earlier work on the mean field theory in single-mode driven systems included a degenerate parametric oscillator without anharmonicity (i.e. $\chi=0$) [3], and a purely anharmonic system without parametric terms [37]. The steady state solutions for the deterministic equations were found and their stability was studied by considering the behavior of small perturbations around these steady state solutions. This work showed that the subharmonic mode has a steady state with zero mean amplitude when the pump amplitude is below a certain critical value. Above this critical value, the zero mean subharmonic mode amplitude is no longer a stable steady state solution. Rather, the mode has a bistable solution where the two steady states have an equal amplitude but with a phase difference of $\pi$. We now show that the situation when there is both anharmonic and detuning terms added is more complex than this.

3.1 Mean field dynamics

In the mean field case, one sets $\alpha_{k}^{+}=\alpha_{k}^{*}$, and treats the dynamical equation of

[TABLE]

where the drift terms are

[TABLE]

Here the single photon damping rates $\gamma_{1}$ and $\gamma_{2}$ are generally complex. For $\gamma_{2}^{(1)}\gg\gamma_{1}^{(1)}$, one can adiabatically eliminating the mode $\alpha_{2}$, just as in the Fokker-Planck approach of equation (2.27).

Next, we set $\mathcal{E}_{1}=0$ in the mean field calculation, to correspond to the bistable situation of most interest here. The mean-field equation for the mode amplitude $\alpha$ is then

[TABLE]

To simplify this, we introduce the dimensionless parameters of (2.35), which normalize the damping by the driving. These equations can also be obtained directly from the scaled Fokker-Planck equation (2.36). On scaling, the second-order derivative terms can be neglected in the mean-field limit of $n\gg 1$, leading to the mean-field dynamical equations:

[TABLE]

3.2 Mean-field stationary states

To obtain the mean-field stationary values, we set $\dot{\beta}=\dot{\beta}^{*}=0$, so that

[TABLE]

This has the solutions $\beta=0$, together with the solutions of the quadratic

[TABLE]

The quadratic has two roots for the intra-cavity intensity:

[TABLE]

where:

[TABLE]

Since this is an intensity, a negative solution or complex solution is not possible. A real solution clearly requires $|\Im(c)|\leq 1$. After calculating the corresponding amplitude in each case of a non-negative intensity, we find that there are in general three types of stationary solutions, which we classify below.

Each may be stable or unstable, as we see in the next subsection:

Vacuum solutions. These have $\beta^{(1)}=0$, which is always a stationary solution. 2. 2.

Positive branch solutions. Taking the upper sign in the intensity equation (3.7), this has a positive solution either if $\Re(c)=\left(c+c^{*}\right)/2<0$ and $|\Im(c)|\leq 1$, or else if there is a large enough driving field so that $|c|<1\,$. The corresponding amplitude is:

[TABLE]

This is the only possible solution in the mean-field limit, if one has no detunings. We show in the next section that this corresponds to a stationary point (4.17) of the quantum P-representation distribution. 3. 3.

Negative branch solutions. Next, investigating the lower sign case of 3.7, this has a positive solution if both $\Re(c)<0$, and there is a small enough driving field so that $|c|>1>|\Im(c)|\,$, with a resulting amplitude of:

[TABLE]

Since we are considering the case of $\Delta_{2}=0$ , we have the equation:

[TABLE]

Thus $\Re(c)<0$ is possible if we take $\chi\Delta_{1}<-\gamma^{(2)}\gamma_{1}^{(1)}$. This can occur with an anharmonic term at large detunings, even if $\gamma_{1}^{(1)}\gg\chi$ as required for the mean-field limit. As we see below, this occurs as an unstable branch in degenerate parametric oscillation. This is analogous to the bistable intensity found in an anharmonic cavity [37] driven at the fundamental. In this case the cavity is driven at its second harmonic.

3.3 Mean-field stability properties

We are interested in the stability of these steady states. To obtain this, we calculate the linearized equations for small perturbations, which reads

[TABLE]

The corresponding eigenvalues $\lambda_{\pm}$ are

[TABLE]

We obtain stability if all eigenvalues are negative, so that there are two conditions for stability:

[TABLE]

and:

[TABLE]

We now analyze the different types of solution:

Type 1 solutions

For $\beta^{(1)}=0$ the stability condition (a) implies that:

[TABLE]

which is always true, and also from condition (b) that the driving field is below threshold, i.e,

[TABLE]

Thus, the vacuum solution is stable below threshold and unstable above threshold.

Type 2 solutions

The $\beta^{(2)}$ solutions can occur for either sign of $\Re\left(c\right)$.

• •

Firstly, consider the case of $\Re\left(c\right)>0$. For the above threshold case, condition (a) implies that

[TABLE]

which is always true in the region of $\left|c\right|<1$ where the solution is valid. Condition (b) takes the form of

[TABLE]

which is equivalent to $\left|c\right|<1$. Thus the solution is stable when $\Re\left(c\right)>0$, $\left|c\right|<1$.

• •

Next, consider the case of $\Re\left(c\right)<0$. We find that both of the stability conditions are always true, provided the solution exists, which is for $|\Im(c)|\leq 1$ .

Type 3 solutions

The $\beta^{(3)}$ solutions only occur for $\Re\left(c\right)<0$. Condition (a) implies that, for stability,

[TABLE]

Condition (b) implies that

[TABLE]

which is only satisfied if $|c|<1$. Since the solution is valid only if $|c|>1$, the type 3 solution is never stable.

3.4 Phase-diagram and discrete time symmetry breaking

We have found three types of stationary solutions. Their behavior changes depending on the complex parameter $c$. This gives a definite phase-diagram, as shown in figure 1, with one phase for $|c|<1$, and two phases for $|c|>1$.

(I)

In the region of $\Re\left(c\right)>0$ we find only one stationary point, the $\beta^{(1)}$ stable vacuum solution, with $|c|>1$. This phase also extends to $\Re\left(c\right)<0$, in which case there is an additional constraint. We find a unique stable vacuum solution only if $1<|\Im(c)|<|c|$.

(II)

There are three stationary points when $|c|<1$, for either sign of $\Re\left(c\right)$. In this phase the $\beta^{(1)}$ vacuum solution exists but is unstable, while there are a pair of $\beta^{(2)}$ stable above threshold solutions. This is an example of bistability.

(III)

For $\Re\left(c\right)<0$, and $|\Im(c)|<1<|c|,$ there are five stationary points with $\beta^{(1)}$ stable, $\beta^{(2)}$ stable, and $\beta^{(3)}$ unstable. This is an example of tristability.

In summary, the vacuum steady state $\beta=0$ is stable only at a driving field below threshold, of $\mathcal{E}<\left|\gamma\right|$, or $\kappa\mathcal{E}_{2}<\gamma_{2}\left|\gamma\right|$, and the bistable states occur at large driving field, when $\kappa\mathcal{E}_{2}>\gamma_{2}\left|\gamma\right|$. There is also a tristable regime below threshold in which the vacuum state is stable, but there are stable solutions of finite amplitude as well. Otherwise, the real part of the eigenvalue will be positive, which means that the solutions become unstable in the corresponding parameter region.

Hence, in general there can be up to five stationary solutions. The nonvanishing solutions correspond to an output electromagnetic field of

[TABLE]

where $t_{1}$ is a time-origin that depends on the details of the input and output coupling.

The stable $\beta^{(2)}$ solutions exist in either bistable (phase II) or tristable (phase III) regimes. These solutions always come in pairs, whose two solutions differ by a time translation, since:

[TABLE]

Here $t_{2}=t_{1}+\pi/\omega$, which defines a second discrete time origin. Thus, the system can relax to a steady-state that corresponds to either time origin. This is called spontaneous discrete time symmetry breaking. However, this behavior is only possible if these are stable solutions.

These results imply that in the mean-field limit, the above-threshold steady-state solutions have a spontaneous discrete time-translation symmetry-breaking. However, so far we have ignored quantum fluctuations, which we turn to next.

4 Quantum steady-state and tunneling

In the mean field picture, once the subharmonic mode reaches a stable state, no further dynamics will be observed. This picture is not accurate once fluctuations are taken into account. In the presence of either classical (thermal noise) [38, 39] or quantum fluctuations [7], a steady state can reach the other steady state of the multistable solutions predicted by mean field theory. Beyond mean field theory, where damping and noise is included, the degenerate parametric oscillator obeys a master equation, where the steady state solutions and the switching rate between these states are obtained either by solving the corresponding Fokker-Planck equation [40, 5, 41] or solving for the zero eigenvalue and its eigenvector of the super-operator dictating the time evolution in a master equation, expanded over a basis representation such as the Fock state basis[42, 43, 41, 26, 27].

Hence, we now turn to the full quantum behavior of this system. This allows us to derive the full steady-state quantum statistics, and to demonstrate that spontaneous symmetry breaking has an exponentially long lifetime at large photon number. We note that, unlike the mean-field solutions, these full quantum solutions require the inclusion of off-diagonal coherent projectors, which increase the phase-space dimensionality compared to the classical, mean-field behaviour.

4.1 Steady-state quantum solutions

In the zero temperature case, the steady-state solution of the scaled Fokker-Planck equation (2.36) can be readily found by using the method of potential equations [44, 45, 46, 47]

[TABLE]

where $N$ is a normalization constant and $\Phi$ satisfies

[TABLE]

These differential equations can be obtained by inserting the form (4.1) into the Fokker-Planck equation (2.36) and setting $\partial P_{1}/\partial\tau=0$.

It is simplest to proceed by introducing a shifted $c$ parameter:

[TABLE]

where $\tilde{c}\rightarrow c$ in the classical limit of $n\rightarrow\infty$. Including quantum noise, the exact steady-state solution is defined by the potential:

[TABLE]

so that in the scaled coherent space:

[TABLE]

This gives the exact zero-temperature solution for the steady-state of the density matrix, provided it is accompanied by a choice of contours that leads to a solution that is bounded, which we treat in detail in the next subsection. While formally similar to previous solutions, we note that all the parameters here can have complex values, which is necessary when treating the physics of recent quantum circuit experiments.

In the tunneling calculation below, we assume that the resonant driving field is only added on mode $a_{2}$ so that $\mathcal{E}_{1}=0$ and $\Delta_{2}=0$, which is the situation in recent experiments [24]. In this case, the exact steady-state potential is:

[TABLE]

and the probability distribution is:

[TABLE]

The solution above is a time-symmetric mixed state, which includes an equal probability of observing either of the two possible output amplitudes. We now wish to calculate the rate at which the system can switch from one phase to the other, i.e., the rate at which time symmetry is restored once it is broken by an observation of one or the other of the two possible output amplitudes.

For simplicity in this discussion of tunneling, we will assume $\Re(\tilde{c})>0$ and in the phase (II) bistable region. We get potential solutions vanishing on boundaries when $\Re(\tilde{c})>0$, as shown in this work. In the case that $\Re(\tilde{c})<0$, the potential diverges on the boundaries. This means that the physics of tunneling may be very different, and requires a different type of manifold for its treatment, which we analyze elsewhere.

To this end, we discuss the physical interpretation of the parameter region $\Re(\tilde{c})>0$. We note that $\tilde{c}=\left(\gamma-g\right)/\left(gn\right)$, $\gamma=\gamma_{1}^{(1)}+i\Delta_{1}$ and $g=\gamma^{(2)}+i\chi$. As a result, we have

[TABLE]

Thus we will have $\gamma^{(2)}(\gamma_{1}^{(1)}-\gamma^{(2)})-\chi(\chi-\Delta_{1})>0$ in the region of $\Re(\tilde{c})>0$, which means that typically one has $\gamma_{1}^{(1)}>\gamma^{(2)}$, although there is also a nonlinear coupling $\chi$ and a detuning $\Delta_{1}$ which can change this relationship. We see that in the case that the detuning $\Delta_{1}$ is negligible, the linear damping rate is larger than the nonlinear coupling when $\Re\left(\tilde{c}\right)>0$.

4.2 Tunneling regime with $\Re(\tilde{c})>0$

Geometrically, we can regard the quantum dynamics as occurring via a distribution function defined on a two-dimensional manifold embedded in a four-dimensional complex space. In this paper, we focus on the tunneling-dominated regime of $\Re(\tilde{c})>0$, i.e. the region of large single-photon loss and small nonlinear coupling, so that the potential function vanishes at the boundaries of a probability domain defined by square boundaries at $\beta=\pm 1$ or $\beta^{+}=\pm 1$. The manifold must also includes the vacuum state at $\beta=\beta^{+}=0$, which is the starting point of any dynamical experiment.

We expect tunneling between minima of the potential as in earlier work [5, 41]. However, in this earlier work, the parameters were real and there was no anharmonicity. In the present case, the probability domain that includes these minima is no longer necessarily a plane with real values of $\beta,\beta^{+}$. We now analyze the locations of these minima in the four dimensional space of coherent amplitudes.

To find the stable points, we will solve two equations analogous to the mean-field stationarity conditions, but generalized to four dimensions:

[TABLE]

If $\beta$ and $\beta^{+}$ are nonzero real numbers, we see that $-2\beta^{+}$ is real, but in general $2\tilde{c}\beta/(1-\beta^{2})$ is complex, so the equations (4.9) can’t be satisfied. This means that there is generally at least one complex number in $\beta$ and $\beta^{+}$ for nontrivial solutions of the stationary points of the potential.

The potential function (4.6) could be complex because $\tilde{c}$, $\beta$ and $\beta^{+}$ are complex numbers. The stationary points obtained by equations (4.9) are divided into three types: the origin solution ($\beta=\beta^{+}=0$), the classical solutions ($\beta^{+}=\beta^{*}$) and the nonclassical solutions ($\beta^{+}=-\beta^{*}$). For all three types of solution, the potential functions (4.6) are always real. This means that we can study the stationary points in their neighborhoods to find whether they are local minima, maxima, or saddle points. In each case, we assume that required relations define a locally planar surface in a neighborhood of the solution, in order to define the derivatives.

4.3 Local stationary points

As mentioned above, we find three types of solutions on solving equations (4.9). It is common to use the Hessian matrix to determine whether the roots are local minima, maxima, or saddle points [48]. The Hessian matrix is defined by the second derivatives of the potential function:

[TABLE]

where

[TABLE]

If the Hessian matrix is positive definite at a stationary point $\vec{\beta}$, $\vec{\beta}$ is an isolated local minimum of the potential function $\Phi(\vec{\beta})$. For a $2\times 2$ matrix, positive definite is equivalent to a positive determinant $\left|M\right|>0$ and a positive trace $\Tr(M)>0$. Similarly, the Hessian matrix is negative definite at a stationary point when it is an isolated local maximum, which is equivalent to a positive determinant $\left|M\right|>0$ and a negative trace $\Tr(M)<0$. At a saddle point, the Hessian matrix has both positive and negative eigenvalues, which leads to a negative determinant $\left|M\right|<0$ [48, 49].

The first type of solution is at the origin, $\beta=\beta^{+}=0$. The potential is simply $\Phi^{(o)}\equiv\Phi(0,0)=0$, and in this case, the second derivatives are

[TABLE]

Therefore, the Hessian matrix determinant is obtained as

[TABLE]

If $|\tilde{c}|<1$, we have $\left|M^{(o)}\right|<0$, which means that the origin point is a saddle point. This is generally stable below threshold, and unstable above threshold, as expected from the previous mean-field analysis.

As we will see in the following, other minima as well as the quantum tunneling only occur when $|\tilde{c}|<1$, which means that this first type of solution plays the role of a saddle point in understanding quantum tunneling. In subharmonic generators, the bistable solutions only take place above threshold. In our calculations, bistable solutions are obtained as double minima in the manifold, occurring in the parameter region of $|\tilde{c}|<1$. It is also possible to have tristability, as explained in the mean-field section.

4.4 Classical stable points

The second type of stable point is $\beta^{+}=\beta^{\ast}=r\exp(-i\varphi)\neq 0$. These conditions would correspond to a coherent state projector, so we term them classical stable points, and they closely match the corresponding stable mean-field solutions for $n\gg 1$. For simplicity, we assume that we have a bistable or phase II situation, rather than the more complex tristable phase III situation.

In this case equations (4.9) can be transformed into

[TABLE]

and therefore,

[TABLE]

Because $r$ is real, we have

[TABLE]

Taking the positive sign to obtain the bistable region, the condition $r^{2}>0$ is equivalent to $|\tilde{c}|<1$. This means that $r^{2}=\left|\beta\right|^{2}<1$ as well. We note that for small values of $n$, the exact phase boundaries are modified due to the fact that $c\neq\tilde{c}.$ In fact it is better to regard this as a somewhat fuzzy criterion at small $n$, since quantum fluctuations tend to broaden these phase distinctions.

Labelling these stationary points as $\beta^{(c)}$, $\beta^{(c)+}$, we finally get

[TABLE]

Here we have set $\Pi(\tilde{c})=\sqrt{1-\Im\left(\tilde{c}\right)^{2}}$ as in the mean-field analysis, while $\beta^{(c)}$ and $\beta^{(c)+}$ are generally complex. These solutions correspond to the mean-field solutions above threshold obtained previously, except that now we use the quantum noise modified value of the coupling $\tilde{c},$ rather than its mean-field value. If $\tilde{c}$ is a real number, the results will reduce to $(\beta^{(c)},\beta^{(c)+})=(\pm\sqrt{1-\tilde{c}},\pm\sqrt{1-\tilde{c}})$, which corresponds to the first line of equation (4.7) in Ref. [41] with $n=\mu/g^{2}$, $\tilde{c}=\sigma/\mu$, $\alpha_{0}=\sqrt{n}\beta$, with no anharmonic term.

In this classical stationary point case, the second derivatives are

[TABLE]

Hence, the determinant of Hessian matrix is obtained,

[TABLE]

We expect that the classical stable points are only valid with the condition $|\tilde{c}|<1$. It is directly checked that $\left|M^{(c)}\right|>0$ and the real part $Re[\Phi_{11}^{(c)}]>0$. We note that $\Phi_{11}^{(c)}=\Phi_{22}^{(c)\ast}$, which leads to $\Tr[M^{(c)}]>0$. Thus, these stable points are minima in the quantum potential, that is, they are local attractors. As the bistable states only occur above threshold, the condition $|\tilde{c}|<1$ corresponds to the threshold value in degenerate parametric oscillators.

4.5 Complex phase-space manifold

As introduced at the beginning of the section 4.2 and in earlier work [5, 41], in order to study quantum tunneling analytically, we need to find a two-dimensional manifold embedded in a four-dimensional complex space. The potential function defined on this manifold vanishes at the boundaries $\beta=\pm 1$ or $\beta^{+}=\pm 1$. The manifold must also include the vacuum state at $\beta=\beta^{+}=0$, as well as the classical stable points (4.17), so that quantum tunneling can take place. In order to define them as local saddle points or minima, we assume that on the manifold there is a locally planar surface in a neighborhood of the solutions.

Given these considerations, we define a curved surface through the stable points, $\beta^{(c)}=e^{i\varphi}\left|\beta^{(c)}\right|$, parameterized as:

[TABLE]

For large $x,y$ we want to include the real boundaries such that $\beta=\pm 1$, $\beta^{+}=\pm 1$, are on the manifold. Therefore, we can modify this as:

[TABLE]

In the limit of $p\rightarrow 0$ this gives the correct behavior of the required manifold, as a tilted plane which is cutoff at the edges to give the square manifold with vanishing boundaries.

This manifold is plotted in the figure 2, for $p=0.01$. One could also modify the cutoff function to give different behavior in the classical and quantum directions, but we only wish to consider the simplest case here.

On this manifold (4.21), we show the potential $\Phi(\vec{\beta})$ (4.6) in figure 3. We can obtain the local minima and the saddle point from the real part of the potential in figure 3 (a) and (b). In the neighborhood of these points, figures 3 (c) and (d) show that $Im(\Phi)=0$. Thus it is valid to define local minima and a saddle point just as with real potentials.

We also show that quantum tunneling will take place between these classical minima through the saddle point in the A, where we generalized the potential-barrier approximation to complex cases and derived the tunneling time for a complex Fokker-Planck equation.

4.6 Quantum stable points

There is another possible stationary value of the potential, which is for$\beta=-\beta^{+\ast}$. These would correspond to a superposition of distinct coherent states, which is a uniquely quantum effect, so we term them quantum stable points. By labelling these points as $\beta^{(q)}$ and $\beta^{(q)+}$, we finally obtain

[TABLE]

with $|\tilde{c}|<1$, i.e. above threshold. If $c$ is a real number, the results reduce to $(\beta^{(q)},\beta^{(q)+})=(\pm\sqrt{1+\tilde{c}},\mp\sqrt{1+\tilde{c}})$, which corresponds to the second line of equation (4.7) in Ref. [41], where there is no anharmonic term.

In this case, the second derivatives are

[TABLE]

Hence, the Hessian determinant is

[TABLE]

Considering that $|\tilde{c}|<1$, it is easily checked that $\left|M^{(q)}\right|>0$ and $Re[\Phi_{11}^{(q)}]>0$. Because $\Phi_{22}^{(q)}$ is the complex conjugate of $\Phi_{11}^{(q)}$, the trace of Hessian matrix is then positive: $\Tr[M^{(q)}]>0$. Thus, these stable points are also expected to be minima of the potential.

We note that $|\beta^{(q)}|>1$ since $|\tilde{c}|<1$, so these are farther from the origin than the classical minima, and in fact outside the boundaries of the stable manifold considered here. Therefore, tunneling occurs between the stable classical minima through the saddle point at the origin. This type of quantum stable point becomes important for extremely strong coupling, and will be treated in detail elsewhere.

4.7 Tunneling rate

In order to calculate the tunneling rate between the classical minima through the saddle point at the origin, we use a transformation to define variables $(u,v)=F(\beta,\beta^{+})$ where the classical minimal points are placed on the axis of $u$. With these new variables, the diffusion coefficient is a constant. This simplifies the calculation of the tunneling time via the potential-barrier approximation [50, 51, 5, 41], which is generalized for complex cases in A.

Here we consider the transformation introduced in A, and combine equations (1.2) and (1.5) to give:

[TABLE]

The inverse transformation of equation (4.25) is

[TABLE]

with the notation that

[TABLE]

As explained in the A, $\phi=\psi-\theta/2$ with $\sin^{-1}\beta^{(c)}=re^{i\psi}$. Here $e^{i\theta}=g/|g|=n/\epsilon$ has been introduced in the section 2.6. In the following, we will express the classical stable points $\beta^{(c)}$ with the newly introduced variables $u$ and $v$. Then we can study the tunneling rate for our system by applying the analytic tunneling time result (1.27) obtained in the A.

Considering the manifold (4.21) where our potential is defined on, the variables $(u,v)$ can be parameterized directly,

[TABLE]

with $p\to 0$. It is straightforward to find that for the $\phi$ defined in the A, the classical minimal points will be placed on the axis of $u$, i.e. $v=0$. In this case, the Fokker-Planck equation (2.36) is transformed to

[TABLE]

Here we have introduce another shifted coupling constant which interpolates between the men-field coupling and the one used to analyse the $\beta$ manifold:

[TABLE]

which effectively defines the relevant region of the phase diagram. Thus the diffusion coefficients are all constant and equal, which means that we can use the potential-barrier approximation [50, 51, 5, 41], generalized in A, to obtain the tunneling rate.

We can obtain the Jacobean in the form of

[TABLE]

Using the relations $P^{\prime}_{ss}(u,v)=JP_{ss}(\beta,\beta^{+})$ and $P^{\prime}_{ss}(u,v)=N^{\prime}\exp(-\Phi(u,v))$, the potential is given by

[TABLE]

The minimal points of the potential in the $u$ and $v$ variables are found where the gradient of the potential is zero. Thus we find that the classical minimal points are located at $(u^{(c)},v^{(c)})=(\pm 2r\cos 2\phi,0)$, where

[TABLE]

with $\phi=\psi-\theta/2$ and hence

[TABLE]

From the manifold (4.21) and the transformation (4.25), we will find that the line of $v=0$ with real $u$ is on this manifold, where we will have $\Upsilon_{+}=\Upsilon_{-}^{\ast}$. Hence the potential (4.33) is proved to be real. In the meanwhile, the classical minimal points and the saddle point at the origin are all on this line. Considering

[TABLE]

the current $\vec{J}$ can be obtained easily via the Fokker-Planck equation (4.30). It is directly checked that the current through this line $J_{u}$ is real. This shows that quantum tunneling mostly occurs through this line. Further analysis is given in the A.

The second derivatives on this line are always real as well, given that

[TABLE]

Therefore, the potentials of the saddle points and the classical minimal points are

[TABLE]

The related second derivatives are therefore

[TABLE]

The tunneling time for a symmetric bistable potential in two dimensions is calculated using an extension of the Kramers method developed by Landauer and Swanson [50, 51], which is called the potential-barrier approximation [5, 41] and generalized for complex cases in A.

In the A, we have obtained the analytic formation of the tunneling time as shown in equation (1.27). Since the potential and the related second derivatives have been calculated in equations (4.38) and (4.39), the analytic result of the tunneling time is:

[TABLE]

which can be expressed in terms of our parameters as:

[TABLE]

Here we have used the relation $(1-B^{2})(1-B^{\ast 2})=|\bar{c}|^{2}$, which can be checked via equations (4.34) and (4.35). If we set $\theta=\phi=0$ and $\bar{c}=\bar{c}^{\ast}$, which means all the parameters are real, the tunneling time (4.41) can be simplified to give

[TABLE]

This simplified form agrees with equation (4.22) of Ref. [41], with parameters are defined as $n=\mu/g^{2}$, $\bar{c}=\bar{\sigma}/\mu$ and $g=g^{2}\gamma_{1}$. Thus, in the situation where all the parameters are real numbers, without anharmonic terms or detunings, the tunneling time (4.41) reduces to previous results Ref. [5, 41].

5 Number-state calculations

The tunneling time can also be obtained by solving the master equation (2.38) numerically in the number-state basis. In this basis the master equation reduces to an infinite matrix equation. Nevertheless, as any physical system has a finite energy, a suitable energy cutoff will reduce the system to a finite matrix equation. While this method is only numerically feasible for small photon number, it allows us to check the accuracy of the approximate analytic calculation given above.

5.1 Number-state basis expansions

We first expand the density operator $\rho$ in terms of its number-state matrix elements $\rho_{kl}$, which are defined by

[TABLE]

Thus, the time evolution is given by

[TABLE]

Then the master equation (2.38) under number-state basis becomes

[TABLE]

Here we have used the Einstein summation convention on identical indices. And $T_{ij}^{kl}$ is a four-dimensional transition matrix describing the rate of transition from the state $\rho_{kl}$ to the state $\rho_{ij}$, which is in the form of

[TABLE]

Here we have used the assumption $\Delta_{2}=0$ so that $\mathcal{E}$ is real, and

[TABLE]

The behavior of the system can be characterized in terms of the eigenvalues and eigenvectors of the transition matrix $T_{ij}^{kl}$. For instance, the eigenvector corresponding to the zero eigenvalue is exactly the steady state of the system. The first negative eigenvalue is related to the quantum tunneling rate [41].

In order to make the matrix finite, we set a photon number cutoff $N$ so that $0\leq i,j,k,l\leq N$. This approximation is valid if the high-photon-number states that are ignored play no significant role in determining the evolution of the system. The four-dimensional matrix $T_{ij}^{kl}$ can be reduced to a two-dimensional one $T_{\bar{\alpha}}^{\bar{\beta}}$ with this truncation that

[TABLE]

with

[TABLE]

where

[TABLE]

Here $\delta_{\bar{\alpha}}^{\bar{\beta}}$ is a Kronecker delta, and $\bar{\alpha}$, $\bar{\beta}$ are in the range of $[1,(N+1)^{2}]$. Note that the transition matrix is not Hermitian because of the single- and two-photon decay process.

We label the $k$-th eigenvalue by $\epsilon_{k}$ and its corresponding eigenvector by $\rho_{\bar{\alpha}}^{(k)}$ so that

[TABLE]

Here the coefficients $A_{k}$ define the initial state. We order the indices $k$ by the size of the real part of the eigenvalues, so that $\Re(\epsilon_{k})\geq\Re(\epsilon_{k+1})$. $\epsilon_{0}$ is the stable eigenvalue ($\epsilon_{0}=0$) and $\rho_{\bar{\alpha}}^{(0)}$ is the stable state. $\epsilon_{1}$ is the tunneling eigenvalue so the tunneling time is obtained [41],

[TABLE]

5.2 Tunneling time calculations and comparisons

We will compare the tunneling times obtained from using the P-representation (4.41) with those obtained using a number-state expansion (5.10). By changing the parameters $\gamma_{1}^{(1)}$, $\gamma^{(2)}$, $\mathcal{E}$ and $\Delta_{1}$, the results for the tunneling time are shown in the figure 4. Here we have noted that parameter $\tilde{c}$ consists of $\gamma_{1}^{(1)}$, $\gamma^{(2)}$, $\chi$, $\Delta_{1}$ and $\mathcal{E}$ as shown in equation (4.8), and $n$ consists of $\gamma^{(2)}$, $\chi$ and $\mathcal{E}$ shown in equation (2.35). The parameters has been chosen to satisfy $\Re(\tilde{c})>0$, with large single-photon loss and small nonlinear couplings, and $|\tilde{c}|<1$, above threshold.

Since there is a potential barrier between the minima, quantum tunneling is a slow effect. Hence, we use $\ln(\gamma_{1}^{(1)}T)$ to compare the results in the limit of large tunneling time $T$. It is expected that the analytic potential-barrier approximation will be most reliable for large tunneling times. We show below that, as expected, the analytic results agree exceptionally well with numerical results for large tunneling times, which is the limit of most interest for understanding spontaneously broken symmetry.

Figures 4 (a) and (b) show that damping speeds up quantum tunneling. Figure 4 (c) shows that an increased driving will increase the tunneling time, since the parameter $\mathcal{E}$ is proportional to the driving $\mathcal{E}_{2}$. The change in tunneling time with the detuning $\Delta_{1}$ is shown in figure 4 (d), where the largest tunneling time occurs at $\Delta_{1}=-0.2$kHz because of the nonzero nonlinearity $\chi$. With large detuning $\Delta_{1}$, the tunneling time will be reduced as shown in the figure 4 (d). It is also clear in the figure 4 (c) that the potential-barrier approximation fails when the quantum tunneling becomes too fast, as expected for this approach.

The behaviors of tunneling time changing with driving can be understood from the analytic equations. At large driving $\mathcal{E}$, the parameter $n$ increases while $n\bar{c}$ remains unchanged. From the formation of the classical stable points (4.17), it is directly checked that the stable points move away from each other in this case. In addition, the potential barrier becomes larger according to equations (4.38). Therefore, the quantum tunneling is suppressed when we increase the driving $\mathcal{E}$ considering the analytic result (4.41) provides $T\sim\exp(\Phi^{(o)}-\Phi^{(c)})$.

A similar analysis can be applied on the results of figure 4 (a) and (b) as well. When the dampings $\gamma_{1}^{(1)}$ and $\gamma^{(2)}$ become larger, we find that the positions of the stable points (4.17) only change slightly, but the potential barrier (4.38) will greatly decrease. Thus, increasing the dampings $\gamma_{1}^{(1)}$ and $\gamma^{(2)}$ speed up quantum tunneling.

The dependence of the tunneling rate changing on the nonlinearities $\chi$ and $\kappa$ are shown in figure 5. From the formations (4.8) and (4.31), $\tilde{c}$ and $\bar{c}$ will be real if $\chi=\gamma^{(2)}\Delta_{1}/\gamma_{1}^{(1)}$ considering there is always a nonzero damping $\gamma_{1}^{(1)}$ in realistic systems. In the case of $\Delta_{1}=0$, we will have $\chi=0$ when $\bar{c}$ is real. This is the case shown in figure 5 (a), where $\chi=0$ corresponds to the largest tunneling time. This shows that the nonlinear coupling $\chi$ will decrease the tunneling time. This result is consistent with the previous ones in similar by simpler systems [26, 27].

This behavior can also be obtained from the analytic result for the tunneling time (4.41). Considering the definitions $\bar{c}=(2\gamma-g)/(2gn)$, $n=\left|\mathcal{E}/g\right|$ and $g=\gamma^{(2)}+i\chi$, it is directly checked that $\bar{c}$ and $n$ will decrease if the nonlinearity $\chi$ becomes larger. Then from the results of the potential for the saddle point and the classical stable points (4.38), we see that the difference of the potentials becomes smaller. Hence, the nonlinearity $\chi$ will reduce the potential barrier height and thus speed up the tunneling.

The effects of the nonlinearity $\kappa$ are more complicated than $\chi$ as shown in figure 5 (b), which were not treated in previous study on single-mode nonlinear resonators [26, 27]. Here we show that $\kappa$ will increase the tunneling time when it is small, and decrease the tunneling time once it is large enough. This can be understood by considering that $\mathcal{E}\propto\kappa$ and $\gamma^{(2)}\propto\kappa^{2}$ if $\gamma_{1}^{(2)}$ is negligible. Thus when $\kappa$ is small, the effect of $\mathcal{E}$ dominates, which will increase the tunneling time. When $\kappa$ is large enough so that the effect of $\gamma^{(2)}$ becomes more important, the quantum tunneling will then be sped up.

We also note that the analytical results are not expected to agree very well with the numerical results with large nonlinearities $\chi$ and $\kappa$ if other parameters are fixed. This is because the dimensionless driving field is smaller, which also decreases the potential barrier height, and therefore reduces the validity of the analytic approximations used for tunneling calculations. Essentially, in this limit there is no real tunneling, and the broken time-translation symmetry is rapidly restored.

6 Conclusion

In this paper, we have studied general quantum subharmonic generation with additional detunings and anharmonicity, which has been experimental achieved [24] in superconducting microwave cavities. With driving, damping and nonlinearity considered, we obtained the steady-state solution of the Fokker-Planck equation using the adiabatic approximation and the zero-temperature limit, in order to understand pure quantum tunneling effects. Because of the nonlinearity, a complex parameter $\tilde{c}$ has been introduced. This means that the potential of the steady state is in general complex, which is different from the previously studied quantum optical subharmonic generation systems where the potentials were always real.

Quantum tunneling has been studied in this non-equilibrium system. This is related to quantum time symmetry breaking, as it defines the maximum time that a spontaneously broken time phase can exist before randomly switching to a different discrete time phase. By studying the manifold of the steady-state potential, we find that quantum tunneling will occur in the parameter region of $\Re(\tilde{c})>0$ and $|\tilde{c}|<1$, i.e the region of large single-photon loss, small nonlinear couplings and above threshold. The tunneling time has been obtained analytically using the potential-barrier approximation. In the expected domain of applicability of large tunneling time, the results agree with numerical calculations using a number-state basis.

These results show that the anharmonicity $\chi$ will enhances quantum tunneling rates compared to previous cases with no anharmonic term. This may have practical applications for escaping a local minimum in quantum neural networks [23, 22], where the global potential minimum is the desired computational solution.

This work is supported by the National Key R&D Program of China (Grants No. 2016YFA0301302 and No. 2018YFB1107200), the National Natural Science Foundation of China (Grants No. 11622428, No. 61475006, and No. 61675007) and the Graduate Academic Exchange Fund of Peking University. PDD and MDR thank the Australian Research Council and the hospitality of the Institute for Atomic and Molecular Physics (ITAMP) at Harvard University, supported by the NSF. This research has also been supported by the Australian Research Council Discovery Project Grants schemes under Grants DP180102470 and DP190101480.

Appendix A Tunneling rate for complex Fokker-Planck equations

Previous research on tunneling rates [51, 41, 5, 50] has treated real potentials. For our system, we obtain a complex potential barrier and a complex Fokker-Planck equation. Thus, we must generalize the potential-barrier approximation [51, 50] to treat such complex Fokker-Planck cases.

Without loss of generality, we will study the complex Fokker-Planck equation

[TABLE]

The notation $h.c.$ has the same meaning as we introduced in section 2.6. Our Fokker-Planck equation (2.36) can be written in this form with the transformation

[TABLE]

The manifold we are concerned with (4.21) is then transformed into:

[TABLE]

with the classical stationary points $\beta^{(c)}=\pm e^{i\varphi}\left|\beta^{(c)}\right|$ transformed into $\zeta^{(c)}=\pm re^{i\psi}$. This manifold gives the correct behavior in the limit of $p\rightarrow 0$, as introduced in the section 4.5. The steady-state solution can be expressed by a potential $\Phi$, which satisfies

[TABLE]

Thus the potential $\Phi(\zeta,\zeta^{+})$ is in general complex for complex variables $(\zeta,\zeta^{+})$. In the situation where $\zeta^{+}=\zeta^{*}$, it is directly checked that the potential $\Phi(\zeta,\zeta^{+})$ is real. The stationary points can be calculated by using the first derivatives, which are divided into three groups: the origin solution ($\zeta=\zeta^{+}=0$), the classical solutions ($\zeta^{+}=\zeta^{*}$) and the nonclassical solutions ($\zeta^{+}=-\zeta^{*}$). Here we are interested in the quantum tunneling between the classical stationary points through the origin.

Given the classical stationary points $\zeta^{(c)}=re^{i\psi}$, we will introduce the transformation

[TABLE]

where $\phi=\psi-\theta/2$. Then the inverse transformation takes the form

[TABLE]

In this case, the Fokker-Planck equation is transformed to

[TABLE]

The notation $h.c.$ indicates hermitian conjugate terms obtained by the replacement of $u\rightarrow u$, $v\to-v$ and the conjugation of all complex parameters. Considered the manifold (1.3), it is directly checked that $v$ is real on this manifold except on the boundaries, while $u$ is in general complex. But for the situation of $\zeta=\zeta^{+}$, we have $v=0$ with $u$ is real. All the points on the line of $v=0$ with real $u$ have real potential $\Phi(u,v)$ as well as the real second derivatives $\Phi_{uu}$ and $\Phi_{vv}$, which is proved in the section 4.7. The classical stationary points and the origin solution are all located on the line of $v=0$ with real $u$. We suppose that the classical stationary points are local minima and the origin is the saddle point, as we have in section 4.3. As discussed in the Ref. [51], the quantum tunneling will take place through the direction of $u$ because of the symmetry. In the following, we will reproduce the analysis of [51] so that the potential-barrier approximation can be generalized for the complex potential cases.

As introduced in Ref. [51], the current flow from a local minimum of the potential, located at negative $u$ and labelled as $C_{1}$, to the another minimum, located at positive $u$ and labelled as $C_{2}$, has the form

[TABLE]

Here $\rho$ is the density and $D$ is the diffusion coefficient. We have taken the zero-temperature limit as we did in the main text. The equilibrium case with Boltzmann distribution,

[TABLE]

leads to $j=0$. Thus we set

[TABLE]

in the non-equilibrium case, then the variation of $\eta$ indicates the extent of the deviation from equilibrium and the flow (1.8) becomes

[TABLE]

In the following, we assume that $D$ is constant in the neighborhood of the saddle point, which is exactly true in our situation where $D=\cos(2\phi)/n$. To obtain the magnitude of current (1.11), this assumption requires a relatively large value of $\nabla\eta$ near the saddle point, where $\exp(-\Phi)$ is small, and much smaller values of $\nabla\eta$ near the minima, where $\exp(\Phi)$ has larger values. Therefore the major departures from equilibrium take place only in the neighborhood of the saddle point.

For our case where the two potential minima are at the same value, Ref. [51] has shown that at the symmetry plane, $j$ is perpendicular to the symmetry plane, and $j$ has the same direction throughout the neighborhood of the saddle points. Here we will use the same assumption where $u$ has been proved to be this direction. Then equation (1.11) tells us that $\eta$ is only a function of $u$ in the neighborhood of the saddle point. Therefore, we can integrate equation (1.11).

[TABLE]

In the saddle-point neighborhood, the potential $\Phi$ depends quadratically on the spatial coordinates

[TABLE]

where the second derivatives $\xi_{u}^{(o)}=-\Phi_{uu}^{(o)}$, $\xi_{v}^{(o)}=\Phi_{vv}^{(o)}$ are both positive due to the manifold of the saddle-point neighborhood. Thus we have

[TABLE]

The continuity of current, $\nabla\cdot j=0$, requires that $j_{u}$ be independent of $u$. Considering that $\eta$ is only dependent on $u$ in the neighborhood of the saddle point, the factor $j_{u}\exp[\xi_{v}^{(o)}v^{2}/2]$ is thus a constant. The only remaining variable is $\exp[-\xi_{u}^{(o)}(u^{\prime})^{2}/2]$. The integrand is then large only at the saddle point $u=0$, and diminishes rapidly. Therefore, at a relatively short distance away from the saddle point, $\eta(u)$ approaches a constant limiting value: a positive value in the minimum $C_{1}$ and a negative value in the minimum $C_{2}$.

Next, we will evaluate the population difference $\Delta$, which is equal to twice the population of the classical minimum $C_{1}$,

[TABLE]

Here we have used an expansion appropriate to the minimum $C_{1}$

[TABLE]

where $\xi_{u}^{(c)}=\Phi_{uu}^{(c)}$, $\xi_{v}^{(c)}=\Phi_{vv}^{(c)}$ are both positive second derivatives due to the manifold of the neighborhood of the minimal point. Then we will get the difference

[TABLE]

In order to obtain the tunneling time, we need to evaluate the total current $J$ crossing the saddle point as well. Since $\eta$ is only dependent on $u$ in the neighborhood of the saddle point, the equation (1.11) is equivalent to

[TABLE]

In the symmetry plane containing the saddle point ($u=0$), we will then obtain

[TABLE]

Integrating over $v$ gives the total current

[TABLE]

The dimensionless tunneling time is then obtained as

[TABLE]

Now we need to evaluate $\eta(C_{1})/(\partial\eta/\partial u)_{u=0}$. From the equation (1.12), we will get

[TABLE]

where $j_{u}$ is the current density at the saddle point, which takes the form via equation (1.11)

[TABLE]

Then we have

[TABLE]

Thus, we get

[TABLE]

The tunneling time in dimensionless units then takes the form

[TABLE]

Noting that we have rescaled the time as $\tau=\mathcal{E}t$ in section 2.6 and $n=\mathcal{E}/|g|$, the dimensional tunneling time of our system is therefore

[TABLE]

which was used to obtain the analytic expression for the tunneling time in (4.41).