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

We discuss general properties of discrete time quantum symmetry breaking in superconducting quantum circuits. 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 above threshold. We also obtain analytic solutions for the tunneling time over which the time symmetry-breaking is lost. We find that additional anharmonic terms found in the superconducting case increase the tunneling rate. Our analytic results are confirmed by number state calculations.


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. 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 [3]. For quantum optical systems, an exact solution for the steady-state quantum density matrix [4] is known for a special case. This is 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.
As well as transient macroscopic superpositions and tunneling, subharmonic generators can generate space-time ordering [11]. These effects are related to time crystals [12], which are recent, similar phenomena. Physically, such devices are above threshold quantum squeezed state generators. Below threshold, they are used to reduce quantum noise in gravity-wave detectors [13,14], while networks of above threshold parametric devices are used for NP-hard optimization [15,16].
The most general case of single mode time symmetry breaking involves arbitrary detunings, nonlinear losses and anharmonic nonlinearities in the fundamental mode. This is characteristic, for example, of recent experimentally realized superconducting cavity experiments [17]. 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. The readout is pumped off-resonantly at a frequency ω p = 2ω − ω r , where ω and ω r are the storage and readout frequencies, respectively. 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 change the physics, and modify both the steady-state quantum behaviour and the tunneling rates compared to previous studies.
In this paper, the exact steady-state solution is derived, which is expressed in form of a 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, which equilibrates 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, which is shown to agree with those numerically obtained in the number-state basis. This opens the way to studying quantum tunneling in novel multi-mode non-equilibrium systems.
The remainder of this 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 make use of the generalized P-representation in section 3 to obtain a novel conditional Fokker-Planck equation. In section 4, the strongly damped pump mode is adiabatically eliminated to obtain a simpler one-mode Fokker-Planck equation. The tunneling case for large single-photon damping is studied in section 5, demonstrating agreement between the analytic results and numerical calculations using number states. Finally, section 6 gives an outlook and our conclusions, with proofs of the complex tunneling rate method given in an Appendix.

General model Hamiltonian
Quantum optical and quantum circuit physics are closely related. The principle 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.
We suppose that a k , a † k are the k-th mode annihilation and creation operators of modes at two different frequencies ω k in coupled resonant cavities. There is an overall quantum Hamiltonian given by is the Hamiltonian for the k-th mode, with driving, damping, and nonlinear terms. We define H (n) = k H (n) k as the sum over modes for the n − th type of interaction. The detailed structure of these terms is as follows: HereΓ k are the external reservoir coupling operators with reservoir Hamiltonians H R k , and E k are the envelope amplitudes of external coherent driving fields at angular frequency kω p for the k-th mode. The nonlinear parameters are κ for subharmonic generation and χ for anharmonicity. We suppose ω 2 2ω 1 , so the system can be externally driven simultaneously at fundamental and subharmonic frequencies.
In many typical experiments [17], the inputs are in the higher frequency mode, a 2 . Thus, the driving field E 1 = 0 for such experiments. We will keep this term in the general Hamiltonian and exact solutions, for completeness, but set it to zero in tunneling calculations below.

Master Equation
The system Hamiltonian H rev is defined as the reversible part of the Hamiltonian without reservoir couplings, so that: (2. 3) 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: and evolve the density matrix using the subtracted Hamiltonian, H SR . From now on we let a † k = a † k (0), and define: This transformation has the effect of changing the mode frequencies in the Hamiltonian to relative detunings, so that ω k → ∆ k = ω k − kω p . The resulting master equation [18,19] for the quantum density matrix ρ, on tracing over all the reservoirs in the Markovian approximation, is: Here γ (j) k are the linear and nonlinear amplitude relaxation rates for j-boson relaxation in the k-th mode, and n (j) k are the corresponding thermal occupations of the reservoirs. These are n (j) where jω 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 [20], but for simplicity we only include nonlinear losses for k = 1. Higher order decoherence is also possible [21], but not treated here.
We setÔ =â j k to describe general j-photon damping in the k-th mode, giving the the master equation super-operator for decoherence as (2.7) In this paper we will put n (1) k = 0, for simplicity, since these reservoirs have at least twice the frequency of the fundamental reservoirs, and will have lower thermal occupations, which we will neglect. The energy relaxation rate in each mode for single-particle decay is 2γ (1) k . 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. We include these fundamental mode 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.

Fokker-Planck and stochastic equations
We now introduce the generalized P-representation developed in Ref. [22]. This expands the quantum density matrix in terms of a complete operator basisΛ (α, α + ), and a P-representation P (α, α + ) so that: The operator basis uses off-diagonal coherent state projection operators, and has the form,Λ (α, α + ) = |α α + * | / α + * | α , where |α is a two-mode coherent state. This includes a normalizing factor so that the P-distribution integrates to unity. The integration measures dµ can be chosen as either a volume measure on a full eightdimensional complex space over the combined vector α = (α, α + ), 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 quantumdominated 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.

Fokker-Planck equation
All normally ordered quantum correlation functions are moments of the distribution, since: Using standard operator identities, the resulting Fokker-Planck equation has the form that: where we introduce the notation that D (n) k is a differential operator acting on the Pfunction. The derivative operators are: ∂ k ≡ ∂/∂α k , ∂ + k ≡ ∂/∂α + k , and the individual terms involved that correspond to each Hamiltonian coupling are:

Stochastic Equations
If all modes have strong damping, so that all boundary terms vanish in the Fokker-Planck equation, there is a corresponding stochastic equation which can be written in a combined vector form as: where the Gaussian noise term ζ has a vanishing mean, and the only nonzero correlations are: In this form it is convenient to introduce complex single particle decay rates γ k and two-particle decay rates g 1 , so that these parameters can be combined into complex rate terms: With these definitions, the combined deterministic or drift term becomes: (3.8) The corresponding combined stochastic coefficient is then: where the individual mode noise sub-matrices are: 1 is the thermal noise coefficient. Although these equations are useful when the damping rates of both modes are similar, they have no known analytic solutions. For this reason, we turn next to a method that allows us to derive a solution for the complex P-representation of the subharmonic mode, which allows us to analytically calculate the tunnelling rate.

Hybrid representation
We 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 define a conditional P 2 distribution for the second harmonic mode that depends on the amplitude of the first mode, so that:ρ Since this α 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: In the limit of γ 1 , the second harmonic mode is rapidly damped to a deterministic solution, There is a similar equation for α + 2 , which means that in the adiabatic limit . (3.14)

Adiabatic Elimination
We now introduce a reduced P-representation obtained by tracing over the secondharmonic mode, so thatρ 1 = Tr 2 (ρ). If we expand the quantum density matrix in terms of a single-mode operator basisΛ 1 ( α), and a single-mode P-representation P ( α), where α = α 1 , one then obtains: The operator basis uses coherent state projection operators, as before, but with a simpler form,Λ 1 ( α) ≡ |α α + * | / α + * | α . All normally ordered single-mode quantum correlation functions are moments of the distribution, since: With this approach, the elimination of the α 2 amplitude results in a single-mode Fokker-Planck equation, We have defined with its conjugate, ε * (α + ), as an input field that allows for depletion. Here γ = γ 1 , = κE 2 /γ 2 and g = g 1 + |κ| 2 /2γ 2 . The notation hc indicates hermitian conjugate terms obtained by the replacement of α → α + , and the conjugation of all complex parameters.
The zero temperature case where Γ 1 = 0 is of most interest here, as it leads to pure quantum tunneling effects. We note that provided the nonlinearity is not too strong, one may convert this into an equivalent stochastic equation, which has the form: The real, Gaussian noise terms dη i have non-vanishing correlations of dη i dη j = δ ij dt. However, this procedure is only valid when the resulting stochastic equation is confined to a bounded manifold, which may not be true in general. In particular, the presence of anharmonic couplings and detunings can lead to an unstable manifold, and more powerful analytic methods become necessary, as we explain below.

Equivalent Hamiltonian
As we see from the previous section, if the second harmonic mode is rapidly decaying, so γ 1 , it can be removed from the dynamical equations. After adiabatically eliminating the second harmonic mode [17], we obtain a new stochastic equation. This can also be interpreted as corresponding to the unitary dynamics for the fundamental quantum system with a ≡ a 1 , governed by the adiabatic Hamiltonian, In addition, there is an effective single-photon loss γ and a two-photon loss γ (2) , where Thus we have g = g 1 + |κ| 2 /2γ 2 = γ (2) + iχ. The master equation of the density matrix ρ 1 that is equivalent to the stochastic equation given above is then obtained as ∂ ∂t where we have used a zero-temperature limit. Similar systems have been studied before where quantum squeezing and bifurcation have been found [23][24][25][26], but the quantum tunneling is still lack of understanding. Here we will try to make it clear on the quantum tunneling. Equation (4.8) 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 (4.3).

Steady-state solution of the Fokker-Planck equation
In the zero temperature case, the steady-state solution of the Fokker-Planck equation (4.3) can be readily found by using the method of potential equations [27][28][29][30] where N is a normalization constant and Φ satisfies Introducing dimensionless parameters c = (γ − g)/g, e c = E 1 / √ g and λ c = /g, the potential solution is: Discrete time symmetry breaking in quantum circuits: exact solutions and tunneling9 Thus, the steady-state probability distribution is given that: If the driving field is only added on the mode a 2 with E 1 = 0 and hence e c = 0, which is the case we use in the tunneling calculation below, the steady-state probability distribution can be simplified as: This is the exact zero-temperature solution for the steady-state of the density matrix. It must be accompanied by a choice of contours that gives a solution that is bounded, which we treat in detail in the next section. While formally similar to previous solutions, we note that all the parameters here can have complex values, which is necessary when treating the parameter values found in recent quantum circuit experiments.
In the following, we make use of a scaled form of the solution obtained on defining To obtain a dimensionless equation, we note that characteristic time-scales are defined by the nonlinearity. Rescaling time such that τ = |g|t, and defining e iθ = g/|g|, gives 14) The exact steady-state solution is defined by the potential: so that in the scaled coherent space: In the tunneling calculation below, we assume that the resonant driving field is only added on the mode a 2 so that E 1 = 0 and ∆ 2 = 0, which is the same as the situation in the experiment [17]. Then the exact steady-state potential is simplified as (4.17) and the probability distribution 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(c) > 0, so that the potential function vanishes at the boundaries of a probability domain defined by square boundaries at β = ±1 or β + = ±1. The manifold must also includes the vacuum state at β = β + = 0, which is the starting point of any dynamical experiment.
We expect tunneling between minima of the potential as in earlier work [5,31]. However, in this earlier work, the parameters were real and there was no anharmonicity. In the present case, however, the probability domain that includes these minima is no longer necessarily a plane with real values of β, β + . We now analyze the locations of these minima in the four dimensional space of coherent amplitudes.
We note that c = (γ − g)/g, γ = γ 1 + i∆ 1 and g = γ (2) + iχ, which means that this is a relatively weak coupling or strong damping regime. As a result, we have In order to study the region of Re(c) > 0, we must have γ (2) (γ To find the stable points, we will solve the equations Here λ = |λ c | is a real parameter, and c is a complex parameter. If β and β + are nonzero real numbers, we see that −2λβ + is real but (2cβ)/(1 − β 2 ) is complex. Then the equations (5.2) can't be satisfied. This means that there is at least one complex number in β and β + for nontrivial solutions of the stationary points of the potential. The potential function (4.17) could be complex because c, β and β + are complex numbers. Fortunately, the stationary points obtained by equations (5.2) are divided into three types: the origin solution (β = β + = 0), the classical solutions (β + = β * ) and the nonclassical solutions (β + = −β * ). For all the three types of solution, the potential functions (4.17) 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.

Local stationary points
As mentioned above, we find three types of solutions on solving equations (5.2). It is common to use the Hessian matrix for determining whether the roots are local minima, maxima, or saddle points [32]. The Hessian matrix is defined by the second derivatives of the potential function If the Hessian matrix is positive definite at a stationary point β, β is an isolated local minimum of the potential function Φ( β). As for a 2 × 2 matrix, positive definite is equivalent to a positive determinant |M | > 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 |M | > 0 and a negative trace Tr(M ) < 0 in a 2 × 2 matrix. As for a saddle point, the Hessian matrix has both positive and negative eigenvalues, which leads to a negative determinant |M | < 0 in a 2 × 2 matrix [32,33].
Hence, the Hessian matrix determinant is obtained as If |c| < λ, we have M (o) < 0, which means that the origin point is a saddle point. As we will see in the following, other minima as well as the quantum tunneling only occur when |c| < λ, which means that the first type of solution plays the role of a saddle point in understanding the quantum tunneling phenomena.

Classical stable points
The second type is β + = β * = r exp(−iϕ) = 0. These conditions would correspond to a coherent state projector, so we term them classical stable points. In this case equations (5.2) can be transformed into and therefore, Because r is real, we have The condition r 2 > 0 is equivalent to |c| < λ. This means that r 2 = |β| 2 < 1 as well. Labelling these stationary points as β (c) , β (c)+ , we finally get Here we have set Λ(λ, c) = 4λ 2 + c 2 + c * 2 − 2|c| 2 for simplicity, while β (c) and β (c)+ are obviously complex numbers because c is a complex parameter. If c is a real number, the results will reduce to (β (c) , β (c)+ ) = (± (λ − c)/λ, ± (λ − c)/λ), which corresponds to the first line of equation (4.7) in Ref. [31] with λ = µ/g 2 , c = σ/g 2 , α 0 = √ λβ, with no anharmonic term. In this case, the second derivatives are , which leads to Tr[M (c) ] > 0. Thus, these stable points are expected to be minima in the potential, that is, they are local attractors. Suppose we define a plane through the stable points, β (c) = e iϕ β (c) . This is easily parameterized as: For large x, y we want to include the real boundaries such that β = ±1, β + = ±1, are on the manifold. Therefore, we can modify this as: β = x + ix tan (ϕ) cos p (xπ/2) cos p (yπ/2) , β + = y − iy tan (ϕ) cos p (xπ/2) cos p (yπ/2) . (5.14) In the limit of p → 0 this gives the correct behavior, as a tilted plane which is cutoff at the edges to give the square manifold with vanishing boundaries. This manifold is plotted in Fig. 1. 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 (5.14), we show the potential Φ( β) (4.17) in Fig. 2. We can figure out the local minima and the saddle point from the real part of the potential in Fig. 2 (a) and (b). In the neighborhood of them, Fig. 2 (c) and (d) shows that Im(Φ) = 0. Thus it is valid to define them as local minima and saddle point as in the real cases. We also show that quantum tunneling will take place between these classical minima through the saddle point in the Appendix Appendix A, where we generalized the potential-barrier approximation into complex cases and derived the tunneling time for the complex Fokker-Planck equation.

Quantum stable points
The last type is β = −β + * . 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 β (q) and β (q)+ , we finally get Discrete time symmetry breaking in quantum circuits: exact solutions and tunneling14 with λ > |c|. If c is a real number, the results will reduce to (β (q) , β (q)+ ) = (± (λ + c)/λ, ∓ (λ + c)/λ), which corresponds to the second line of equation (4.7) in Ref. [31], where there is no anharmonic term. In this case, the second derivatives are Hence, the Hessian determinant is Considering that λ > |c|, it is easily checked that M (q) > 0 and Re[Φ 22 is the complex conjugate of Φ (q) 11 , 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 also note that |β (q) | > 1 since λ > |c|, so these are farther from the origin than the classical minima, and in fact outside the boundaries of the stable manifold considered here. Therefore, the tunneling is expected to occur between the stable classical minima through the saddle point at the origin.

Tunneling rate
In order to calculate the tunneling rate between the classical minima through the saddle point at the origin, we need to find a transformation to define variables (u, v) = F (β, β + ) where the classical minimal points are placed on the axis of u. With these new variables, the diffusion coefficient should be a constant number in order to obtain the tunneling time via the potential-barrier approximation [5,31,34,35], which is generalized for complex cases in Appendix Appendix A.
Here we consider the transformation introduced in Appendix Appendix A u = e −iθ/2 e iφ sin −1 β + e iθ/2 e −iφ sin −1 β + , v = e −iθ/2 e −iφ sin −1 β − e iθ/2 e iφ sin −1 β + , (5.18) and the inverse transformation with the notation that where φ is a variable to guarantee that the classical minimal points will be placed on the axis of u. The manifold we have parameterized (5.14) are then transformed into the variables (u, v). It is straightforward to find that for a proper φ as in Appendix Appendix A, v is real on this manifold except on the boundaries, while u is in general complex. In this case, the Fokker-Planck equation (4.14) is transformed to with φ = ψ − θ/2 and hence From the manifold (5.14) and the transformation (5.18), we will find that the line of v = 0 with real u is on this manifold, where we will have Υ + = Υ * − . Hence the potential (5.23) 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 ∂P 1 ∂τ + ∇ · J = 0, (5.26) the current J can be obtained easily via the Fokker-Planck equation (5.21). It is directly checked that the current through this line J u is real. This could be a hint that the quantum tunneling could occur through this line. Further analysis can be found in the Appendix Appendix A. The second derivatives on this line are always real as well, given that (5.27) Therefore, we have the potential of the saddle points and the classical minimal points (5.28) The related second derivatives are The tunneling time for a symmetric bistable potential in two dimensions can be calculated using an extension of the Kramers method developed by Landauer and Swanson [34,35], which is called the potential-barrier approximation [5,31] and generalized for complex cases in Appendix Appendix A. Thus we have the tunneling time according to equation (A.27) Here we have used the relation (1 − B 2 )(1 − B * 2 ) = |c| 2 /λ 2 , which can be checked directly via equations (5.24) and (5.25). If we directly set θ = φ = 0 and c = c * , which means all the parameters are real, the tunneling time (5.30) can be simplified as which is consistent with the equation (4.22) of Ref. [31], where the parameters λ = µ/g 2 ,c =σ/g 2 and g = g 2 γ 1 . This means that in the situation where all the parameters are real numbers, without anharmonic terms, the tunneling time (5.30) reduces to the previous results Ref. [5,31].

Number-state expansion
The tunneling time is also able to be obtained by solving the master equation (4.8) 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. We first expand the density operator ρ in terms of its number-state matrix elements ρ mn , which are defined by ρ mn = m|ρ|n . Here we have used the Einstein summation convention on identical indices. And T mn ij is a four-dimensional transition matrix describing the rate of transition from the state ρ mn to the state ρ ij , which is in the form of The behavior of the system can be characterized in terms of the eigenvalues and eigenvectors of the transition matrix T mn ij . For instance, the eigenvector corresponding to the zero eigenvalue is exactly the steady state of the system. And the first negative eigenvalue is related to the quantum tunneling rate [31]. In order to make the matrix finite, we set a photon number cutoff N so that 0 ≤ i, j, m, n ≤ N . This type of approximation will be valid if the high-photonnumber states that are ignored play no significant role in determining the evolution of the system. What's more, the four-dimensional matrix T mn ij can be reduced to a two-dimensional one Tβ α with this truncation that d dt ρᾱ = Tβ α ρβ, Here δβ α is a Kronecker delta, andᾱ,β 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 λ k and its corresponding eigenvector by ρ (k) α so that (5.40) 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(λ k ) ≥ Re(λ k+1 ). λ 0 is the stable eigenvalue (λ 0 = 0) and ρ (0) α is the stable state. λ 1 is the tunneling eigenvalue so the tunneling time is obtained [31], We will compare the tunneling times obtained from using the Prepresentation (5.30) with those by number-state expansion (5.41). By changing the parameters γ (1) 1 , γ (2) , χ and ∆ 1 , the results of tunneling time are shown in Fig. 3.
The figures of ln(γ 1 T ) changing with the real and image parts of are shown in Fig. 4. Since the parameter is proportional to the driving E 2 , Fig. 4 shows that the driving will increase the tunneling time. It is expected that the potential-barrier approximation is reliable for the large tunneling time limit, which is demonstrated by both Fig. 3 and Fig. 4. When the quantum tunneling becomes fast, the potentialbarrier approximation fails. This is clear in Fig. 4 (a).

Conclusion
In this paper, we have studied general quantum subharmonic generation with additional detunings and anharmonicity, which has been experimental achieved [17] 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 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(c) > 0 and |c| < λ. 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 will greatly 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 [15,16], where the global potential minimum is the desired computational solution. Previous research on tunneling rates [5,31,34,35] 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 [34,35]  Thus the potential Φ(ζ, ζ + ) is in general complex for complex variables (ζ, ζ + ). In the situation where ζ + = ζ * , it is directly checked that the potential Φ(ζ, ζ + ) is real. The stationary points can be calculated by using the first derivatives, which are divided into three groups: the origin solution (ζ = ζ + = 0), the classical solutions (ζ + = ζ * ) and the nonclassical solutions (ζ + = −ζ * ). Here we are interested in the quantum tunneling between the classical stationary points through the origin. Given the classical stationary points ζ (c) = re iψ , we will introduce the transformation u = e −iθ/2 e iφ ζ + e iθ/2 e −iφ ζ + , v = e −iθ/2 e −iφ ζ − e iθ/2 e iφ ζ + , (A.5) where φ = ψ − θ/2. Then the inverse transformation takes the form In this case, the Fokker-Planck equation is transformed to The notation h.c. indicates hermitian conjugate terms obtained by the replacement of u → u, v → −v and the conjugation of all complex parameters. Considered the manifold (A.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 ζ = ζ + , we have v = 0 with u is real. All the points on the line of v = 0 with real u have real potential Φ(u, v) as well as the real second derivatives Φ uu and Φ vv , which is proved in the section 5.2. 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 5.1. As discussed in the Ref. [35], the quantum tunneling will take place through the direction of u because of the symmetry. In the following, we will reproduce the analysis of [35] so that the potential-barrier approximation can be generalized for the complex potential cases.
As introduced in Ref. [35], 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 j = −nD∇Φ − D∇n. (A.8) Here n is the density of the systems 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, 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φ). To obtain the magnitude of current (A.11), the assumption requires a relatively large value of ∇η near the saddle point, where exp(−Φ) is small, and much smaller values of ∇η near the minima, where exp(Φ) 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. [35] 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 (A.11) tells us that η is only a function of u in the neighborhood of the saddle point. Therefore, we can integrate equation. (A.11) In the saddle-point neighborhood, the potential Φ depends quadratically on the spatial coordinates where the second derivatives ξ vv are both positive due to the manifold of the saddle-point neighborhood. Thus we have 14) The continuity of current, ∇ · j = 0, requires that j u be independent of u. Considered that η is only dependent on u in the neighborhood of the saddle point, the factor j u exp[ξ u (u ) 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, η(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 ∆, which is equal to twice the population of the classical minimum C 1 , Here we have used an expansion appropriate to the minimum C 1 where ξ In the symmetry plane containing the saddle point (u = 0), we will then obtain Integrating over v gives the total current