Unraveling the switching dynamics in the double well of a Kerr-cat qubit

Qile Su [email protected] Department of Physics and Applied Physics, Yale University, New Haven, CT 06511, USA Rodrigo G. Cortiñas Department of Physics and Applied Physics, Yale University, New Haven, CT 06511, USA Jayameenakshi Venkatraman Department of Physics, University of California Santa Barbara, Santa Barbara, 93106, CA, U.S.A Shruti Puri Department of Physics and Applied Physics, Yale University, New Haven, CT 06511, USA

Abstract

A quantum particle trapped on one side of a symmetric double-well potential can spontaneously switch between wells by tunneling or by interacting with its fluctuating environment. It is generally believed that the switching rate decreases smoothly and exponentially with the size of the energy barrier. This view was challenged by a recent experiment performed on a driven superconducting nonlinear oscillator (often called the Kerr-cat qubit or the Kerr parametric oscillator), where the combination of nonlinearity and the drive stabilizes two symmetric states of oscillation by activating an electrical version of the double-well potential in the oscillator’s phase space. Specifically, it was found that as the drive amplitude and, correspondingly, the size of the energy barrier increases, the switching rate between the wells decreases in a step-like fashion that have been colloquially dubbed the “staircase”. In this work, we explain the staircase by deriving a semi-analytical formula for the switching rate as a function of the drive amplitude. The formula attributes each step in the staircase to a pair of energy levels connected via resonant tunneling. Our main insight behind the formula is that dissipation due to the environment dephases inter-well superpositions and therefore limits the tunneling rate via the quantum Zeno effect. Such competition between tunneling and dissipation generates a tunneling-dominated regime and a dissipation-dominated regime that are respectively observed in the flat part and the steep part of each step in the staircase. Our formula allows us to determine the critical drive amplitudes at which abrupt decreases in the switching rate occur. It also helps us identify the physical processes that move population from the ground state manifold to the excited state manifold dominating the switching rate. These processes are cascaded and direct thermal heating, rather than quantum heating as previously thought. Our theory provides a finer level of detail on switching processes in a double well potential in the quantum regime.

^†^†preprint: APS/123-QED

I Introduction

The quantum mechanical double-well potential is ubiquitous in physics. One example is the two-photon driven Kerr nonlinear oscillator milburn_quantum_1991; wielinga_quantum_1993; dykman_fluctuating_2012, which has two metastable states of oscillation. The bistability of this oscillator inspires its applications as quantum information processors cochrane_macroscopically_1999; goto_bifurcation-based_2016; goto_universal_2016; puri_engineering_2017; royer_qubit_2018; puri_stabilized_2019; puri_bias-preserving_2020; darmawan_practical_2021; bhandari_symmetrically_2024 and quantum sensors heugel_quantum_2019. In the first application, information is stored in the manifold formed by the two metastable states. In this context the oscillator is sometimes called the Kerr-cat qubit grimm_stabilization_2020; frattini_squeezed_2022; xu_engineering_2022; hajr_high-coherence_2024; ding_quantum_2024, or the Kerr parametric oscillator (KPO) qubit wang_quantum_2019; chono_two-qubit_2022; yamaguchi_spectroscopy_2024; iyama_observation_2024. In the second application, the oscillator couples to a weak perturbation that breaks the symmetry between the metastable states, and the unknown perturbation amplitude can be deduced from the critical oscillator frequency at which an abrupt transition from one metastable state to the other occurs.

A fundamental property of the double well is the spontaneous switching rate between the two metastable states. This property is relevant to the applications of the two-photon driven Kerr nonlinear oscillator because spontaneous switching (a) gradually degrades stored information when the oscillator operates as an information processor, and (b) introduces noise into the measured signal when the oscillator operates as a sensor. A recent experiment frattini_squeezed_2022 showed that in the quantum regime of a well-resolved energy levels inside the wells and in the absence of external drives that break the symmetry between the two wells, the spontaneous switching rate in the oscillator displays a step-like dependence as a function of the depth of the double well, a behaviour which has been referred to as the “staircase”. A detailed explanation of this behavior is lacking in previous theoretical works where the spontaneous switching rate has either been predicted to be constant gautier_combined_2022 or have a smooth exponential dependence marthaler_switching_2006; leggett_dynamics_1987; puri_bias-preserving_2020 on the depth of the double well.

In this work we bridge the gap between theory and experiments by finding a semi-analytical formula for the spontaneous switching rate of the two-photon driven Kerr-nonlinear oscillator. Each step in the staircase is explained by the transition from tunneling-dominated dynamics to dissipation-dominated dynamics in one of many two-level manifolds that are formed by quasi-degenerate states within the wells. Our formula predicts the depths of the double well, controlled by the normalized two-photon drive amplitude, at which steps in the staircase occur. The formula shows that the drive amplitude that is required for an overall exponential suppression of the spontaneous switching rate is such that the tunneling splitting of the first excited state manifold becomes small compared to the rate of dissipation. Independently from the semi-analytical formula, we also find numerically that the population of the manifold that contributes the most to switching is supplied by direct and cascaded thermal activation. Although in previous work the spontaneous switching has been attributed to quantum activation marthaler_switching_2006; vijay_invited_2009; ong_quantum_2013, where the system heats despite that the temperature of the environment approaches zero, we find this mechanism to not be dominant in low temperature limit.

The paper is organized as follows. In section II we provide background on the driven Kerr nonlinear oscillator by introducing its effective Hamiltonian, energy levels, and energy eigenstates. We introduce the Lindbladian that models the effect of the environment, and then define what we mean by the spontaneous switching rate. In section III we make approximations which are justified in the regime of interest, and then derive an expression for the spontaneous switching rate base on these equations and compare it to numerical simulations. In section IV we predict the critical values of the drive strength that controls the depth of the double well, at which the spontaneous switching rate changes abruptly. In section V we show that the equilibrium population inside a well is established by direct and cascaded thermal activation rather than quantum activation marthaler_switching_2006. In section VI we compare our work to previous works. Finally we conclude in section VII.

Table 1: Table of symbols.

Symbol	Definition / Result
Section II: Background
$K$	The size of the Kerr nonlinearity in a Kerr nonlinear oscillator, expressed in units of frequency.
$\epsilon_{2}$	The two-photon drive amplitude, expressed in units of frequency.
$\hat{H}$	The Hamiltonian of a two-photon driven Kerr nonlinear oscillator under the rotating-wave approximation and in a frame rotating at the oscillator frequency. It is written in units of frequency and with an overall minus sign factored out: $\hat{H}=-K\hat{a}^{\dagger 2}\hat{a}^{2}+\epsilon_{2}(\hat{a}^{\dagger 2}+\hat{a}^{2})$ [Eq. (1)], where $[\hat{a},\hat{a}^{\dagger}]=1$ .
$\ket{\pm\alpha}$	The two coherent states that spans the highest-energy manifold of $\hat{H}$ . The amplitude is $\alpha=\sqrt{\epsilon_{2}/K}$ [Eq. (2)].
$E_{0}$	The energy eigenvalue of the highest-energy manifold of $\hat{H}$ . Its value is $E_{0}=\epsilon_{2}^{2}/K$ [Eq. (2)].
$\ket{\mathcal{C}_{\alpha}^{\pm}}$	The even and odd parity Schrödinger cat states, defined as $\ket{\mathcal{C}_{\alpha}^{\pm}}=\mathcal{N}_{\alpha}^{\pm}(\ket{+\alpha}\pm\ket{-\alpha}$ ) [Eq. (3)], and $\mathcal{N}_{\alpha}^{\pm}=(2\pm 2e^{-2\|\alpha\|^{2}})^{-1/2}$ .
$\ket{\psi_{n}^{\pm}}$	Even and odd parity energy eigenstates of $\hat{H}$ [Eq. (6)]. The $n=0$ case reduces to $\ket{\psi_{0}^{\pm}}=\ket{\mathcal{C}_{\alpha}^{\pm}}$ .
$E_{n}^{\pm}$	Energy eigenvalues of $\hat{H}$ corresponding to the eigenstates $\ket{\psi_{n}^{\pm}}$ [Eq. (6)]. The $n=0$ case reduces to $E_{0}^{\pm}=E_{0}$ .
$U(\beta)$	The double-well potential in the phase space of the nonlinear oscillator obtained from the negative expectation value of $\hat{H}$ on coherent states $\ket{\beta}$ [Eq. (7)].
$\delta_{n}$	The frequency of the tunnel splitting between the even and odd parity states in manifold $n$ , defined as $\delta_{n}=E_{n}^{+}-E_{n}^{-}$ [Eq. (8)].
$E_{\text{gap},n}$	The energy gap between neighboring energy manifolds, written in units of frequency [Eq. (10)].
$\kappa$	The damping rate of the nonlinear oscillator due the environment.
$n_{\text{th}}$	The thermal occupation number of the environment.
$C$	A set of collapse operators [Eq. (15)]. It includes single photon loss and single photon gain: $C=\{\sqrt{\kappa(n_{\text{th}}+1)}\hat{a},\sqrt{\kappa n_{\text{th}}}\hat{a}^{\dagger}\}$ .
$\mathcal{D}[\hat{O}]$	The dissipator. It is defined as $\mathcal{D}[\hat{O}]\hat{\rho}=\hat{O}\hat{\rho}\hat{O}^{\dagger}-\frac{1}{2}\hat{O}^{\dagger}\hat{O}\hat{\rho}-\frac{1}{2}\hat{\rho}\hat{O}^{\dagger}\hat{O}$ . $\hat{O}$ is a collapse operator.
$\hat{\rho}$	The oscillators density operator.
$\mathcal{L}$	The Lindbladian in the master equation [Eq. (14)]. It is defined as $\mathcal{L}(\hat{\rho})=-i[\hat{H},\hat{\rho}]+\sum_{\hat{O}\in C}\mathcal{D}[\hat{O}]\hat{\rho}$ .
$\Gamma$	The spontaneous switching rate between the wells, defined as the rate of exponential dacay of the inter-well population difference after the intra-well quasi-equilibrium state has been achieved [Eq. (16)].
Section III: Solving for the spontaneous switching rate
$\ket{\psi_{n}^{p}}\bra{\psi_{m}^{q}}$	Terms in the density operator $\hat{\rho}$ when it is written in the basis of $\ket{\psi_{n}^{\pm}}$ . When $m=n$ and $p=q$ , it is referred to as population. When $m=n$ and $p\neq q$ , it is referred to as in-manifold coherence. When $m\neq n$ , it is referred to as cross-manifold coherence.
$\hat{I}_{n}$	Projector on to the manifold $n$ , defined as $\hat{I}_{n}=\ket{\psi_{n}^{+}}\bra{\psi_{n}^{+}}+\ket{\psi_{n}^{-}}\bra{\psi_{n}^{-}}$ [Eq. (18)].
$\mathcal{P}$	Projector that removes cross-manifold coherences from the density operator, defined as $\mathcal{P}(\hat{\rho})=\sum_{n}\hat{I}_{n}\hat{\rho}\hat{I}_{n}$ [Eq. (17)].
$\hat{\rho}_{\text{proj}}$	The density operator after the cross-manifold coherences have been removed. It is defined as $\hat{\rho}_{\text{proj}}=\mathcal{P}(\hat{\rho})$ . [Eq. (17)]
$\ket{\psi_{n}^{R/L}}$	Right- and left-well states defined as $\ket{\psi_{n}^{R/L}}=(\ket{\psi_{n}^{+}}\pm\ket{\psi_{n}^{-}})/\sqrt{2}$ .
$\hat{X}_{n}$	Pauli $x$ operator in the manifold $n$ . It is defined as $\hat{X}_{n}=\ket{\psi_{n}^{R}}\bra{\psi_{n}^{R}}-\ket{\psi_{n}^{L}}\bra{\psi_{n}^{L}}$ . [Eq. (21)]
$\hat{Z}_{n}$	Pauli $z$ operator in the manifold $n$ . It is defined as $\hat{Z}_{n}=\ket{\psi_{n}^{+}}\bra{\psi_{n}^{+}}-\ket{\psi_{n}^{-}}\bra{\psi_{n}^{-}}$ . [Eq. (22)]
$\hat{Y}_{n}$	Pauli $y$ operator in the manifold $n$ . It is defined as $\hat{Y}_{n}=i\hat{X}_{n}\hat{Z}_{n}$ . [Eq. (23)]
$\mathcal{L}_{\text{eff}}$	The effective Lindbladian after ignoring the cross-manifold coherences. It is defined as $\mathcal{L}_{\text{eff}}=\mathcal{P}\mathcal{L}\mathcal{P}$ . [Eq. (24)]
$\mathcal{L}_{\text{H}}$	The Hamiltonian part of the effective Lindbladian, defined as $\mathcal{L}_{\text{H}}(\hat{\rho})=-i\mathcal{P}[\hat{H},\mathcal{P}(\hat{\rho})]$ . [Eq. (25)]
$\mathcal{L}_{\text{D}}$	The dissipative part of the effective Lindbladian, defined as $\mathcal{L}_{\text{D}}(\hat{\rho})=\sum_{\hat{O}\in C}\mathcal{D}[\hat{O}]\mathcal{P}(\hat{\rho}$ ). [Eq. (26)]
$\hat{O}_{pq}$	The projected collapse operator, defined as $\hat{O}_{pq}=\sum_{i,j=R/L}\ket{\psi_{p}^{i}}\bra{\psi_{p}^{i}}\hat{O}\ket{\psi_{q}^{j}}\bra{\psi_{q}^{j}}$ [Eq. (27)]
$\hat{X}_{pq}$	Transition operator from manifold $q$ to manifold $p$ that simultaneously flips the Bloch vector around the $x$ -axis. It is defined as $\hat{X}_{pq}=\ket{\psi_{p}^{R}}\bra{\psi_{q}^{R}}-\ket{\psi_{p}^{L}}\bra{\psi_{q}^{L}}$ . [Eq. (31)]
$\hat{Y}_{pq}$	Transition operator from manifold $q$ to manifold $p$ that simultaneously flips the Bloch vector around the $y$ -axis. It is defined as $\hat{Y}_{pq}=i(\ket{\psi_{p}^{R}}\bra{\psi_{q}^{L}}-\ket{\psi_{p}^{L}}\bra{\psi_{q}^{R}})$ . [Eq. (32)]
$V_{pq}$	Inter-well transition rate, defined as $V_{pq}=\kappa(1+n_{\text{th}})\|\bra{\psi_{p}^{L}}\hat{a}\ket{\psi_{q}^{R}}\|^{2}+\kappa n_{\text{th}}\|\bra{\psi_{p}^{L}}\hat{a}^{\dagger}\ket{\psi_{q}^{R}}\|^{2}$ . [Eq. (45)]
$W_{pq}$	Intra-well transition rate, defined as $V_{pq}=\kappa(1+n_{\text{th}})\|\bra{\psi_{p}^{R}}\hat{a}\ket{\psi_{q}^{R}}\|^{2}+\kappa n_{\text{th}}\|\bra{\psi_{p}^{R}}\hat{a}^{\dagger}\ket{\psi_{q}^{R}}\|^{2}$ . [Eq. (47)]
$\hat{X}$	The total population difference operator between the wells, defined as $\hat{X}=\sum_{n}\hat{X}_{n}$ . [Eq. (53)]
$\Gamma_{n}$	The contribution to spontaneous switching rate $\Gamma$ from the manifold $n$ . It is defined by Eq. (56).
$\mu_{n}$	The rate at which the environment measures “which-well” information. It is defined as $\mu_{n}=2W_{nn}+\sum_{f\neq n}(W_{fn}+V_{fn})$ [Eq. (59)].
$\lambda_{n}$	The decay rate of population in the manifold $n$ due to transitions to other manifolds in the same well or to the other well (but excluding the effect of tunneling). $\lambda_{n}=2\sum_{f}V_{fn}+\sum_{f\neq n}(W_{fn}-V_{fn})$ [Eq. (69)].
$\delta_{n}^{2}/\mu_{n}$	The tunneling rate in manifold $n$ after taking into account of the quantum Zeno effect. It appears in the expression of the switching rate $\Gamma$ [Eq. (65)] and the effective equations of motion [Eq. (67)].
$R_{mn}$	A shorthand notation used to express the quasi-equilibrium distribution $\braket{\hat{X}_{n}}$ , $R_{mn}=(W_{mn}-V_{mn})/(\lambda_{m}+\delta_{m}^{2}/\mu_{m})$ [Eq. (74)].
$K_{n}$	The ratio between $\braket{\hat{X}_{n}}$ and $\braket{\hat{X}_{0}}$ , as in Eq. (76). It can be expressed analytically as in Eq. (77), $K_{n>0}=\sum_{\text{paths}}R_{ni_{s-1}}\cdots R_{i_{2}i_{1}}R_{i_{1}0}.$
$f_{n}$	Branching ratio in manifold $n$ , defined in Eq. (80), $f_{n}=(\delta_{n}^{2}/\mu_{n}+2\sum_{f}V_{fn})/(\delta_{n}^{2}/\mu_{n}+\lambda_{n})$ .
$J_{n}$	The inter-well difference of the total probability current entering manifold $n$ , defined in Eq. (81), $J_{n}=\sum_{i\neq n}(W_{ni}-V_{ni})(K_{i}/\sum_{m\geq 0}K_{m})$ .
$\Gamma_{n}(\delta_{k},W_{pq},V_{pq})$	The semi-analytical formula for the switching rate contributed by manifold $n$ in terms of the tunnel splittings and the intra-/inter-well transitions rates. $\Gamma_{0}(\delta_{k},W_{pq},V_{pq})=\left(\delta_{0}^{2}/\mu_{0}+2\sum_{f}V_{f0}\right)(K_{0}/\sum_{m\geq 0}K_{m})$ , $\Gamma_{n>0}(\delta_{k},W_{pq},V_{pq})=f_{n}J_{n}$ [Eq. (83)].
$\Gamma(\delta_{k},W_{pq},V_{pq})$	The semi-analytical formula for the switching rate in terms of the tunnel splittings and the intra-/inter-well transitions rates. $\Gamma(\delta_{k},W_{pq},V_{pq})=\sum_{n}\Gamma_{n}(\delta_{k},W_{pq},V_{pq})$ [Eq. (82)].
Section IV: Location of the steps in the staircase
$\delta_{n}^{2}/\mu_{n}=\lambda_{n}$	The condition Eq. (86) that determines the critical value of $\alpha$ at which there is a abrupt decrease in the switching rate. It is also the condition that determins the excited state manifold that dominates switching.
Section V: Application of the semi-analytical formula: activation mechanism within the well
$V_{pq}^{(\text{c.t.})},W_{pq}^{(\text{c.t.})}$	The inter-/intra-well transition rates [Eqs. (88) and (89)] that contains cascaded thermal heating but excludes quantum heating or any thermal heating that raises the manifold index by more than some integer $d$ .
$V_{pq}^{(\text{d.t.})},W_{pq}^{(\text{d.t.})}$	The inter-/intra-well transition rates [Eqs. (90) and (91)] that contains direct thermal heating from the ground state to an excited state but excludes quantum heating and any thermal heating from the excited states.
$\Gamma^{(\text{c.t.})}$	The switching rate predicted by the semi-analytical formula when only cascaded thermal heating is allowed. $\Gamma^{(\text{c.t.})}=\Gamma(\delta_{n},V_{pq}^{(\text{c.t.})},W_{pq}^{(\text{c.t.})})$ [Eq. (92)].
$\Gamma^{(\text{d.t.})}$	The switching rate predicted by the semi-analytical formula when only direct thermal heating is allowed. $\Gamma^{(\text{d.t.})}=\Gamma(\delta_{n},V_{pq}^{(\text{d.t.})},W_{pq}^{(\text{d.t.})})$ [Eq. (93)].

II Background

Refer to caption — Figure 1: The double-well shapes the spectrum of the Hamiltonian $\hat{H}$ defined in Eq. (1). (a) The double-well $U(\beta)=-\bra{\beta}\hat{H}\ket{\beta}$ , defined in Eq. (7), in the phase space of the oscillator in the rotating frame. The depth of the well (black arrows) is $\epsilon_{2}^{2}/K$ . Spontaneous switching (red arrows) reduces the inter-well population difference. The 2D projection below the double-well shows the energy contours. Black dots indicate the minimum points of the double-well, $\beta=\pm\alpha$ . (b) The energies $E_{n}^{+}$ of even parity excited states $\ket{\psi_{n}^{+}}$ (solid lines) and $E_{n}^{-}$ of odd parity excited states $\ket{\psi_{n}^{-}}$ (dashed lines) as functions of the two-photon drive amplitude $\epsilon_{2}$ , where $E_{n}^{\pm}$ are defined by Eq. (6). The shaded region indicates bound states for which $U(\pm\alpha)\leq-\bra{\psi_{n}^{\pm}}\hat{H}\ket{\psi_{n}^{\pm}}\lesssim U(0)$ . (c) The excited state manifold tunnel splitting $\delta_{n}$ , defined in Eq. (8), as a function of the two-photon drive amplitude $\epsilon_{2}$ . $\delta_{n}$ becomes exponentially small when the manifold $n$ drops into the double well. The states inside the well and the values of the tunnel splittings are important to determining the spontaneous switching rate.

In this section we introduce the system Hamiltonian and the Lindblad master equation that describes environmental decoherence. We start from the Hamiltonian of a two-photon driven Kerr nonlinear oscillator under the rotating-wave approximation and in a frame rotating at the oscillator frequency (using units where $\hbar=1$ ) puri_engineering_2017,

\displaystyle\hat{H}=-K\hat{a}^{\dagger 2}\hat{a}^{2}+\epsilon_{2}(\hat{a}^{\dagger 2}+\hat{a}^{2}),

(1)

where $K$ is the Kerr nonlinearity, $\epsilon_{2}$ is the amplitude of the two-photon drive, and $\hat{a}$ , $\hat{a}^{\dagger}$ are harmonic oscillator ladder operators satisfying $[\hat{a},\hat{a}^{\dagger}]=1$ . The negative sign in front of the Kerr term follows the convention of recent works where the nonlinearity originates from the Josephson junction. For $K>0$ , which we assume throughout this paper without loss of generality, the free oscillation frequency of the nonlinear oscillator decreases with amplitude. As a consequence of the sign choice, the energy eigenvalues of the Hamiltonian in the rotating frame is bounded above rather than below. The highest-energy manifold of the Hamiltonian is spanned by two coherent states $\ket{+\alpha}$ and $\ket{-\alpha}$ that have the degenerate energy eigenvalue $E_{0}$ puri_engineering_2017, where $\alpha$ and $E_{0}$ are defined as

\displaystyle\alpha=\sqrt{\epsilon_{2}/K},\quad E_{0}=\frac{\epsilon_{2}^{2}}{K}=K\alpha^{4}.

(2)

The two coherent states are not orthogonal but have an overlap of $\braket{-\alpha}{+\alpha}=e^{-2|\alpha|^{2}}$ which approaches zero in the limit $|\alpha|\to\infty$ . From the coherent states we can always construct an orthonormal basis for the ground state manifold consisting of the even- and odd-parity Schrödinger cat states cochrane_macroscopically_1999,

\displaystyle\ket{\mathcal{C}_{\alpha}^{\pm}}=\mathcal{N}_{\alpha}^{\pm}(\ket{+\alpha}\pm\ket{-\alpha}),

(3)

where the normalization constants are defined as $\mathcal{N}_{\alpha}^{\pm}=(2\pm 2e^{-2|\alpha|^{2}})^{-1/2}$ . We can use the Schrödinger cat states to define another orthonormal basis, $\ket{\psi_{0}^{R}}$ and $\ket{\psi_{0}^{L}}$ as follows,

$\displaystyle\Ket{\psi_{0}^{R/L}}$	$\displaystyle=\frac{\ket{\mathcal{C}_{\alpha}^{+}}\pm\ket{\mathcal{C}_{\alpha}^{-}}}{\sqrt{2}}$	(4)
	$\displaystyle\approx\left(1+\frac{3}{8}e^{-4\|\alpha\|^{2}}\right)\ket{\pm\alpha}-\frac{1}{2}e^{-2\|\alpha\|^{2}}\ket{\mp\alpha}$	(5)
	$\displaystyle\quad\quad\text{ for }\|\alpha\|^{2}\gg 1,$

which are asymptotic to the coherent states $\ket{\pm\alpha}$ in the limit of large $|\alpha|$ . We remark that the highest-energy manifold of $\hat{H}$ functions as the computational subspace when the system is used as a qubit grimm_stabilization_2020.

Because the Hamiltonian conserves parity, we can find a basis of simultaneous eigenstates of energy and parity, but unlike the highest energy eigenvalue, the rest of the energy eigenvalues are not degenerate. Letting $+$ and $-$ indicate even and odd parity and letting $n=0,1,2,...$ indicate decreasing energy in each parity sector, we denote the eigenstates and their energies by $\ket{\psi_{n}^{\pm}}$ and $E_{n}^{\pm}$ . As usual they satisfy the eigenvalue equation

\displaystyle\hat{H}\ket{\psi_{n}^{\pm}}=E_{n}^{\pm}\ket{\psi_{n}^{\pm}}.

(6)

The $n=0$ case reduces to $\ket{\psi_{0}^{\pm}}=\ket{\mathcal{C}_{\alpha}^{\pm}}$ and $E_{0}^{\pm}=E_{0}$ . Each $n$ represents a two-level manifold consisting of exactly one even-parity and one odd-parity eigenstate.

The energy spectrum of $\hat{H}$ and the wavefunctions of the eigenstates are molded by a symmetric double-well potential in phase space. The simplest way to reveal the double-well is by taking the negative expectation value of $\hat{H}$ on the coherent states $\ket{\beta}$ ,

\displaystyle U(\beta)=-\bra{\beta}\hat{H}\ket{\beta}

\displaystyle=K\beta^{*2}\beta^{2}-\epsilon_{2}(\beta^{*2}+\beta^{2}),

(7)

and varying $\beta$ in the complex plane (Fig. 1a). The reason we choose to define $U(\beta)$ with the added negative sign is to recover a conventional potential energy function whose graph opens upwards. Consistent with this choice, we refer to the $n=0$ manifold as the ground state manifold, and all other manifolds as excited state manifolds. We will also say a manifold $n$ lies above (below) another manifold $m$ if $n>m$ ( $n<m$ ), and we will from this point forward refrain from comparing energy eigenvalues.

$U(\beta)$ has two global minima located symmetrically at $U(\pm\alpha)=-\epsilon_{2}^{2}/K$ and a saddle point located at $U(0)=0$ ¹¹1The number of saddle points depend on the phase-space representation we choose for $\hat{H}$ . Here, the representation used in Eq. (7) is proportional to the $Q$ -function. Other phase space representations can also be used. The Wigner function, for example, has two saddle points for some values of $\epsilon_{2}$ venkatraman_driven_2024.. $\epsilon_{2}^{2}/K$ is the depth of the double-well and, equivalently, the height of the energy barrier, which we can control by changing the two-photon drive amplitude $\epsilon_{2}$ . The bound states of the double well are the eigenstates that satisfy $U(\pm\alpha)\leq-\bra{\psi_{n}^{\pm}}\hat{H}\ket{\psi_{n}^{\pm}}\lesssim U(0)$ , which include the two ground states and a subset of excited states. An estimate of the number $n_{\text{bound}}$ of quantized energy levels in each well is $n_{\text{bound}}\approx\epsilon_{2}/(\pi K)=\alpha^{2}/\pi$ given by the Bohr quantization condition frattini_squeezed_2022. As is commonly the case for double-well potentials, two types of energy gaps of distinct scales (Fig. 1b) characterize the energy spectrum for bound states puri_stabilized_2019; frattini_squeezed_2022. The first type consists of tunnel splittings $\delta_{n}$ with the following definition and scaling,

	$\displaystyle\delta_{n}$	$\displaystyle=E_{n}^{+}-E_{n}^{-}$		(8)
		$\displaystyle\sim f_{n}(\alpha^{2})e^{-2\alpha^{2}},\quad\text{ for }n\ll\alpha^{2}/\pi\text{ and }\alpha\gg 1,$		(9)

where $f_{n}$ is a polynomial function (Fig. 1c). The second type consists of gaps between consecutive manifolds with the following definition and scaling,

	$\displaystyle E_{\mathrm{gap},n}$	$\displaystyle=\|\bar{E}_{n+1}-\bar{E}_{n}\|$		(10)
		$\displaystyle\sim 4K\alpha^{2},\quad\text{ for }n\ll\alpha^{2}/\pi\text{ and }\alpha\gg 1,$		(11)

where $\bar{E}_{n}$ denotes the average energy in the $n$ th manifold.

The decoherence of the oscillator due to its weak coupling to an environment with a flat noise power spectrum can be modeled in a standard way by the Lindblad master equation, which includes the effect of single photon loss and single photon gain,

$\displaystyle\frac{\mathrm{d}\hat{\rho}}{\mathrm{d}t}$	$\displaystyle=\mathcal{L}(\hat{\rho})$	(12)
	$\displaystyle=-i[\hat{H},\hat{\rho}]+\kappa(1+n_{\mathrm{th}})\mathcal{D}[\hat{a}]\hat{\rho}+\kappa n_{\mathrm{th}}\mathcal{D}[\hat{a}^{\dagger}]\hat{\rho}$	(13)
	$\displaystyle=-i[\hat{H},\hat{\rho}]+\sum_{\hat{O}\in C}\mathcal{D}[\hat{O}]\hat{\rho}.$	(14)

$\hat{\rho}$ is the system density operator, $\mathcal{D}[\hat{O}]$ is a dissipator defined as $\mathcal{D}[\hat{O}]\hat{\rho}=\hat{O}\hat{\rho}\hat{O}^{\dagger}-\frac{1}{2}\hat{O}^{\dagger}\hat{O}\hat{\rho}-\frac{1}{2}\hat{\rho}\hat{O}^{\dagger}\hat{O}$ , the set of collapse operators is

\displaystyle C=\left\{\sqrt{\kappa(n_{\mathrm{th}}+1)}\hat{a},\sqrt{\kappa n_{\mathrm{th}}}\hat{a}^{\dagger}\right\},

(15)

$\kappa$ is the damping rate of the environment, and $n_{\text{th}}$ is the thermal occupation number of the environment. We work in the regime of $n_{\text{th}}\ll 1$ .

We define the spontaneous switching rate $\Gamma$ between the wells to be

\displaystyle\Gamma=-\frac{1}{\braket{\hat{X}}}\frac{\mathrm{d}}{\mathrm{d}t}\braket{\hat{X}},

(16)

where $\hat{X}$ is an observable that measures the population difference between the wells, to be defined later in Eq. (53) of section III.3, and the expectation value is evaluated at a time after the intra-well quasi-equilibrium state has been attained. Numerical solution of Eq. (14) qualitatively reproduces the experimental staircase, that is, as a function of the amplitude of the two-photon drive $\epsilon_{2}$ , which controls the depth of the double well, the spontaneous switching rate $\Gamma$ between the wells shows step-like decreases at specific values of $\epsilon_{2}$ (Fig. 2). See appendix A for details of the numerical solution johansson_qutip_2012; johansson_qutip_2013. It is the goal of the subsequent sections to explain these steps by finding a semi-analytical formula for the spontaneous switching rate $\Gamma$ .

III Solving for the spontaneous switching rate

We solve the Lindblad master equation for the spontaneous switching rate $\Gamma$ by making a series of approximations. These approximations allow us to obtain a semianalytical expression of $\Gamma$ in terms of the tunnel splittings $\delta_{n}$ , the single-photon loss rate $\kappa$ , the thermal occupation $n_{\text{th}}$ , and the matrix elements of $\hat{a}$ and $\hat{a}^{\dagger}$ .

III.1 Projection of the density operator onto two-level manifolds

We start with a perturbation theory argument showing that given the Lindblad master equation of Eq. (14) and the condition $\kappa,\kappa n_{\text{th}}\ll\min(E_{\text{gap},n})$ , coherent dynamics is largely confined within two-level manifolds. First we divide the terms of $\hat{\rho}$ in the energy-parity eigenbasis $\ket{\psi_{n}^{\pm}}$ into three categories: populations refer to terms of the form $\ket{\psi_{n}^{\pm}}\bra{\psi_{n}^{\pm}}$ ; in-manifold coherences are terms of the form $\ket{\psi_{n}^{\pm}}\bra{\psi_{n}^{\mp}}$ , and cross-manifold coherences are those of the form $\ket{\psi_{n}^{p}}\bra{\psi_{m}^{q}}$ for all $n\neq m$ and $p,q\in\{+,-\}$ . In the absence of dissipation ( $\kappa,\kappa n_{\text{th}}=0$ ), Hamiltonian evolution causes the cross-manifold coherences to advance in phase at the angular frequency of $|E_{n}^{p}-E_{m}^{q}|$ , the energy gap between manifolds $n$ and $m$ . In contrast, Hamiltonian evolution leaves populations invariant, and only causes the in-manifold coherences for manifolds under the energy barrier to advance slowly in phase, at the angular frequency of the tunnel splittings $\delta_{n}$ . Thus there is a frequency gap of order $\min(E_{\text{gap},n})$ from the populations and in-manifold coherences to cross-manifold coherences. Now we regard the dissipators $\mathcal{D}[\hat{O}]$ as a perturbation and perform time-dependent perturbation theory. The perturbation only creates a static coupling of order $\kappa$ and $\kappa n_{\text{th}}$ from the slow components to the fast components. Therefore, if the initial $\rho$ contains only populations and in-manifold coherences, then the build-up of cross-manifold coherences are off-resonantly suppressed on the order of $O[\max(\kappa,\kappa n_{\text{th}})/\min(E_{\text{gap},n})]$ . We remark that the confinement of coherent dynamics to within two-level manifolds is a direct extension of the fact that energy eigenstates become pointer states of decoherence under the combination of a strong non-degenerate system self-Hamiltonian and a weak interaction with the environment paz_quantum_1999.

Given our preceding perturbation theory argument, in the limit $\kappa,\kappa n_{\mathrm{th}}\ll E_{\mathrm{gap},n}$ the system density operator $\hat{\rho}$ is well approximated by an incoherent mixture across the manifolds within which coherent dynamics may occur. Equivalently speaking, the manifold index $n$ can be thought of as a classical random variable. This motivates the first approximation under which the density operator is projected onto two-level manifolds (Fig. 3a),

\displaystyle\hat{\rho}_{\text{proj}}

\displaystyle=\mathcal{P}(\hat{\rho})=\sum_{n}\hat{I}_{n}\hat{\rho}\hat{I}_{n},

(17)

where

\displaystyle\hat{I}_{n}=\ket{\psi_{n}^{+}}\bra{\psi_{n}^{+}}+\ket{\psi_{n}^{-}}\bra{\psi_{n}^{-}}

(18)

is the projector onto the $n$ th manifold. The projection preserves the populations and the in-manifold coherences while projecting out the cross-manifold coherences. We remark that in the limit of very small $n_{\text{th}}$ we must take into account the small correction to the energy eigenstates $\ket{\psi_{n}^{\pm}}$ due to the non-hermitian part of the effective Hamiltonian $\hat{H}_{\text{eff}}=\hat{H}-i\kappa\hat{a}^{\dagger}\hat{a}/2$ , but we defer this discussion to appendix B, where we numerical compute the effect of this projection on the switching rate. The projected density operator $\hat{\rho}_{\text{proj}}$ can be visualized as a collection of Bloch spheres (Fig. 3c). The north and south poles of each Bloch sphere represent the even- and odd-parity eigenstates in a particular manifold, and the equator represents a family of equal-weight superpositions between the parity eigenstates parameterized by their phase difference. Although $\ket{\psi_{n}^{+}}$ and $\ket{\psi_{n}^{-}}$ each has equal probability of occupying either well due to their being parity eigenstates, they can interfere with each other constructively or destructively. We define the right-well state $\ket{\psi_{n}^{R}}$ (left-well state $\ket{\psi_{n}^{L}}$ ) to be the particular equal-weight superposition of $\ket{\psi_{n}^{+}}$ and $\ket{\psi_{n}^{-}}$ that maximizes (minimizes) $\mathrm{Re}\braket{\hat{a}}$ (Fig. 3b). In particular this definition ensures $\bra{\psi_{n}^{R}}\hat{a}\ket{\psi_{n}^{R}}=-\bra{\psi_{n}^{L}}\hat{a}\ket{\psi_{n}^{L}}>0$ . Once the phase difference used to define the right- and left-well states are absorbed into $\ket{\psi_{n}^{+}}$ and $\ket{\psi_{n}^{-}}$ , $\ket{\psi_{n}^{R}}$ and $\ket{\psi_{n}^{L}}$ can be expressed as

	$\displaystyle\ket{\psi_{n}^{R}}$	$\displaystyle=\frac{1}{\sqrt{2}}\left(\ket{\psi_{n}^{+}}+\ket{\psi_{n}^{-}}\right),$		(19)
	$\displaystyle\ket{\psi_{n}^{L}}$	$\displaystyle=\frac{1}{\sqrt{2}}\left(\ket{\psi_{n}^{+}}-\ket{\psi_{n}^{-}}\right),$		(20)

hence fixing them to the $x$ -axis of the Bloch sphere. The components of the Bloch vector along the cardinal axes are given in the usual way by the expectation values of the following Pauli operators,

$\displaystyle\hat{X}_{n}$	$\displaystyle=\ket{\psi_{n}^{R}}\bra{\psi_{n}^{R}}-\ket{\psi_{n}^{L}}\bra{\psi_{n}^{L}}$	(21)
$\displaystyle\hat{Z}_{n}$	$\displaystyle=\ket{\psi_{n}^{+}}\bra{\psi_{n}^{+}}-\ket{\psi_{n}^{-}}\bra{\psi_{n}^{-}}$	(22)
$\displaystyle\hat{Y}_{n}$	$\displaystyle=i\hat{X}_{n}\hat{Z}_{n}$	(23)
	$\displaystyle=i(\ket{\psi_{n}^{R}}\bra{\psi_{n}^{L}}-\ket{\psi_{n}^{L}}\bra{\psi_{n}^{R}}).$

Our approximation of ignoring the cross-manifold coherence in Eq. (17) simplifies the Lindbladian dynamics because it reduces the degrees of freedom in the system. The effective Lindbladian restricted to the remaining degrees of freedom is

\displaystyle\mathcal{L}_{\text{eff}}

\displaystyle=\mathcal{P}\mathcal{L}\mathcal{P}=\mathcal{L}_{\text{H}}+\mathcal{L}_{\text{D}},

(24)

where at the second equality sign line we separated the effective Lindbladian into a Hamiltonian part $\mathcal{L}_{\text{H}}$ and a dissipative part $\mathcal{L}_{\text{D}}$ . After some algebra, we can show that the explicit expression for $\mathcal{L}_{\text{H}}$ is

\displaystyle\mathcal{L}_{\text{H}}

\displaystyle=-i\mathcal{P}[\hat{H},\mathcal{P}\ \cdot\ ]=-i\sum_{n}\frac{\delta_{n}}{2}[\hat{Z}_{n},\ \cdot\ ],

(25)

and explicit expression for $\mathcal{L}_{\text{D}}$ is

\displaystyle\mathcal{L}_{\text{D}}

\displaystyle=\sum_{\hat{O}\in C}\mathcal{P}\mathcal{D}[\hat{O}]\mathcal{P}=\sum_{\hat{O}\in C}\sum_{pq}\mathcal{D}[\hat{O}_{pq}],

(26)

where we define $\hat{O}_{pq}=\hat{I}_{p}\hat{O}\hat{I}_{q}$ , the projection of the collapse operator onto given initial and final manifolds $q$ and $p$ . As a reminder, the collapse operators are drawn from the set $C$ defined in Eq. (15) and includes single photon loss and gain. Consistent with the approximation of Eq. (17), the dissipators $\mathcal{D}[\hat{O}_{pq}]$ do not create superpositions between different manifolds.

Physically, the new dissipators $\mathcal{D}[\hat{O}_{pq}]$ cause random transitions across the Bloch spheres as well as random flips of the Bloch vector within each Bloch sphere. To illustrate this, we write the projected collapse operator $\hat{O}_{pq}$ in terms of its matrix elements

\displaystyle\hat{O}_{pq}

\displaystyle=\sum_{i,j=L/R}\ket{\psi_{p}^{i}}\bra{\psi_{p}^{i}}\hat{O}\ket{\psi_{q}^{j}}\bra{\psi_{q}^{j}}.

(27)

Since the operator $\hat{O}\in C$ is proportional to either $\hat{a}$ or $\hat{a}^{\dagger}$ , it is has odd parity symmetry,

\displaystyle e^{i\pi\hat{a}^{\dagger}\hat{a}}\hat{O}e^{-i\pi\hat{a}^{\dagger}\hat{a}}=-\hat{O}.

(28)

The states $\ket{\psi_{n}^{R}}$ and $\ket{\psi_{n}^{L}}$ transform into each other under the parity operation, so the expression for $\hat{O}_{pq}$ can be simplified into,

$\displaystyle\hat{O}_{pq}$	$\displaystyle=\bra{\psi_{p}^{R}}\hat{O}\ket{\psi_{q}^{R}}\left(\ket{\psi_{p}^{R}}\bra{\psi_{q}^{R}}-\ket{\psi_{p}^{L}}\bra{\psi_{q}^{L}}\right)$
	$\displaystyle\quad+\bra{\psi_{p}^{L}}\hat{O}\ket{\psi_{q}^{R}}\left(\ket{\psi_{p}^{L}}\bra{\psi_{q}^{R}}-\ket{\psi_{p}^{R}}\bra{\psi_{q}^{L}}\right)$	(29)
	$\displaystyle=\bra{\psi_{p}^{R}}\hat{O}\ket{\psi_{q}^{R}}\hat{X}_{pq}+i\bra{\psi_{p}^{L}}\hat{O}\ket{\psi_{q}^{R}}\hat{Y}_{pq},$	(30)

where we define new operators

	$\displaystyle\hat{X}_{pq}$	$\displaystyle=\ket{\psi_{p}^{R}}\bra{\psi_{q}^{R}}-\ket{\psi_{p}^{L}}\bra{\psi_{q}^{L}}$		(31)
	$\displaystyle\hat{Y}_{pq}$	$\displaystyle=i\left(\ket{\psi_{p}^{R}}\bra{\psi_{q}^{L}}-\ket{\psi_{p}^{L}}\bra{\psi_{q}^{R}}\right).$		(32)

The operators $\hat{X}_{pq}$ and $\hat{Y}_{pq}$ generate transitions across Bloch spheres when $p\neq q$ and flips the Bloch vector within the $p$ th Bloch sphere when $p=q$ . The notation for the two operators indicate the axis around which the Bloch vector is flipped when a transition happens.

III.2 Effective equations of motion of the density operator components

In this section we derive the effective equations of motion of the density operator components and, as we do so, develop an intuitive understanding of the terms involved. Under the approximation of Eq. (17), the density operator can be written as

\displaystyle\hat{\rho}_{\text{proj}}=\frac{1}{2}\sum_{\hat{S}\in\mathcal{S}}\braket{\hat{S}}\hat{S},

(33)

where $\mathcal{S}=\bigcup_{n}\{\hat{I}_{n},\hat{X}_{n},\hat{Y}_{n},\hat{Z}_{n}\}$ . Then the equation of motion of the expectation value of any $\hat{S}\in\mathcal{S}$ can be written as

	$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\braket{\hat{S}}$	$\displaystyle=\mathrm{Tr}[\hat{S}\mathcal{L}_{\text{eff}}(\hat{\rho}_{\text{proj}})]$		(34)
		$\displaystyle=\frac{1}{2}\sum_{\hat{S}^{\prime}\in\mathcal{S}}\mathrm{Tr}[\hat{S}\mathcal{L}_{\text{eff}}(\hat{S}^{\prime})]\braket{\hat{S}^{\prime}}.$		(35)

The following coefficients vanish,

	$\displaystyle\mathrm{Tr}[\hat{X}_{n}\mathcal{L}_{\text{eff}}(\hat{I}_{m})]=\mathrm{Tr}[\hat{X}_{n}\mathcal{L}_{\text{eff}}(\hat{Z}_{m})]=0,$		(36)
	$\displaystyle\mathrm{Tr}[\hat{Y}_{n}\mathcal{L}_{\text{eff}}(\hat{I}_{m})]=\mathrm{Tr}[\hat{Y}_{n}\mathcal{L}_{\text{eff}}(\hat{Z}_{m})]=0,$		(37)

indicating that the equations of motion for $\braket{\hat{X}_{n}}$ and $\braket{\hat{Y}_{n}}$ do not depend on $\braket{\hat{I}_{n}}$ or $\braket{\hat{Z}_{n}}$ . This is convenient because for the purpose of finding the spontaneous switching rate between the wells, we are only interested in any processes that change the $x$ component of the Bloch vector. Hence we only need the equations of motion for $\braket{\hat{X}_{n}}$ and $\braket{\hat{Y}_{n}}$ ,

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\braket{\hat{X}_{n}}$	$\displaystyle=\sum_{m}\left\{\frac{1}{2}\mathrm{Tr}[\hat{X}_{n}\mathcal{L}_{\text{eff}}(\hat{X}_{m})]\braket{\hat{X}_{m}}\right.$
	$\displaystyle\quad+\left.\frac{1}{2}\mathrm{Tr}[\hat{X}_{n}\mathcal{L}_{\text{eff}}(\hat{Y}_{m})]\braket{\hat{Y}_{m}}\right\},$	(38)
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\braket{\hat{Y}_{n}}$	$\displaystyle=\sum_{m}\left\{\frac{1}{2}\mathrm{Tr}[\hat{Y}_{n}\mathcal{L}_{\text{eff}}(\hat{Y}_{m})]\braket{\hat{Y}_{m}}\right.$
	$\displaystyle\quad+\left.\frac{1}{2}\mathrm{Tr}[\hat{Y}_{n}\mathcal{L}_{\text{eff}}(\hat{X}_{m})]\braket{\hat{X}_{m}}\right\}.$	(39)

To proceed, we must evaluate and interpret the action of the effective Lindbladian $\mathcal{L}_{\text{eff}}=\mathcal{L}_{\text{H}}+\mathcal{L}_{\text{D}}$ (defined in Eqs. (24) to (26)) on $\hat{X}_{n}$ and $\hat{Y}_{n}$ . For the Hamiltonian part $\mathcal{L}_{\text{H}}$ , we have

	$\displaystyle\mathcal{L}_{\text{H}}(\hat{X}_{n})$	$\displaystyle=+\delta_{n}\hat{Y}_{n},$		(40)
	$\displaystyle\mathcal{L}_{\text{H}}(\hat{Y}_{n})$	$\displaystyle=-\delta_{n}\hat{X}_{n}.$		(41)

This shows that the Hamiltonian part causes inter-well coherent tunneling in the $n$ th manifold. Now we evaluate and interpret the action of the dissipative part $\mathcal{L}_{\text{D}}$ . Because $\mathcal{L}_{\text{D}}$ can be expressed as a sum over a collection of dissipators following Eq. (26), it suffices to look at the action of a single dissipator $\mathcal{D}[\hat{O}_{pq}]$ on $\hat{X}_{n}$ and $\hat{Y}_{n}$ . After some algebra, we obtain

	$\displaystyle\sum_{\hat{O}\in C}\mathcal{D}[\hat{O}_{pq}]\hat{X}_{n}$	$\displaystyle=\left\{V_{pq}\mathcal{D}[\hat{Y}_{pq}]+W_{pq}\mathcal{D}[\hat{X}_{pq}]\right\}\hat{X}_{n},$		(42)
	$\displaystyle\sum_{\hat{O}\in C}\mathcal{D}[\hat{O}_{pq}]\hat{Y}_{n}$	$\displaystyle=\left\{V_{pq}\mathcal{D}[\hat{Y}_{pq}]+W_{pq}\mathcal{D}[\hat{X}_{pq}]\right\}\hat{Y}_{n},$		(43)

where the rates $V_{pq}$ and $W_{pq}$ are defined as

$\displaystyle V_{pq}$	$\displaystyle=\sum_{\hat{O}\in C}\left\|\bra{\psi_{p}^{L}}\hat{O}\ket{\psi_{q}^{R}}\right\|^{2}$	(44)
	$\displaystyle=\kappa(1+n_{\text{th}})\left\|\bra{\psi_{p}^{L}}\hat{a}\ket{\psi_{q}^{R}}\right\|^{2}+\kappa n_{\text{th}}\left\|\bra{\psi_{p}^{L}}\hat{a}^{\dagger}\ket{\psi_{q}^{R}}\right\|^{2},$	(45)
$\displaystyle W_{pq}$	$\displaystyle=\sum_{\hat{O}\in C}\left\|\bra{\psi_{p}^{R}}\hat{O}\ket{\psi_{q}^{R}}\right\|^{2}$	(46)
	$\displaystyle=\kappa(1+n_{\text{th}})\left\|\bra{\psi_{p}^{R}}\hat{a}\ket{\psi_{q}^{R}}\right\|^{2}+\kappa n_{\text{th}}\left\|\bra{\psi_{p}^{R}}\hat{a}^{\dagger}\ket{\psi_{q}^{R}}\right\|^{2},$	(47)

and represent the inter-well and intra-well transition rates, respectively.

Eqs. (42) to (45) shows that the action of the dissipator $\mathcal{D}[\hat{O}_{pq}]$ has two components, $V_{pq}\mathcal{D}[\hat{Y}_{pq}]$ and $W_{pq}\mathcal{D}[\hat{X}_{pq}]$ , with distinct physical meanings:

•

The component $V_{pq}\mathcal{D}[\hat{Y}_{pq}]$ describes direct inter-well transitions. This is because the collapse operator $\hat{Y}_{pq}$ (defined in Eq. (32)) flips the $x$ component of the Bloch vector, which measures the inter-well population difference. Corroborating our interpretation is the fact that the rate $V_{pq}$ defined in Eq. (45) only contains matrix element between states in different wells.

•

The component $W_{pq}\mathcal{D}[\hat{X}_{pq}]$ describes dephasing of superpositions between the right and left wells due to continuous measurement of “which-well” information by the environment (Fig. 4). To understand this statement, we trace the origin of $W_{pq}\mathcal{D}[\hat{X}_{pq}]$ back to the typical system-environment coupling Hamiltonian that gives rise to the Lindblad master equation [Eq. (14)]. The coupling Hamiltonian is

\displaystyle\hat{H}_{\text{int}}=\hat{a}\hat{B}^{\dagger}+\hat{a}^{\dagger}\hat{B}

(48)

where $\hat{B}$ and $\hat{B}^{\dagger}$ are environment operators that can be mathematically realized as bosonic raising and lowering operators. The projection of $\hat{H}_{\text{int}}$ onto transitions within the right- or left-wells are

	$\displaystyle\bra{\psi_{m}^{R}}\hat{H}_{\text{int}}\ket{\psi_{n}^{R}}$
	$\displaystyle\quad=\bra{\psi_{m}^{R}}\hat{a}\ket{\psi_{n}^{R}}\hat{B}^{\dagger}+\bra{\psi_{m}^{R}}\hat{a}^{\dagger}\ket{\psi_{n}^{R}}\hat{B}$		(49)
	$\displaystyle\bra{\psi_{m}^{L}}\hat{H}_{\text{int}}\ket{\psi_{n}^{L}}$
	$\displaystyle\quad=\bra{\psi_{m}^{L}}\hat{a}\ket{\psi_{n}^{L}}\hat{B}^{\dagger}+\bra{\psi_{m}^{L}}\hat{a}^{\dagger}\ket{\psi_{n}^{L}}\hat{B}$		(50)

Due to parity symmetry, the matrix elements $\bra{\psi_{m}^{R}}\hat{a}\ket{\psi_{n}^{R}}$ and $\bra{\psi_{m}^{L}}\hat{a}\ket{\psi_{n}^{L}}$ differ by a sign (and similarly for $\hat{a}^{\dagger}$ ). This means that the environment is driven in opposite directions depending on which well the state is in and slowly becomes entangled with the system. Therefore we can say that the environment is continuously measuring “which-well” information. Continuous measurement becomes dephasing in the right- and left-well state basis when we trace out the environment degrees of freedom. This explains our term of interest, $W_{pq}\mathcal{D}[\hat{X}_{pq}]$ , where the dissipator $\mathcal{D}[\hat{X}_{pq}]$ indeed dephases superpositions of right- and left-well states as per the definition of $\hat{X}_{pq}$ in Eq. (31). The dephasing rate $W_{pq}$ is related to matrix elements of $\hat{a}$ and $\hat{a}^{\dagger}$ between intra-well states.

Although the dephasing term $W_{pq}\mathcal{D}[\hat{X}_{pq}]$ does not directly cause any inter-well transitions, it interferes with coherent inter-well tunneling and results in the quantum Zeno effect misra_zenos_1977; itano_quantum_1990; facchi_quantum_2008. That is, in the limit where the environment measures “which-well” information much faster than the rate of coherent inter-well tunneling, the induced dephasing rate on the right- and left-well state basis is high and the system becomes frozen in either the right- or left-well.

At this point we have enough information to qualitatively understand the origin of the staircase as follows: the staircase is created by a competition between coherent tunneling driven by $\mathcal{L}_{\text{H}}$ , which generates inter-well superpositions, and dissipation driven by the continuous measurement from the environment, which dephases those superpositions. Depending on the relative rates of these two processes, each two-level manifold can be in two regimes, A and B, defined as follows:

•

In regime A, the tunnel splitting $\delta_{n}$ of the manifold is much larger than the intra-well transition rates, and a state in this manifold oscillates freely between the two wells. The contribution to the spontaneous switching rate $\Gamma$ from a manifold in this regime is limited by the transition rate into the manifold, and is insensitive to the value of the tunnel splitting $\delta_{n}$ . Therefore in this regime the switching rate changes slowly with the drive amplitude.
•

In regime B, the tunnel splitting $\delta_{n}$ is much smaller than the intra-well transition rates and thus the rate at which the environment measures “which-well” information. The quantum Zeno effect freezes the population in the initial well and the contribution to the spontaneous switching rate is suppressed. Switching in the manifold is limited by the residual rate of tunneling, which decreases as the tunnel splitting becomes exponentially suppressed. Therefore in this regime the switching rate decreases rapidly as the drive amplitude increases.

As a result of these two disparate regimes, at any fixed value of the two-photon drive amplitude $\epsilon_{2}$ and therefore depth of the double well, only the manifolds in regime A contribute significantly to the total spontaneous switching rate $\Gamma$ . Such manifolds must have energies close to or above the top of the energy barrier between the wells because the tunnel splitting is exponentially suppressed below the energy barrier. As $\epsilon_{2}$ increases, the double well becomes deeper and captures more excited state manifolds, and the tunnel splittings in all captured manifolds decrease exponentially. So long as a manifold stays in regime A, its contribution to the total switching rate is insensitive to the decreasing tunnel splitting. Eventually, the lowest lying manifold in regime A have such a small tunnel splitting that it abruptly transitions into regime B (Fig. 5), which is dominated by the quantum Zeno effect, and cease to contribute significantly to the spontaneous switching rate. As a result we see a step-like decrease in the total spontaneous switching rate $\Gamma$ . Each step-like decrease in the staircase therefore reflects one additional pair of levels that fall sufficiently below the energy barrier.

The physics contained in the qualitative description above should be captured by the complete set of equations of motion necessary for finding the switching rate. Carrying out the traces in Eqs. (38) and (39) produces the following equations of motion,

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\braket{\hat{X}_{n}}$	$\displaystyle=-\delta_{n}\braket{\hat{Y}_{n}}-\sum_{f}(W_{fn}+V_{fn})\braket{\hat{X}_{n}}$
	$\displaystyle\quad+\sum_{i}(W_{ni}-V_{ni})\braket{\hat{X}_{i}},$	(51)
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\braket{\hat{Y}_{n}}$	$\displaystyle=+\delta_{n}\braket{\hat{X}_{n}}-\sum_{f}(W_{fn}+V_{fn})\braket{\hat{Y}_{n}}$
	$\displaystyle\quad+\sum_{i}(V_{ni}-W_{ni})\braket{\hat{Y}_{i}}.$	(52)

We solve these equations of motion in the next subsection.

III.3 Perturbative solution of the effective equations of motion

In this subsection we extract the switching rate from the solution of Eqs. (51) and (52). As long as switching is the slowest dynamics, it will dominate in the long time limit where it manifests as an exponential decay in the population difference between the right and left wells. If we define the inter-well population difference as the expectation value of the following operator,

\displaystyle\hat{X}=\sum_{n}\hat{X}_{n},

(53)

then the spontaneous switching rate $\Gamma$ can be expressed as

\displaystyle\Gamma=-\frac{1}{\braket{\hat{X}}}\frac{\mathrm{d}}{\mathrm{d}t}\braket{\hat{X}},

(54)

where the expectation values are evaluated in the long time limit (identical as Eq.(16)). Substituting in Eq. (51) results in the following expression for $\Gamma$ ,

\displaystyle\Gamma=\sum_{n}\Gamma_{n},

(55)

where

\displaystyle\Gamma_{n}=\frac{1}{\sum_{m}\braket{\hat{X}_{m}}}\left[\delta_{n}\braket{\hat{Y}_{n}}+\left(2\sum_{f}V_{fn}\right)\braket{\hat{X}_{n}}\right].

(56)

The term $-\delta_{n}\braket{\hat{Y}_{n}}$ describes tunneling and the term $\left(2\sum_{f}V_{fn}\right)\braket{\hat{X}_{n}}$ describes direct inter-well transitions (as discussed in section III.2). To calculate the spontaneous switching rate, it remains to find $\braket{\hat{X}_{n}}$ and $\braket{\hat{Y}_{n}}$ via the equations of motion Eqs. (51) and (52).

III.3.1 Exponential ansatz in the long time limit

The first step towards solving for $\braket{\hat{X}_{n}}$ and $\braket{\hat{Y}_{n}}$ is to assume an ansatz for their time dependence. In the long time limit of interest, all expectation values decay exponentially at the switching rate $\Gamma$ towards their steady state values. Assuming a unique steady state, we must have $\braket{\hat{X}_{n}}\to 0$ and $\braket{\hat{Y}_{n}}\to 0$ in the long time limit by parity symmetry considerations. Therefore we can adopt the exponential ansatz, $\braket{\hat{X}_{n}},\braket{\hat{Y}_{n}}\propto e^{-\Gamma t}$ . This ansatz converts the equation of motion into an eigenvalue problem,

$\displaystyle-\Gamma\braket{\hat{X}_{n}}$	$\displaystyle=-\delta_{n}\braket{\hat{Y}_{n}}-\sum_{f}(W_{fn}+V_{fn})\braket{\hat{X}_{n}}$
	$\displaystyle\quad+\sum_{i}(W_{ni}-V_{ni})\braket{\hat{X}_{i}},$	(57)
$\displaystyle-\Gamma\braket{\hat{Y}_{n}}$	$\displaystyle=+\delta_{n}\braket{\hat{X}_{n}}-\sum_{f}(W_{fn}+V_{fn})\braket{\hat{Y}_{n}}$
	$\displaystyle\quad+\sum_{i}(V_{ni}-W_{ni})\braket{\hat{Y}_{i}}.$	(58)

We can solve the eigenvalue equations numerically to get $\Gamma$ , but to get an analytic expression we need to make additional approximations.

III.3.2 Solving for the quantum Zeno dynamics via the adiabatic elimination of $\braket{\hat{Y_{n}}}$

To obtain a solution that reflects the quantum Zeno effect that we qualitatively discussed in section III.2, now we solve the equation of motion for $\braket{\hat{Y}_{n}}$ , the coherence between the right- and left-well states. Recall that, while writing down the expression for the spontaneous switching rate $\Gamma$ in Eq. (54) and imposing the exponential ansatz for $\braket{\hat{X}}_{n}$ and $\braket{\hat{Y}}_{n}$ , we have assumed that spontaneous switching is the slowest dynamics and therefore the only process occuring in the long time limit. Based on this assumption, the system should be in a quasi-equilibrium state that changes slowly in time. Therefore, in Eq. (58), we can invoke the adiabatic approximation and ignore the term proportional to $\Gamma$ , which arose from the $\frac{\mathrm{d}}{\mathrm{d}t}\braket{\hat{Y}_{n}}$ (see appendix B for numerical verification). The condition for adiabaticity can be found by comparing $\Gamma\braket{\hat{Y}_{n}}$ to terms on the right-hand side of Eq. (58) as follows. If we denote the coefficient in front of $\braket{\hat{Y}_{n}}$ on the right-hand side of Eq. (58) as $\mu_{n}$ , i.e.,

\displaystyle\mu_{n}=2W_{nn}+\sum_{f\neq n}(W_{fn}+V_{fn}),

(59)

then Eq. (58) becomes

\displaystyle-\Gamma\braket{\hat{Y}_{n}}

\displaystyle=+\delta_{n}\braket{\hat{X}_{n}}-\mu_{n}\braket{\hat{Y}_{n}}+\sum_{i\neq n}(V_{ni}-W_{ni})\braket{\hat{Y}_{i}}.

(60)

A sufficient condition for the adiabaticity of $\braket{\hat{Y}_{n}}$ is then $\Gamma\ll\mu_{n}$ . For a sense of scale, the rate $\mu_{n}$ can be estimated to be $\mu_{n}\sim 2\kappa\alpha^{2}$ in the limit of $\alpha\gg 1$ and $n_{\text{th}}\ll 1$ (see appendix C [fix it so it makes sense]). We note that because $\mu_{n}$ controls the rate at which the coherence between the right- and the left-well states decays, it plays the central role of the measurement rate in the quantum Zeno effect.

As a result of the adiabatic approximation, the equation for $\braket{\hat{Y}_{n}}$ becomes independent of $\Gamma$ ,

\displaystyle 0

\displaystyle=\mu_{n}\braket{\hat{Y}_{n}}+\sum_{i\neq n}(W_{ni}-V_{ni})\braket{\hat{Y}_{i}}-\delta_{n}\braket{\hat{X}_{n}}.

(61)

We can derive an analytic approximate solution by rewriting Eq. (61) into the following recursive form,

\displaystyle\braket{\hat{Y}_{n}}

\displaystyle=\frac{1}{\mu_{n}}\left(\delta_{n}\braket{\hat{X}_{n}}-\sum_{i\neq n}(W_{ni}-V_{ni})\braket{\hat{Y}_{i}}\right).

(62)

If we treat transitions across manifolds (the second term in the parenthesis in the above expression) as a small perturbation, then the zeroth-order solution is

\displaystyle\braket{\hat{Y}_{n}}

\displaystyle\approx\frac{\delta_{n}}{\mu_{n}}\braket{\hat{X}_{n}}.

(63)

Physically, the zeroth-order solution describes a quasi-static balance between the generation of coherence due to tunneling and the dephasing of coherence due to measurement, and the manifold are decoupled. We can evaluate the error committed by the zeroth-order approximation and therefore its validity by recursively inserting the solution back into Eq. (62), yielding

	$\displaystyle\braket{\hat{Y}_{n}}$
	$\displaystyle=\frac{1}{\mu_{n}}\delta_{n}\braket{\hat{X}_{n}}$
	$\displaystyle-\sum_{i\neq n}\frac{1}{\mu_{n}}(W_{ni}-V_{ni})\frac{1}{\mu_{i}}\delta_{i}\braket{\hat{X}_{i}}$
	$\displaystyle+\sum_{i\neq n}\sum_{j\neq i}\frac{1}{\mu_{n}}(W_{ni}-V_{ni})\frac{1}{\mu_{i}}(W_{ij}-V_{ij})\frac{1}{\mu_{j}}\delta_{j}\braket{\hat{X}_{j}}$
	$\displaystyle+\cdots.$		(64)

In this series solution, the perturbative parameter $(W_{ni}-V_{ni})/\mu_{i}$ compares importance of cross-manifold transitions to that of tunneling and controls the size of the higher-order corrections. The estimates in Eqs. (111) and (112) lead to a bound on this perturbative parameter $(W_{ni}-V_{ni})/\mu_{i}=O(\alpha^{-2|i-n|})\ll 1$ for all $i\neq n$ in the limit of $\alpha\gg 1$ and $n_{\text{th}}\ll 1$ . This means that for a given manifold $n$ , the only way for the first-order correction to become significant compared to the zeroth-order approximation is when $\delta_{n}\braket{\hat{X}_{n}}\ll\delta_{i}\braket{\hat{X}_{i}}$ for at least some $i\neq n$ , that is, there is more tunneling happening in some other manifold. Therefore, we expect Eq. (63) to be a good approximation for, at the very least, the manifold that experiences the strongest tunneling, and errors committed in a manifold where tunneling is slow can be neglected. In appendix Bc we show numerically that the contribution $\Gamma_{n}$ to switching rate from each manifold is not significantly altered by the approximation, no matter if the manifold is the dominant contributor to switching. This allows us to use Eq. (63) uniformly for all manifolds.

With Eq. (63) we can re-express $\Gamma_{n}$ as written in Eq. (56) in terms of $\braket{\hat{X}_{n}}$ only,

\displaystyle\Gamma_{n}=\frac{1}{\sum_{m}\braket{\hat{X}_{m}}}\left[\frac{\delta_{n}^{2}}{\mu_{n}}\braket{\hat{X}_{n}}+\left(2\sum_{f}V_{fn}\right)\braket{\hat{X}_{n}}\right].

(65)

The term $(\delta_{n}^{2}/\mu_{n})\braket{\hat{X}_{n}}$ arises from the term $-\delta_{n}\braket{\hat{Y}_{n}}$ in Eq. (56), which describes tunneling. The rate $\delta_{n}^{2}/\mu_{n}$ can be regarded as a consequence of the quantum Zeno effect for the following reason. Suppose that we uncouple the system from the environment and instead directly perform projective measurements of $\hat{X}_{n}$ on it periodically at the interval of $\tau_{n}=2/\mu_{n}$ , and let the initial state be $\ket{\psi_{n}^{R}}$ in the right well. This creates a competition between the coherent tunneling rate and the incoherent measurement rate. In the segment of free evolution under the Hamiltonian $\hat{H}$ , the state coherently tunnels from the right well to the left well at the rate $\delta_{n}$ , creating an inter-well superposition. This superposition is periodically collapsed onto the right- or left-well states by the projective measurements. In the limit of $\delta_{n}\tau_{n}\ll 1$ , the transition probability between $\ket{\psi_{n}^{R}}$ and $\ket{\psi_{n}^{L}}$ is small, and the outcome of the projective measurements will show long sequences of the same outcome, randomly interrupted by infrequent transitions. On average, the expectation value of $\braket{\hat{X}_{n}}$ decays by a factor of $\cos(\delta_{n}\tau_{n})\approx\exp(-\delta_{n}^{2}\tau_{n}^{2}/2)$ per round of Hamiltonian evolution and projective measurement. Over time, the expectation value of $\hat{X}_{n}$ evolves as

\displaystyle\braket{\hat{X}_{n}}

\displaystyle\approx e^{-\frac{1}{2}\delta_{n}^{2}\tau_{n}^{2}\cdot\frac{t}{\tau_{n}}}=e^{-\frac{\delta_{n}^{2}}{\mu_{n}}t}.

(66)

We thus recover an incoherent tunneling rate of $\delta_{n}^{2}/\mu_{n}$ . The limit of $\mu_{n}\gg\delta_{n}$ corresponds to the well-known quantum Zeno effect where tunneling is frozen by frequent measurement.

III.3.3 Solving for the quasi-equilibrium population distribution via the adiabatic approximation on $\braket{\hat{X}_{n>0}}$

Our expression for the switching rate, Eq. (65), still depends on the unknown inter-well population difference, $\langle\hat{X}_{n}\rangle$ , in the long time limit. To solve for $\braket{\hat{X}_{n}}$ , we substitute Eq. (63) into Eq. (57) and eliminate $\braket{\hat{Y}_{n}}$ to get,

\displaystyle-\Gamma\braket{\hat{X}_{n}}

\displaystyle=-\frac{\delta_{n}^{2}}{\mu_{n}}\braket{\hat{X}_{n}}-\lambda_{n}\braket{\hat{X}_{n}}+\sum_{i\neq n}(W_{ni}-V_{ni})\braket{\hat{X}_{i}},

(67)

where we use $\lambda_{n}$ to denote the following rate

	$\displaystyle\lambda_{n}$	$\displaystyle=2V_{nn}+\sum_{f\neq n}(W_{fn}+V_{fn})$		(68)
		$\displaystyle=2\sum_{f}V_{fn}+\sum_{f\neq n}(W_{fn}-V_{fn}).$		(69)

$\lambda_{n}$ has the interpretation of the decay rate of population in the $n$ th manifold due to transitions to other manifolds in the same well or to the other well (but excluding the effect of tunneling).

To analytically solve Eq. (67), we can invoke the adiabatic approximation again for $\braket{\hat{X}_{n>0}}$ , leaving $\braket{\hat{X}_{0}}$ as the only independent degree of freedom (see appendix B for numerical verification). A sufficient condition for the adiabaticity of $\braket{\hat{X}_{n>0}}$ is $\Gamma\ll\lambda_{n}$ , i.e., the switching rate that we compute must also be much slower than the decay rate of populations in the well, which returns the system into an intra-well quasi-equilibrium state. Using the estimates for $W_{mn}$ and $V_{mn}$ from Eqs. (110) to (112) in the appendix, we arrive at an estimate of the decay rate $\lambda_{n}\sim\kappa n$ in the limit of $\alpha\gg 1$ and $n_{\text{th}}\ll 1$ , so our updated criterion for adiabaticity is $\Gamma\ll\kappa$ .

As a result of the adiabatic approximation, we can rewrite Eq. (67) for $\braket{\hat{X}_{n>0}}$ into the following form that separates upward transitions from downward transitions:

	$\displaystyle\braket{\hat{X}_{n>0}}$	$\displaystyle=\frac{\sum_{i<n}(W_{ni}-V_{ni})\braket{\hat{X}_{i}}}{\lambda_{n}+\frac{\delta_{n}^{2}}{\mu_{n}}}$
		$\displaystyle\quad+\frac{\sum_{i>n}(W_{ni}-V_{ni})\braket{\hat{X}_{i}}}{\lambda_{n}+\frac{\delta_{n}^{2}}{\mu_{n}}}.$		(70)

The first term consists of transitions coming from lower manifolds and the second coming from higher manifolds. To ascertain which term is more important, we should estimate the magnitude of each term. The estimates for $W_{mn}$ and $V_{mn}$ from Eqs. (110) to (112) shows that between two manifolds $m>n$ , the upward transition rate is always a factor of $n_{\text{th}}$ smaller than the downward transition rate. Also, the rate of upward transitions decreases exponentially with the degree of separation, $|m-n|$ , between the manifolds. Therefore, the quasi-equilibrium population differences should satisfy $\braket{\hat{X}_{0}}\gg\braket{\hat{X}_{1}}\gg\braket{\hat{X}_{2}}\gg\cdots$ . This means conditioned on the system being in some final state in a given manifold, it is more likely to have been excited from a lower manifold than decayed from an even higher manifold.

Based on this observation, we approximate Eq. (70) by ignoring downward transitions from higher manifolds as follows,

\displaystyle\braket{\hat{X}_{n>0}}\approx\frac{\sum_{i<n}(W_{ni}-V_{ni})\braket{\hat{X}_{i}}}{\lambda_{n}+\frac{\delta_{n}^{2}}{\mu_{n}}}.

(71)

(See appendix Bd for numerical verification of this approximation.) This gives us a triangular set of equations that we can solve one by one analytically, starting from $n=1$ ,

$\displaystyle\braket{\hat{X}_{1}}$	$\displaystyle=\frac{(W_{10}-V_{10})\braket{\hat{X}_{0}}}{\lambda_{1}+\frac{\delta_{1}^{2}}{\mu_{1}}},$	(72)
$\displaystyle\braket{\hat{X}_{2}}$	$\displaystyle=\frac{(W_{20}-V_{20})\braket{\hat{X}_{0}}}{\lambda_{2}+\frac{\delta_{2}^{2}}{\mu_{2}}}+\frac{(W_{21}-V_{21})\braket{\hat{X}_{1}}}{\lambda_{2}+\frac{\delta_{2}^{2}}{\mu_{2}}}.$	(73)
	$\displaystyle\cdots$

The first term of each equation captures excitations from the ground state manifold, and subsequent terms capture excitations from excited state manifolds. Using the following shorthand notation,

\displaystyle R_{mn}=\frac{W_{mn}-V_{mn}}{\lambda_{m}+\frac{\delta_{m}^{2}}{\mu_{m}}},\quad(m\neq n),

(74)

the explicit solution can be written down formally as the finite sum

\displaystyle\braket{\hat{X}_{n}}=\sum_{s=1}^{n}\left(\sum_{\begin{subarray}{c}0<i_{1}<i_{2}<\\ ...<i_{s-1}<s\end{subarray}}R_{ni_{s-1}}\cdots R_{i_{2}i_{1}}R_{i_{1}0}\right)\braket{\hat{X}_{0}}

(75)

where the inner sum goes over all $s$ -step paths connecting the ground-state manifold to the $n$ th manifold, and the outer sum goes over all possible number of steps.

We then compute the corrections due to downward transitions from higher manifolds perturbatively by substituting the solution Eq. (75) into the second term of Eq. (70) and solving the triangular set of equations again. The solution that contains corrections from downward transitions to order $c$ has the following form,

\displaystyle\braket{\hat{X}_{n}}=K_{n}\braket{\hat{X}_{0}}=\frac{K_{n}}{\sum_{m\geq 0}K_{m}},

(76)

where $K_{0}=1$ and $K_{n>0}$ is defined as

\displaystyle K_{n>0}=\sum_{\text{paths}}R_{ni_{s-1}}\cdots R_{i_{2}i_{1}}R_{i_{1}0}.

(77)

The sum now goes over all possible paths $0\mapsto i_{1}\mapsto i_{2}\mapsto\cdots\mapsto i_{s-1}\mapsto n$ that begins from the ground state manifold, does not return to the ground state manifold, ends at the $n$ th manifold, and involves at most 1 downward step. See appendix Bd for a numerical comparison between the approximations of Eq. (75) and Eq. (76).

Plugging the solutions Eq. (76) for the quasi-equilibrium population differences into Eq. (65) yields the following explicit expression for the contribution to the spontaneous switching rate from manifold $n$ ,

\displaystyle\Gamma_{n}=\left(\frac{\delta_{n}^{2}}{\mu_{n}}+2\sum_{f}V_{fn}\right)\frac{K_{n}}{\sum_{m\geq 0}K_{m}}.

(78)

The bracket has the dimension of rate and represents the switching rate for the population in the $n$ th manifold. The two terms in the parenthesis captures the effective tunneling rate due to the quantum Zeno effect and the transition rate due to quantum jumps, respectively. The term to the right of the square bracket weighs the rates for all the manifolds by the population in each manifold.

III.3.4 Final solution and summary of the derivation

While we have found the explicit solution for the contribution $\Gamma_{n}$ to the spontaneous switching rate from manifold $n$ [Eq. (78)], its current form obscures the physics that gives rise to the staircase. It isn’t clear which term is responsible for the step-like decrease in the total spontaneous switching rate $\Gamma$ as a function of the two-photon drive amplitude $\epsilon_{2}/K=\alpha^{2}$ , because the two factors in $\Gamma_{n}$ change simultaneously and in opposite directions. In particular, when we increase $\alpha^{2}$ past a critical value that allows one additional excited state manifold to fall sufficiently into the double well, $\delta_{n}$ decreases exponentially, but the normalized population difference in this manifold $\braket{\hat{X}_{n}}=K_{n}/[\sum_{m\geq 0}K_{m}]$ sharply increases precisely because of the exponential suppression of tunneling.

To elucidate the physics of the staircase contained in the expression of $\Gamma$ , we rewrite the expression for $\Gamma_{n}$ Eq. (78) into the following form for all excited state manifolds $n>0$ (see appendix D for derivation),

$\displaystyle\Gamma_{n}$	$\displaystyle=f_{n}J_{n}\quad\text{for }n>0,$	(79)
$\displaystyle f_{n}$	$\displaystyle=\frac{\frac{\delta_{n}^{2}}{\mu_{n}}+2\sum_{f}V_{fn}}{\frac{\delta_{n}^{2}}{\mu_{n}}+\lambda_{n}},$	(80)
$\displaystyle J_{n}$	$\displaystyle=\sum_{i\neq n}(W_{ni}-V_{ni})\frac{K_{i}}{\sum_{m\geq 0}K_{m}},$	(81)

where $f_{n}$ is the branching ratio into transitions that cause switching between the wells, and $J_{n}$ is (the inter-well difference of) the total probability current entering manifold $n$ . The difference between the two expressions for the switching rate contribution $\Gamma_{n}$ , Eqs. (78) and (79), is that the latter is expressed in terms of the populations in manifolds $i\neq n$ rather than the population of manifold $n$ itself. For a manifold $n$ that is being captured by the double well, the probability current $J_{n}$ entering it predominantly depends on the quasi-equilibrium population distribution in the double well, which is insensitive to changes near the top of the energy barrier (Fig. 6c). We therefore attribute the sharp drop in the spontaneous switching rate $\Gamma$ to the abrupt changes in the behavior of the branching ratios $f_{n}$ as functions of the two-photon drive amplitude $\epsilon_{2}/K=\alpha$ (Fig. 6b).

Thus we have expressed the spontaneous switching rate $\Gamma$ in terms of the tunnel splittings $\delta_{k}$ , the inter-well and intra-well transition rates $W_{pq}$ and $V_{pq}$ induced by $\hat{a}$ and $\hat{a}^{\dagger}$ , with the options to control the order of perturbation used in the computation of the quasi-equilibrium population distribution using an interger $c$ ,

	$\displaystyle\Gamma=\Gamma(\delta_{k},W_{pq},V_{pq})=\sum_{n}\Gamma_{n}(\delta_{k},W_{pq},V_{pq})$		(82)
	$\displaystyle\Gamma_{n}(\delta_{k},W_{pq},V_{pq})$
	$\displaystyle\quad=\begin{cases}\left(\frac{\delta_{0}^{2}}{\mu_{0}}+2\sum_{f}V_{f0}\right)\frac{K_{0}}{\sum_{m\geq 0}K_{m}}&n=0\\ f_{n}J_{n}&n>0\end{cases},$		(83)

Here is a summary of the assumptions we have made to obtain this result from the Lindblad master equation [Eq. (14)]:

1.

We projected the density operator onto two-level manifolds;
2.

We assumed that spontaneous switching is the slowest process in the system and therefore is the only dynamics in the long-time limit;
3.

We ignored transitions between manifolds when we calculated the effective tunneling rate due to the Quantum Zeno effect for each manifold.
4.

We invoked the adiabatic approximation to eliminate all other degrees of freedom apart from the inter-well population difference.
5.

Finally, when finding the quasi-equilibrium population differences $\braket{\hat{X}_{n>0}}$ , we ignored contributions that have multiple downward transitions from highly excited manifolds.

IV Location of the steps in the staircase

To locate the critical values of the normalized two-photon drive amplitude $\alpha=\epsilon_{2}/K$ where a sharp decrease in the spontaneous switching rate $\Gamma$ occurs, we plot the branching ratios $f_{n}$ , defined in Eq. (80), as functions of $\epsilon_{2}/K$ for excited state manifolds $n>0$ (Fig. 6b). At small $\alpha$ , tunneling is dominant and $(\delta_{n})^{2}/\mu_{n}\gg\lambda_{n}$ , so $f_{n}\approx 1$ , meaning all probability currents that enters an excited state manifold leads to spontaneous switching between wells. As $\alpha$ increases, one by one the branching ratio decreases. For a given manifold $n$ , we may define the critical normalized two-photon drive amplitude $\alpha_{\text{crit},n}$ to be the value of $\alpha$ at which the inter-well transition rate and the intra-well transition rate become equal, or equivalently, the value of $\alpha$ at which $f_{n}=1/2$ . That is, $\alpha_{\text{crit},n}$ is such that

\displaystyle f_{n}=\frac{\frac{\delta_{n}^{2}}{\mu_{n}}+2\sum_{f}V_{fn}}{\frac{\delta_{n}^{2}}{\mu_{n}}+\lambda_{n}}=\frac{1}{2}.

(84)

Comparing Fig. 6a and 6b shows that the location of the steps in the staircase agrees with this condition. To obtain a simpler expression for the condition on $\alpha_{\text{crit},n}$ , we notice from Fig.6 that when $f_{n}=\frac{1}{2}$ , the branching ratio can be approximated by the expression

\displaystyle f_{n}\approx\frac{\frac{\delta_{n}^{2}}{\mu_{n}}}{\frac{\delta_{n}^{2}}{\mu_{n}}+\lambda_{n}}.

(85)

Using this approximate expression, we find, after some manipulation, that $\alpha_{\text{crit},n}$ is the solution to the following equation

\displaystyle\frac{\delta_{n}^{2}(\alpha)}{\mu_{n}(\alpha)}=\lambda_{n}(\alpha).

(86)

By using the estimates of appendix C (and potentially WKB estimates of $\delta_{n}$ in appendix E), we can rewrite this condition as

\displaystyle\delta_{n}^{2}(\alpha)\approx 2\kappa^{2}\alpha^{2}.

(87)

In particular, $\alpha_{\text{crit},1}$ , determined by the special case $n=1$ of the condition above, increases as the damping rate $\kappa$ decreases. Therefore, in this limit, the critical $\alpha$ beyond which there is an exponential suppression of the spontaneous switching rate $\Gamma$ will be delayed, as seen in Fig. 2.

V Application of the semi-analytical formula: activation mechanism within the well

In this section, we use our semi-analytical formula to investigate which heating transitions are the dominant ones that populate the excited-state manifolds. We refer to the collection of the dominant heating transitions as the activation mechanism. The activation mechanism can be deduced from the sensitivity of our semi-analytical formula to the matrix elements of $\hat{a}$ and $\hat{a}^{\dagger}$ , which govern the rates of in- and cross-manifold transitions. Each jump operator can cause transitions within the same manifold, towards lower manifolds (cooling), or towards higher manifolds (heating), and some are dominant over others. First, $\hat{a}$ dominates over $\hat{a}^{\dagger}$ for in-manifold transitions because $n_{\text{th}}\ll 1$ . Next, among the cooling transitions, we can neglect all those induced by $\hat{a}^{\dagger}$ because $n_{\text{th}}\ll 1$ and the relevant matrix elements are much smaller than the corresponding ones for $\hat{a}$ -induced cooling transitions and (see estimates of appendix C). The remaining cooling transitions are those induced by $\hat{a}$ , among which we can also neglect the ones that lower the manifold index by more than 1, for the following reason. The probability current flowing from some initial manifold $q$ into a fixed final manifold $p$ due to $\hat{a}$ -induced cooling transitions is the product of the matrix element squared $|\bra{\psi_{p}^{R}}\hat{a}\ket{\psi_{q}^{R}}|^{2}$ and the population in the state $\ket{\psi_{q}^{R}}$ . Both the matrix element squared (see estimates in appendix C) and the population are exponentially suppressed with increasing $q$ . Therefore the dominant $\hat{a}$ -induced cooling transitions only occur between consecutive manifolds.

Among the heating transitions, those induced by $\hat{a}$ persist in the limit $n_{\text{th}}\to 0$ and have been referred to as quantum heating marthaler_switching_2006. They contrast with heating due to the presence of thermal photons in the environment, that is, the $\hat{a}^{\dagger}$ -induced heating transitions. To verify the assumptions we made about the cooling transitions and to identify the which heating transitions are dominant, we consider two scenarios where we set different matrix elements in the inter-well and intra-well transition rates $V_{pq}$ , $W_{pq}$ to zero and see what effect the change has on the switching rate $\Gamma$ computed using the semianalytical formula. Neither of the two scenarios will allow quantum heating. In the first scenario, the modified inter-well and intra-well transition rates have the following non-zero elements,

	$\displaystyle V_{pq}^{(\text{c.t.})}=\begin{cases}\kappa(1+n_{\text{th}})\left\|\bra{\psi_{p}^{L}}\hat{a}\ket{\psi_{q}^{R}}\right\|^{2}&q-1\leq p\leq q\\ \kappa n_{\text{th}}\left\|\bra{\psi_{p}^{L}}\hat{a}^{\dagger}\ket{\psi_{q}^{R}}\right\|^{2}&q<p\leq q+d\\ 0&\text{otherwise}\\ \end{cases},$		(88)
	$\displaystyle W_{pq}^{(\text{c.t.})}=\begin{cases}\kappa(1+n_{\text{th}})\left\|\bra{\psi_{p}^{R}}\hat{a}\ket{\psi_{q}^{R}}\right\|^{2}&q-1\leq p\leq q\\ \kappa n_{\text{th}}\left\|\bra{\psi_{p}^{R}}\hat{a}^{\dagger}\ket{\psi_{q}^{R}}\right\|^{2}&q<p\leq q+d\\ 0&\text{otherwise}\\ \end{cases},$		(89)

where $d$ is some positive integer. This set of transition rates excludes all quantum heating and any thermal heating that raises the manifold index by more than $d$ . This means that to reach a highly excited state, multiple thermal heating events are required if $d$ is small. We refer to this set of rates as the cascaded thermal heating scenario, hence the superscript “(c.t.)”. The second set of modified inter-well and intra-well transition rates we define have the following non-zero elements,

	$\displaystyle V_{pq}^{(\text{d.t.})}=\begin{cases}\kappa(1+n_{\text{th}})\left\|\bra{\psi_{p}^{L}}\hat{a}\ket{\psi_{q}^{R}}\right\|^{2}&q-1\leq p\leq q\\ \kappa n_{\text{th}}\left\|\bra{\psi_{p}^{L}}\hat{a}^{\dagger}\ket{\psi_{0}^{R}}\right\|^{2}&q=0\text{ and }p>0\\ 0&\text{otherwise}\\ \end{cases},$		(90)
	$\displaystyle W_{pq}^{(\text{d.t.})}=\begin{cases}\kappa(1+n_{\text{th}})\left\|\bra{\psi_{p}^{R}}\hat{a}\ket{\psi_{q}^{R}}\right\|^{2}&q-1\leq p\leq q\\ \kappa n_{\text{th}}\left\|\bra{\psi_{p}^{R}}\hat{a}^{\dagger}\ket{\psi_{0}^{R}}\right\|^{2}&q=0\text{ and }p>0\\ 0&\text{otherwise}\\ \end{cases}.$		(91)

This set of transition rates excludes all quantum heating and any thermal heating that does not originate from the ground state manifold. This means that reaching a highly excited state requires a single thermal heating event from the ground state to the target manifold. We refer to this set of rates as the direct thermal heating scenario, hence the superscript “(d.t.)”. Then the switching rates predicted by the semianalytical formula under cascaded thermal heating or under direct thermal heating are defined as follows:

	$\displaystyle\Gamma^{(\text{c.t.})}$	$\displaystyle=\Gamma(\delta_{n},V_{pq}^{(\text{c.t.})},W_{pq}^{(\text{c.t.})}),$		(92)
	$\displaystyle\Gamma^{(\text{d.t.})}$	$\displaystyle=\Gamma(\delta_{n},V_{pq}^{(\text{d.t.})},W_{pq}^{(\text{d.t.})}).$		(93)

In Fig. 7 we compare $\Gamma^{(\text{c.t.})}$ and $\Gamma^{(\text{d.t.})}$ to the numerically obtained switching rate. Fig. 7(a) shows that to account for the switching rate at $n_{\text{th}}=10^{-3}$ up to the $k^{\text{th}}$ step in the staircase, it is necessary and sufficient to select a the cutoff $d$ in the definition for cascaded thermal heating that is at least $d=k$ . That is to say, thermal heating events must be allowed to raise the manifold index by at least $k$ , the exact number required to directly excite population from the ground state manifold to the $k^{\text{th}}$ excited state manifold. We conclude that at low $n_{\text{th}}$ , the activation mechanism is direct thermal heating rather than cascaded thermal heating or quantum heating. Fig. 7(a) further shows that to account for the switching rate at $n_{\text{th}}=10^{-1}$ , it is necessary and sufficient to choose a cutoff $d$ of 3. We therfore conclude that at moderate $n_{\text{th}}$ , the activation mechanism are cascaded thermal heating events in which the manifold index is raised by 3 or less at a time, rather than direct thermal heating or quantum heating. The difference in activation mechanism for different $n_{\text{th}}$ is expected because of the population in intermediate excited states are higher at larger $n_{\text{th}}$ . Our conclusions are supported by Fig. 7(b), where the semianalytical switching rate $\Gamma^{(d.t.)}$ under only direct thermal heating agrees with exact numerics for low $n_{\text{th}}$ but fails for moderate $n_{\text{th}}$ .

In Fig. 7(c), we examine the scaling of the spontaneous switching rate as a function of the ratio $n_{\text{th}}/(1+n_{\text{th}})$ . This particular ratio measures the relative strength of thermal heating and $\hat{a}$ -induced cooling, which comprise the dominant cross-manifold transitions in the limit $n_{\text{th}}\ll 1$ . By choosing values of $\epsilon_{2}/K$ that sit at the $k^{\text{th}}$ step in the staircase, we can fix the dominant manifold in which switching occurs to be manifold $k$ . Then, for fixed $\epsilon_{2}/K$ , if the switching rate scales as $\Gamma\propto[n_{\text{th}}/(1+n_{\text{th}})]^{s}$ , then there are roughly $s$ degrees of separation between the ground state manifold and manifold $k$ , measured by the average number of thermal heating events required to complete the excitation. An alternative way to interpret the scaling is as follows. If $n_{\text{th}}$ is taken to be the average thermal photon number in the environment at the frequency $\omega_{0}$ , i.e. $n_{\text{th}}=(e^{\hbar\omega_{0}/k_{B}T}-1)^{-1}$ , then the ratio is $[n_{\text{th}}/(n_{\text{th}}+1)]^{s}=\exp(-s\hbar\omega_{0}/k_{B}T)$ . This scaling implies an activation energy of $s\hbar\omega_{0}$ .

In Fig. 7(c) shows that for $n_{\text{th}}/(1+n_{\text{th}})\lesssim 10^{-9}$ and $\epsilon_{2}/K=5$ , chosen such that the dominant switching manifold is $k=1$ , the degree of separation is $s=0$ . This is because $n_{\text{th}}$ is so low that the switching occurs in the ground state manifold and no thermal heating events are required. For $10^{-9}\lesssim n_{\text{th}}/(1+n_{\text{th}})\lesssim 10^{-3}$ , the slope of the red curves suggests there is only one degree of separation between the ground state manifold and any of the first seven excited state manifolds probed by the chosen values of $\epsilon_{2}/K$ . This confirms our conclusion that direct thermal heating is dominant for low $n_{\text{th}}$ . Finally, for $10^{-3}\lesssim n_{\text{th}}/(1+n_{\text{th}})<1$ , the degree of separation $s$ increases beyond 1 but stays below the manifold index $k$ (with the exception of $\epsilon_{2}/K=5$ , for which $s>k$ likely because of corrections higher order in $n_{\text{th}}$ ). This means that the path of excitation from the ground state manifold to the $k^{\text{th}}$ excited state manifold is broken down to smaller steps of average size larger than one. This confirms that cascaded thermal heating is dominant for moderate $n_{\text{th}}$ .

VI Comparison to previous work

In this section we compare our work to previous works where theoretical predictions about the switching rate were obtained. We will first focus the comparison to theoretical works that deal specifically with the two-photon driven Kerr nonlinear oscillator at zero detuning under the influence of single photon loss and gain but no dephasing, noting that the situation with detuning or dephasing has been studied in marthaler_switching_2006; frattini_squeezed_2022; ruiz_two-photon_2022; venkatraman_driven_2024. Then, we broaden our scope and compare our results to previous theory of resonant tunneling in the magnetization of molecular magnets garanin_thermally_1997. Finally, we will discuss how our result relates to the famous Kramers problem (see review in hanggi_reaction-rate_1990).

The most important distinction of this work from previous work on the two-photon driven Kerr nonlinear oscillator is our identification of the role played by the measurement-induced dephasing rate $\mu_{n}$ in determining the switching-rate $\Gamma_{n}$ . This was previously overlooked in marthaler_quantum_2007, gautier_combined_2022 and appendix I of putterman_stabilizing_2022. The scaling $\mu_{n}\sim 2\kappa\alpha^{2}$ means that for increasing $\alpha$ the neglect of measurement-induced dephasing will lead to an increasing error in the switching rate $\Gamma$ .

A related distinction of this work from previous work is that our condition for the critical manifold $n_{\text{crit}}$ that dominates switching depends on the measurement-induced dephasing rate $\mu_{n}$ and the decay rate $\lambda_{n}$ , both of which depend on the parameters of dissipation. In putterman_stabilizing_2022 and marthaler_quantum_2007 however, $n_{\text{crit}}$ was asserted to be the manifold that lies at the top of the energy barrier between the double wells. From our semi-analytics we see that the assertion need not be the case. Using Eq. (87) we find that for a fixed $\alpha$ , the critical manifold $n_{\text{crit}}$ decreases with $\kappa$ , so it is possible for a manifold significantly below the energy barrier to dominate switching. When this manifold sits between the top and bottom of the energy well, the mechanism is typically known as thermally-assisted tunneling (see for example harris_thermally_1991). Indeed in gautier_combined_2022, putterman_stabilizing_2022 and in our Fig. 2(b), when $\kappa$ takes on the relatively small value of $10^{-3}K$ or $10^{-4}K$ , the switching rate $\Gamma$ remains on the first plateau due to switching in the first excited state for $\epsilon_{2}/K=\alpha^{2}$ as large as 8 or 9 despite an estimated the number of states in the well $\alpha^{2}/\pi\approx 2.5$ , and the onset of exponential suppression of the switching rate is delayed.

A particular difference between our work and gautier_combined_2022 is that, much like marthaler_switching_2006, we compute the switching rate after the intra-well quasi-equilibrium population distribution has already been reached through thermalization and relaxation. This distribution holds at all $t$ if the initial state is already the thermalized state, or at $t\gg\kappa^{-1}$ for a general initial state. However, gautier_combined_2022 computes the average switching rate for an initial state confined to the ground state in the right-well, and only in the the transient limit, i.e. $t\ll\kappa^{-1}$ , when the quasi-equilibrium has yet to be reached. The latter choice has interesting consequences. The short-time average switching rate increases as the excited state population builds up, and at a fixed time $t=t_{\text{eval}}$ the average switching rate also has a staircase dependence as a function of the two-photon drive amplitude. We emphasize that this alternative staircase is not the same as the one we consider in our work. In particular, the condition that determines the critical $\alpha$ at which the steps occur is

\displaystyle\delta_{n}(\alpha)=\frac{\pi}{t_{\text{eval}}},

(94)

and is independent of parameters of the dissipation. We provide two justifications for our choice of the switching rate to compute. First, in the experiments frattini_squeezed_2022; hajr_high-coherence_2024; ding_quantum_2024, the switching rate was extracted from fitting the exponential decay of the inter-well population difference over a wait time that is much longer than $\kappa^{-1}$ , so the measured switching rate is dominated by the dynamics after the intra-well quasi-equilibrium has been reached. Second, because the measurement protocol (known as cat-quadrature readout grimm_stabilization_2020) used in these experiments projects the system onto either the right- or left-well without disturbing the intra-well population distribution, the quasi-equilibrium distribution will persist over many rounds of measurement and play the role of the initial state once it has been reached (assuming constant heating and cooling rates over time).

The activation mechanisms we find responsible for exciting the population to the critical manifold $n_{\text{crit}}$ are different from previous work. gautier_combined_2022 and putterman_stabilizing_2022 only considered direct thermal heating, whereas marthaler_quantum_2007 finds quantum heating to be dominant for general values of the detuning between half the two-photon drive frequency and the oscillator resonant frequency (see also experimental observation of quantum heating in vijay_invited_2009; ong_quantum_2013 in the semi-classical regime of many levels in the well). Our work shows that in our regime, quantum heating from the ground and excited states can be neglected, and both direct and cascaded thermal heating are responsible for activation.

Finally, we note some minor distinctions of this work from previous work that allow us to identify the staircase in the switching rate. Unlike puri_bias-preserving_2020, we do not include a strong two-photon loss term $\kappa_{2}\mathcal{D}[\hat{a}^{2}]$ in the Lindbladian, which otherwise would strongly suppress the population in all excited state manifolds and hide staircase. Another distinction is that just like gautier_combined_2022, we do not make a continuous approximation on the critical manifold $n_{\text{crit}}$ as a function of $\epsilon_{2}/K=\alpha^{2}$ . This continuous approximation typically smooths out the staircase and only preserves the overall trend in the switching rate $\Gamma$ as a function of $\epsilon_{2}/K$ as was the case in marthaler_switching_2006; putterman_stabilizing_2022.

We now compare our result to previous theory on resonant magnetization tunneling in molecular magnets garanin_thermally_1997. The mathematical technique used in this work is similar to that in garanin_thermally_1997, where the analytically computed switching rate between two magnetisation states $\ket{m=-S}$ and $\ket{m=+S}$ of a molecular spin $S$ displayed staircase-like dependence on the transverse magnetic field strength (see also experimental observation in adams_geometric-phase_2013). When the longitudinal magnetic field is zero, the Hamiltonian of the spin is $\hat{H}=-D\hat{S}_{z}^{2}-H_{x}\hat{S}_{x}$ , whose energy levels in the zero transverse field limit are similar to the energy levels of a symmetric double well. The transverse field $H_{x}$ induces resonant tunneling between pairs of opposite magnetization states $m^{\prime}=-m$ at rates rates $\Omega_{m,m^{\prime}}$ that are exponentially suppressed with decreasing $m$ . As $H_{x}$ increases, fast resonant tunneling between a pair of states $m_{\text{top}},m_{\text{top}}^{\prime}$ below the top of the original energy barrier effectively reduces the depth of the double well friedman_quantum_1998, so $H_{x}$ in the spin system plays the role of $\epsilon_{2}$ in our work, albeit in opposite directions. The distinction between the switching dynamics of the spin and that found in our work is as follows. The coupling between the spin and the environment does not cause transitions that change the manifold index by more than 1 or self-transitions within the same two-level manifold. As a consequence, the activation mechanism consists only of cascaded thermal activation, leading to a rather simple expression for the intra-well quasi-equilibrium distribution. Furthermore, the condition that determines the critical manifold $m_{b}$ is simply

\displaystyle\Omega_{m_{b},m_{b}^{\prime}}=\Gamma_{m_{b},m_{b}^{\prime}},

(95)

[c.f. Eq. (86)] where $\Omega_{m_{b},m_{b}^{\prime}}$ is the resonant tunnel splitting between magnetization states $m_{b}$ and $m_{b}^{\prime}=-m_{b}$ , and $\Gamma_{m_{b},m_{b}^{\prime}}$ is the decay rate from the two-level manifold consisting of $m_{b}$ and $m_{b}^{\prime}$ .

Lastly, due to its historical relevance, let us close this section by discussing how our result relates to the famous Kramers problem kramers_brownian_1940; hanggi_reaction-rate_1990. One phenomenon of fundamental importance is the crossover from the classical regime to the quantum regime where the effect of discrete energy levels and quantum tunneling is visible (see for example andersen_quantum_2020; yamaji_spectroscopic_2022). Fig. 8a shows that increasing $\kappa/K$ to $1$ smooths out the staircase compared to $\kappa/K=0.025$ (see also frattini_squeezed_2022), so we can consider the $\kappa\approx K$ to be the crossover value of damping. Using our semi-analytical formula for the switching rate, we identify that the smoothing is due to the following two effects: compared to Fig. 6, the branching ratio $f_{n}$ decreases slower in the dissipation-dominated regime at $\kappa/K=1$ , and the excitation current $J_{n}$ decreases faster in the tunneling-dominated regime.

The dependence of the switching dynamics on the damping rate $\kappa$ reveals another effect that is similar to the classical effect known as Kramers turnover (see Fig. 9). At fixed $\epsilon_{2}$ , which holds fixed the size of the energy barrier, our semi-analytical switching rate as a function of the damping rate $\kappa$ exhibits many local maxima. Consider one such maxima in switching rate, dominated by the contribution from a certain manifold $n$ . In the limit of small $\kappa$ , the switching rate contributed by this manifold is limited by the small thermal excitation current entering the manifold, so increasing $\kappa$ increases the switching rate, much like in the underdamped regime of Kramers problem. In the limit of large $\kappa$ , the correspondingly larger measurement-induced dephasing rate $\mu_{n}$ and the decay rate $\lambda_{n}$ hinders tunneling in manifold $n$ , so increasing $\kappa$ lower the switching rate, much like in the overdamped regime of Kramers problem. This turnover repeats in each excited state manifold.

Finally, a well-known effect from the research on the Kramers problem is the existence of a crossover temperature below which quantum tunneling from the ground state is important. Our results (Fig. 7) shows that in our system, the temperature dependence of the switching rate has two crossovers instead of one. In the limit of low $n_{\text{th}}$ , switching remains dominated by quantum jumps from the ground state on one well to any state in the other well. As $n_{\text{th}}$ increases, switching dominantly happens in an excited state manifold populated via direct thermal heating. Further increases in $n_{\text{th}}$ leads to another crossover, where the same excited state manifold is being populated by cascaded thermal heating rather than direct thermal heating.

VII Conclusions

We have derived a semi-analytical formula [Eq. (82)] for the spontaneous switching rate in the double well of a two-photon driven Kerr nonlinear oscillator. From this formula we find that tunneling dynamics in each excited state manifold contributes to one step in the staircase pattern of the switching rate as a function of the driven amplitude. There is no oscillatory behaviour in the switching rate because via adiabatic elimination, we focus on the slow dynamics after the system has relaxed to a quasi-equilibrium intra-well state.

The origin of the staircase is as follows. First, at large enough drive amplitudes, switching dominantly occurs in one excited state manifold rather the ground state manifold. Next, the competition between tunneling and dissipation within this excited state manifold leads to a tunneling-dominated regime and a dissipation-dominated regime, and the two regimes are traversed in this order as the drive amplitude increases. In the tunneling-dominated regime, the tunneling splitting is large compared to the rate of dissipation, and limiting factor of the switching rate is the rate at which the excited manifold is populated. Therefore the switching rate is insensitive to the size of the tunnel splitting and varies slowly with the drive amplitude. In the dissipation-dominated regime, the tunnel splitting is small compared to the rate of dissipation, and the limiting factor for the switching rate is the effective tunneling rate $\delta_{n}^{2}/\mu_{n}$ due to the quantum Zeno effect freezing out the tunneling dynamics. Therefore the switching rate rapidly decreases as a function of the drive amplitude due to the exponential suppression of the tunneling splitting when the manifold is captured in the double well. The two regimes give rise to a flat section and a steep section in the switching rate as a function of drive amplitude. This pattern repeats when switching in a higher excited manifold becomes the dominant at larger drive amplitudes, producing a characteristic staircase pattern.

We also deduced the following using our semi-analytical formula. For each excited state manifold, the critical drive amplitude that separates the two aforementioned regimes is determined from a condition [Eq. (86)] relating the tunneling splittings, the measurement induced dephasing rates, and the decay rates of the manifolds. At fixed drive amplitude, the same condition can also be used to determine the excited state manifold that dominates switching. Finally, in the regime we have examined, the activation mechanism that transports population from the ground state to the excited state that dominates switching is direct and cascaded thermal heating rather than quantum heating. Our theory provides a fine level of detail on the switching dynamics in a double well potential in the quantum regime.

The analysis in the work applies largely to parametric oscillators with bi-stability, and with some modifications, to coherent dissipative tunneling in fluxonium qubits and chemical reactions. One direction for future work is to apply this analysis to when the two-photon drive is detuned ruiz_two-photon_2022; gravina_critical_2023; venkatraman_driven_2024 and when a single-photon drive breaks the inversion symmetry of the double well bones_resonant-force_2024.

VIII Acknowledgments

We thank Max Schäfer, Andy Z. Ding, Benjamin L. Brock, Zhaoyou Wang, and Jahan Claes for fruitful discussions. [Funding acknowledgments]

Appendix A Numerical solution of the spontaneous switching rate

In this section we describe the numerical procedure used to find the total spontaneous switching rate $\Gamma$ and the contribution from each manifold $\Gamma_{n}$ . First, notice that the Lindbladian $\mathcal{L}$ defined in Eq. (14) obeys the symmetry $\hat{\Pi}\mathcal{L}(\cdot)\hat{\Pi}^{\dagger}=\mathcal{L}(\hat{\Pi}(\cdot)\hat{\Pi}^{\dagger})$ where $\hat{\Pi}=\exp(i\pi\hat{a}^{\dagger}\hat{a})$ is the parity transformation operator. This means that the density operator can be written as the sum of two independently evolving components,

\displaystyle\hat{\rho}(t)=\hat{\rho}_{+}(t)+\hat{\rho}_{-}(t),

(96)

where $\hat{\Pi}\hat{\rho}_{+}(t)\hat{\Pi}=\hat{\rho}_{+}(t)$ and $\hat{\Pi}\hat{\rho}_{-}(t)\hat{\Pi}=-\hat{\rho}_{-}(t)$ . That is, $\hat{\rho}_{+}(t)$ is even and $\hat{\rho}_{-}(t)$ is odd under the parity transformation. We are only interested in the evolution of $\hat{\rho}_{-}(t)$ because the inter-well population difference is inverted by the parity transformation. This evolution is generated by $\mathcal{L}_{-}$ , the restriction of $\mathcal{L}$ onto the subspace of odd parity operators. Next, rather than following the time evolution, we extract the switching rate by looking at the eigenvalues of $\mathcal{L}_{-}$ . We only consider the real eigenvalues because we expect spontaneous switching to cause exponential decay of the inter-well population difference with no oscillations. Finally, when spontaneous switching is the slowest dynamics in the double well, we can identify $\Gamma$ as the negative of the real eigenvalue of $\mathcal{L}_{-}$ that is closest to 0.

Now, we numerically compute the contribution $\Gamma_{n}$ to the spontaneous switching rate from each manifold in the following way. First we define a quasi-equilibrium state $\hat{\rho}_{\text{eq}}$ as

\displaystyle\hat{\rho}_{\text{eq}}=\hat{\rho}_{\text{ss}}+c\hat{\rho}_{\text{ex}},

(97)

where $\hat{\rho}_{\text{ss}}$ is the steady state of $\mathcal{L}$ and $c$ is some arbitrary real number. Next we substitute the quasi-equilibrium state into the definition for $\Gamma$ Eq. (54),

	$\displaystyle\Gamma$	$\displaystyle=-\frac{1}{\mathrm{Tr}(\hat{X}\hat{\rho}_{\text{eq}})}\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{Tr}(\hat{X}\hat{\rho}_{\text{eq}})$		(98)
		$\displaystyle=-\frac{\mathrm{Tr}[\hat{X}\mathcal{L}(\hat{\rho}_{\text{eq}})]}{\mathrm{Tr}(\hat{X}\hat{\rho}_{\text{eq}})}$		(99)

If we approximate the equilibrium state as an incoherent mixture between different manifolds, that is

\displaystyle\hat{\rho}_{\text{eq}}\approx\sum_{n}\hat{P}_{n}\hat{\rho}_{\text{eq}}\hat{P}_{n},

(100)

then the total spontaneous switching rate can be approximated as a sum over $\Gamma_{n}$ , defined as

\displaystyle\Gamma_{n}=-\frac{\mathrm{Tr}[\hat{X}\mathcal{L}(\hat{P}_{n}\hat{\rho}_{\text{eq}}\hat{P}_{n})]}{\mathrm{Tr}(\hat{X}\hat{\rho}_{\text{eq}})}

(101)

Because the steady state is symmetric under parity transformation, it does not contribute to the inter-well population difference, and $\Gamma_{n}$ can be equivalently expressed as

\displaystyle\Gamma_{n}=-\frac{\mathrm{Tr}[\hat{X}\mathcal{L}(\hat{P}_{n}\hat{\rho}_{\text{ex}}\hat{P}_{n})]}{\mathrm{Tr}(\hat{X}\hat{\rho}_{\text{ex}})}.

(102)

Appendix B Validity of approximations used in the derivation

Here we numerically demonstrate the validity of some of the approximations used in the main text.

Projection onto two-level manifolds

The projection onto two-level manifolds performed in Eq. (17) results in an effective Lindbladian $\mathcal{L}_{\text{eff}}$ defined in Eq. (24). Fig. 10 compares the switching rate $\Gamma$ obtained using $\mathcal{L}$ and $\mathcal{L}_{\text{eff}}$ . The discrepancy in the $n_{\text{th}}=0$ limt can be fixed using the eigenstates of the non-orthonormal eigenstates of the non-hermitian effective Hamiltonian, as in Fig. 11.

An additional limitation comes from the assumptions we make when we solve for the switching rate $\Gamma$ . One of the first assumptions we make to derive the formula is the validity of Eq. (17), the projection of the density operator $\hat{\rho}$ onto two-level manifolds. This is generally a good approximation for $n_{\text{th}}<1$ , but breaks down in the limit of extremely small $n_{\text{th}}$ (see appendix Ba). This is because although coherent states are eigenstates of $\hat{a}$ , the ground states of the effective Hamiltonian $\hat{H}_{\text{eff}}=\hat{H}-\frac{i\kappa}{2}\hat{a}^{\dagger}\hat{a}$ are perturbed away from coherent states in the presence of non-zero $\kappa$ . If the right-well ground state of the effective Hamiltonian is denoted by $\ket{\phi}$ , then by perturbation theory, we can estimate

\displaystyle\ket{\phi}=\ket{\alpha}+O(\kappa).

(103)

Under the collapse operator $\hat{O}=\sqrt{\kappa}\hat{a}$ corresponding to single photon loss at $n_{\text{th}}=0$ , the scaling of the total heating rate out of the ground state as a function of $\kappa$ can be found as follows

	$\displaystyle\lVert(1-\ket{\phi}\bra{\phi})\hat{O}\ket{\phi}\rVert^{2}$
	$\displaystyle=\lVert[1-\ket{\alpha}\bra{\alpha}+O(\kappa)]\,[\alpha\sqrt{\kappa}\ket{\alpha}+O(\kappa^{3/2})]\rVert^{2}$		(104)
	$\displaystyle=\|O(\kappa^{3/2})\|^{2}$		(105)
	$\displaystyle=O(\kappa^{3}).$		(106)

Our semianalytical formula can be extented to instead use the perturbed eigenstates of the effective Hamiltonian

Adiabatic approximation

The adiabatic approximation happens in Eqs. (61) and (70) in the main text, interspersed with other approximations. To demonstrate its validity, here we examine the adiabatic approximation in isolation. Our starting point is Eq. (51) for $n>0$ and Eq. (52) for all $n$ , which are obtained through projection onto two-level manifolds. Now, by setting the time derivative of all $\hat{Y}_{n}$ and $\hat{X}_{n>0}$ to zero, we can numerically solve for $\hat{Y}_{n}$ and $\hat{X}_{n>0}$ in terms of $\hat{X}_{0}$ . Finally, we can plug the result into Eq. (56), where $\hat{X}_{0}$ will cancel, and obtain $\Gamma_{n}$ under the adiabatic approximation. Fig. 12 compares the numerical $\Gamma_{n}$ and the $\Gamma_{n}$ obtained using the adiabatic approximation.

Perturbative elimination of $\hat{Y}_{n}$ to zeroth order

In the main text immediately after the adiabatic approximation, we perturbatively solve the $\hat{Y}_{n}$ variables in terms of $\hat{X}_{n}$ variables [Eq. (63)]. Here (Fig. 13) we numerically examine the validity of this approximation.

Perturbative solution of $\hat{X}_{n>0}$

The validity of the perturbative solution of $\hat{X}_{n>0}$ [Eqs. (75) and (76)] is examined in Fig. 14.

Finally in Fig. 15 we compare the semianalytical formula Eq. (79) to the numerical value.

Appendix C Matrix element estimates

This section provides estimates of matrix elements and quantities derived from those matrix elements in the limit of $\alpha\gg 1$ and $n_{\text{th}}\ll 1$ .

The norm square matrix elements $|\bra{\psi_{m}^{R}}\hat{a}\ket{\psi_{n}^{R}}|^{2}$ between states in the same well can be estimated via perturbation theory in the limit of $\alpha\gg 1$ because the neighborhoods of the two symmetric potential minima are nearly harmonic. Defining a ladder operator $\hat{b}$ centered at one of the minima, i.e. $\hat{a}=\alpha+\hat{b}$ , we can write the Hamiltonian in Eq. (1) as

\displaystyle\hat{H}=-4\epsilon_{2}\left[\hat{b}^{\dagger}\hat{b}+\frac{1}{4\alpha^{2}}\hat{b}^{\dagger 2}\hat{b}^{2}+\frac{1}{2\alpha}(\hat{b}^{\dagger 2}\hat{b}+\hat{b}^{\dagger}\hat{b}^{2})\right]+\frac{\epsilon_{2}^{2}}{K}.

(107)

Standard state-based perturbation theory provides a crude asymptotic bound on norm square matrix element of $\hat{a}$ as

\displaystyle\left|\bra{\psi_{n}^{R}}\hat{a}\ket{\psi_{m}^{R}}\right|^{2}

\displaystyle=O\left(\alpha^{2-2|m-n|}\right).

(108)

Performing Schrieffer-Wolff perturbation theory up to eleventh order with respect to the small parameter $\alpha^{-1}$ (via the python computer algebra package SymPy meurer_sympy_2017) shows that, for all $m-12\leq n\leq m+4$ , the norm square matrix elements of $\hat{a}$ are bounded by the following estimates,

\displaystyle\left|\bra{\psi_{n}^{R}}\hat{a}\ket{\psi_{m}^{R}}\right|^{2}

\displaystyle=\begin{cases}O\left(\alpha^{2-2|m-n|}\right)&n\leq m\\ O\left(\alpha^{2-6|m-n|}\right)&n>m\end{cases},

(109)

and in particular, $\left|\bra{\psi_{n}^{R}}\hat{a}\ket{\psi_{n}^{R}}\right|^{2}\sim\alpha^{2}$ . We will assume that this estimate continues to apply for those pairs of $m$ and $n$ that we have not verified.

The norm square matrix elements $|\bra{\psi_{m}^{L}}\hat{a}\ket{\psi_{n}^{R}}|^{2}$ between states in different wells are exponentially suppressed in the limit of $\alpha\gg 1$ because the right- and left-well states have nearly disjoint support.

Using the estimates of the matrix elements of $\hat{a}$ , we can assemble estimates for $W_{mn}$ and $V_{mn}$ in the limit of $n_{\text{th}}\ll 1$ and $\alpha\gg 1$ :

$\displaystyle 2W_{nn}$	$\displaystyle\sim 2\kappa\|\bra{\psi_{n}^{R}}\hat{a}\ket{\psi_{n}^{R}}\|^{2}\sim 2\kappa\alpha^{2},$	(110)
$\displaystyle W_{fn}$	$\displaystyle=\begin{cases}O(\kappa\alpha^{2-2\|f-n\|})&f<n\\ O(\kappa\alpha^{2-6\|f-n\|})+O(\kappa n_{\text{th}}\alpha^{2-2\|f-n\|})&f>n\end{cases}$	(111)
$\displaystyle V_{fn}$	$\displaystyle=O\left(h(\alpha^{2})e^{-2\alpha^{2}}\right),$	(112)

where $h(\alpha^{2})$ is a polynomial of $\alpha^{2}$ .

Based on the estimates of $W_{mn}$ and $V_{mn}$ , the measurement rate scales as $\mu_{n}\sim 2\kappa\alpha^{2}$ in the limit of $n_{\text{th}}\ll 1$ and $\alpha\gg 1$ .

Appendix D Rewriting the solution

To elucidate the physics of the staircase contained in the expression of $\Gamma$ , we rewrite Eq. (78) into the following form (see appendix D for steps)

must perform some manipulations. We first notice that in the long time and adiabatic limit of $\Gamma\ll\lambda_{n}$ , Eq. (67) provides a conservation law for all manifolds $n>0$ ,

\displaystyle\left(\frac{\delta_{n}^{2}}{\mu_{n}}+\lambda_{n}\right)\braket{\hat{X}_{n}}=\sum_{i\neq n}(W_{ni}-V_{ni})\braket{\hat{X}_{i}}.

(113)

The equality holds exacts for our solutions of $\braket{\hat{X}_{n}}$ because the same adiabatic approximation was employed to obtain this equality and the solutions. Although the equation deals with population differences $\braket{\hat{X}_{n}}$ , the physical interpretation of this equation becomes clear if we temporarily express $\braket{\hat{X}_{n}}$ in terms of absolute populations. In the long time limit, the populations of the states $\ket{\psi_{n}^{R/L}}$ are close to a quasi-equilibrium distribution,

	$\displaystyle p_{n,R}(t)$	$\displaystyle=c_{n}\ p_{R}(t)$		(114)
	$\displaystyle p_{n,L}(t)$	$\displaystyle=c_{n}\ p_{L}(t),$		(115)

where $p_{R}(t)$ , $p_{L}(t)$ are the total populations in the right- and left-well, and $c_{n}$ is a constant that determines how population is distributed within a well. Then

\displaystyle\braket{\hat{X}_{n}}=p_{n,R}(t)-p_{n,L}(t)=(p_{R}(t)-p_{L}(t))c_{n}.

(116)

Substituting this value into the conservation law, dividing by the factor $(p_{R}(t)-p_{L}(t))/p_{R}(t)$ and then identifying $c_{n}p_{R}(t)$ as the population $p_{n,R}(t)$ gives

\displaystyle\left(\frac{\delta_{n}^{2}}{\mu_{n}}+\lambda_{n}\right)p_{n,R}(t)=\sum_{i\neq n}(W_{ni}-V_{ni})p_{i,R}(t).

(117)

As a reminder, $\frac{\delta_{n}^{2}}{\mu_{n}}$ is the effective tunneling rate due to the quantum Zeno effect, and $\lambda_{n}$ is the total decay rate of population in $\ket{\psi_{n}^{R}}$ except due to tunneling. Therefore the left-hand side of the equality describes the total probability current exiting the state $\ket{\psi_{n}^{R}}$ . The right hand side of the equality describes the total probability current entering the state $\ket{\psi_{n}^{R}}$ (See footnote ²²2The minus $V_{ni}$ is because this part of the probability current entering the state $\ket{\psi_{n}^{R}}$ , together with the probability current entering the state $\ket{\psi_{n}^{L}}$ , forms a equal-weight mixture between the left and right wells and do not contribute to the inter-well population difference for a discussion about the presence of $V_{ni}$ ). Thus the physical interpretation of the original equality Eq. (113) is conservation of total probability.

The conservation law [Eq. (113)] will help us isolate the terms that change quickly with $\alpha$ from those that change slowly with $\alpha$ in the solution of the spontaneous switching rate (78). This is because the incoming probability current on the right-hand side of Eq. (113) is dominated by contributions from manifolds that are well below the energy barrier and therefore changes slowly with $\alpha^{2}$ . If we can factor out such a term from Eq. (65), then we can identify the remaining factor as being responsible for the staircase. To perform the factorization, we return to the expression Eq. (65), and adjust the normalization of $\braket{\hat{X}_{n}}$ such that $\sum_{m}\langle\hat{X}_{m}\rangle=1$ ,

\displaystyle\Gamma=\sum_{n}\left[\frac{\delta_{n}^{2}}{\mu_{n}}+2\sum_{f}V_{fn}\right]\braket{\hat{X}_{n}}.

(118)

Then, for each term under the sum, we divide by the left hand side of the conservation law Eq. (113) and multiply by the right hand side of the same, yielding

\displaystyle\Gamma=\sum_{n}\frac{\frac{\delta_{n}^{2}}{\mu_{n}}+2\sum_{f}V_{fn}}{\frac{\delta_{n}^{2}}{\mu_{n}}+\lambda_{n}}\left(\sum_{i\neq n}(W_{ni}-V_{ni})\braket{\hat{X}_{i}}\right),

(119)

with the understanding that $\braket{\hat{X}_{i}}=K_{n}/(\sum_{m\geq 0}K_{m})$ as before.

Appendix E WKB approximation of the tunnel splitting

The tunnel splittings $\delta_{n}$ of the Hamiltonian in Eq. (1) has been estimated via overlap of shifted Fock states puri_stabilized_2019, empirical fits putterman_stabilizing_2022, and the WKB approximation marthaler_switching_2006; marthaler_quantum_2007; venkatraman_driven_2024. The existing estimates in these previous works only apply to states that are sufficiently inside the well. We are instead interested in the tunnel splittings of the states lying near the top of the energy barrier because of their overwhelming contribution to the spontaneous switching rate. In this section we use the WKB approximation to find the tunnel splittings of these states, taking particular care to use a good approximation of the wavefunction in the region under the central barrier. We start from the Hamiltonian of Eq. (1) but defined with an additional detuning term $\Delta\hat{a}^{\dagger}\hat{a}$ that we will later set to zero,

\displaystyle\hat{H}=\Delta\hat{a}^{\dagger}\hat{a}-K\hat{a}^{\dagger 2}\hat{a}^{2}-\epsilon_{2}(\hat{a}^{\dagger 2}+\hat{a}^{2}),

(120)

where $[\hat{a},\hat{a}^{\dagger}]=1$ . To connect with the WKB calculation performed in marthaler_switching_2006; marthaler_quantum_2007; peano_sharp_2012, we non-dimensionalize the problem by defining rescaled Hamiltonian $\hat{g}$ and rescaled phase space coordinates $\hat{Q}$ and $\hat{P}$ as follows,

	$\displaystyle\hat{H}$	$\displaystyle=-\frac{4\epsilon_{2}^{2}}{K}\hat{g}-\frac{2\Delta+3K}{4},$		(121)
	$\displaystyle\hat{a}$	$\displaystyle=\frac{\hat{Q}+i\hat{P}}{\sqrt{2\lambda}},\quad\text{with }[\hat{Q},\hat{P}]=i\lambda.$		(122)

The scale conversion between eigenvalues of $\hat{H}$ and $\hat{g}$ is

\displaystyle\delta H=-\frac{4\epsilon_{2}^{2}}{K}\,\delta g.

(123)

and the reference energy when $g=0$ is $E_{\text{ref}}=-(2\Delta+3K)/4$ . Here the conversion factor $4\epsilon_{2}^{2}/K$ is four times the height of the energy barrier (or four times the depth of the double well), and the free parameter $\lambda$ defines the scale of the new phase space coordinates $\hat{Q}$ and $\hat{P}$ as well as function as an effective Planck constant. Choosing $\lambda=K/(2\epsilon_{2})$ absorbs most of the $\epsilon_{2}$ dependence of the Hamiltonian into the effective Planck constant and results in the following expression for the rescaled Hamiltonian $\hat{g}$ that was used in marthaler_switching_2006; marthaler_quantum_2007; peano_sharp_2012,

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

(124)

where $\mu=(\Delta+2K)/(2\epsilon_{2})$ is a rescaled detuning parameter. In the limit of zero detuning and large two-photon drive amplitude, we have $\Delta=0$ , $\epsilon_{2}=K\alpha^{2}$ , and $\alpha^{2}\gg 1$ . Consequently the effective Plank constant scales as $\lambda=1/(2\alpha^{2})\ll 1$ , justifying the WKB approximation, and the rescaled detuning parameter scales as $\mu=1/\alpha^{2}\ll 1$ , only slightly perturbing the functional form of $\hat{g}$ .

E.1 Characterizing classical orbits in phase space

We characterize the classical orbits in phase space in preparation of a WKB calculation. Removing hats from Eq. (124) produces the following classical version of the rescaled Hamiltonian,

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

(125)

The energy landscape in phase space has a double well structure for $|\mu|<1$ . In particular, for each energy $g$ within the range $-(1+\mu)^{2}/4<g<0$ , there are two orbits bound to the double well, one in the $Q>0$ half-plane and one in the $Q<0$ half-plane. For $g>0$ , the two orbits merge into a single unbound orbit.

The WKB solution expresses the wavefuntion in terms of the momentum and the velocity along the orbit, and uses the action enclosed by the orbit to determine the quantization condition. We define these quantities here. The momentum along the orbit is the function $P(Q,g;\mu)$ defined as the positive real solution of Eq. (125). The velocity along the orbit is obtained from the Hamilton’s equations

\displaystyle\dot{Q}(Q,g;\mu)=\left.\frac{\partial g}{\partial P}\right|_{P=P(Q,g;\mu)}.

(126)

Finally the action enclosed by the orbit (restricted to the $Q>0$ half-plane) is defined as

	$\displaystyle S(g,\mu)$	$\displaystyle=\oint_{Q>0}P\,\mathrm{d}Q$		(127)
		$\displaystyle=2\int_{Q_{\min}}^{Q_{\max}}P(Q,g;\mu)\,\mathrm{d}Q,$		(128)

where $Q_{\min}=Q_{\min}(g,\mu)$ and $Q_{\max}=Q_{\max}(g,\mu)$ denote the minimum and maximum extent of the classical orbit in the $Q>0$ half-plane (thus setting $Q_{\min}=0$ for $g>0$ ).

To slightly simplify the integrand that defines the enclosed action, we apply the following canonical transformation, akin to polar coordinates, to Eq. (125),

	$\displaystyle Q$	$\displaystyle=\sqrt{2J}\cos(\theta)$		(129)
	$\displaystyle P$	$\displaystyle=-\sqrt{2J}\sin(\theta),$		(130)

and arrive at the following equation for the classical orbit in terms of the new coordinates $\theta,J$ ,

\displaystyle g=J^{2}-J(\cos 2\theta+\mu).

(131)

Solving for $J$ in terms of $\theta$ gives

\displaystyle J(\theta,g;\mu)=\frac{1}{2}\left(\cos 2\theta+\mu\pm\sqrt{(\cos 2\theta+\mu)^{2}+4g}\right).

(132)

Using this solution, we explicitly write down the enclosed action in terms of the $\theta,J$ variables as follows,

	$\displaystyle\quad S(g,\mu)$
	$\displaystyle=\oint_{\|\theta\|<\pi/2}J(\theta,g;\mu)\mathrm{d}\theta$		(133)
	$\displaystyle=\begin{cases}2\int_{0}^{\frac{1}{2}\arccos(\sqrt{-4g}-\mu)}\sqrt{(\cos 2\theta+\mu)^{2}+4g}\ \mathrm{d}\theta&g\leq 0\\ 2\int_{0}^{\pi/2}\frac{\cos 2\theta+\mu+\sqrt{(\cos 2\theta+\mu)^{2}+4g}}{2}\ \mathrm{d}\theta&g>0\end{cases}.$		(134)

We are interested in the enclosed action for orbits near the barrier top $g\to 0$ . A first approximation for the enclosed action is the value at $g=0$ , corresponding to the phase space area occupied by each well,

\displaystyle S(0,\mu)=\frac{\pi\mu}{2}+\sqrt{1-\mu^{2}}+\mu\arcsin(\mu).

(135)

Away from $g=0$ , it is useful to find an asymptotic expression of $S(g,\mu)$ in the limit of $g\to 0$ . After a tedious calculation breaking the integral $S(g<0,\mu)$ into several pieces, computing the asymptotic behavior of each piece in the $g\to 0$ limit and summing them back up, we find that $S(g,\mu)$ has the following asymptotics in the $g\to 0$ limit,

	$\displaystyle\quad S(g,\mu)$
	$\displaystyle=S(0,\mu)+\frac{1}{\sqrt{1-\mu^{2}}}\left[1-\ln\left(\frac{\|g\|}{4(1-\mu^{2})^{2}}\right)\right]g$
	$\displaystyle\quad+\frac{3+14\mu^{2}+(2+4\mu^{2})\ln\left(\frac{\|g\|}{4(1-\mu^{2})^{2}}\right)}{4(1-\mu^{2})^{5/2}}g^{2}$
	$\displaystyle\quad+O(\|g\|^{3}\ln\|g\|),$		(136)

where the use of the absolute value function allows for $g>0$ as well as $g<0$ .

For problems involving tunneling under the barrier, we must additionally consider imaginary solutions $P(Q,g;\mu)$ , which characterizes the exponential tail of the WKB wavefunction, as well as real solutions. Therefore we let $P(Q,g;\mu)=i\Pi(Q,g;\mu)$ . Then Eq. (125) gives the following implicit equation for the imaginary solutions,

\displaystyle g=\frac{1}{4}(Q^{2}-\Pi^{2})^{2}-\frac{1}{2}(1-\mu)\Pi^{2}-\frac{1}{2}(1+\mu)Q^{2}.

(137)

The WKB tunnel splitting is related to the tunnel action integral, which is defined for $-(1-\mu)^{2}/4<g\leq 0$ (a more restrictive range than simply requiring bound states) as the integral of $\Pi(Q,g;\mu)$ in the classically forbidden region,

\displaystyle I(g,\mu)=\int_{-Q_{\min}}^{Q_{\min}}\Pi(Q,g;\mu)\,\mathrm{d}Q.

(138)

The integrand can again be simplified by performing the following canonical transformation on Eq. (137),

	$\displaystyle Q$	$\displaystyle=\sqrt{2F}\cos(\theta)$		(139)
	$\displaystyle\Pi$	$\displaystyle=-\sqrt{2F}\sin(\theta),$		(140)

which produces the equation

\displaystyle g=F^{2}\cos^{2}2\theta-F(1+\mu\cos 2\theta).

(141)

The relevant branch of solution in this case is

\displaystyle F(\theta,g;\mu)=\frac{1+\mu\cos\theta-\sqrt{(1+\mu\cos 2\theta)^{2}+4g\cos^{2}2\theta}}{2\cos^{2}2\theta}.

(142)

We thus write the tunnel action explicitly as

	$\displaystyle\quad I(g,\mu)$
	$\displaystyle=\int_{0}^{\pi}F(\theta,g;\mu)\,\mathrm{d}\theta$		(143)
	$\displaystyle=\int_{0}^{\pi}\frac{1+\mu\cos\theta-\sqrt{(1+\mu\cos 2\theta)^{2}+4g\cos^{2}2\theta}}{2\cos^{2}2\theta}\mathrm{d}\theta.$		(144)

We remark that the above integral is in fact defined for $g>0$ as well as $g\leq 0$ and analytic at $g=0$ . With the help of Mathematica, the asymptotics of $I(g,\mu)$ near $g\to 0$ is found to be

	$\displaystyle I(g,\mu)$	$\displaystyle=-\frac{\pi}{\sqrt{1-\mu^{2}}}g+\frac{\pi(1+2\mu^{2})}{2(1-\mu^{2})^{5/2}}g^{2}$
		$\displaystyle\quad+\frac{\pi(-3-8\mu^{2}(3+\mu^{2}))}{4(1-\mu^{2})^{9/2}}g^{3}+O(g^{4}).$		(145)

We also need its asymptotics near $g\to-\frac{1}{4}(1-\mu)^{2}$ . This is because for the $n$ th excited state in the well, we have that $E_{n}=\frac{\epsilon_{2}^{2}}{K}-4n\epsilon_{2}+O(K)$ , so

$\displaystyle g_{n}$	$\displaystyle=-\frac{K}{\epsilon_{2}^{2}}\left(E_{n}+\frac{3K}{4}\right)$	(146)
	$\displaystyle=-\frac{1}{4}+n\frac{1}{\alpha^{2}}+O\left(\frac{1}{\alpha^{4}}\right)$	(147)
	$\displaystyle=-\frac{1}{4}(1-\mu)^{2}+\left(n-\frac{1}{2}\right)\mu+O(\mu^{2})$	(148)

Let $g=-\frac{1}{4}(1-\mu)^{2}+\tilde{g}\mu$ where $\tilde{g}=n-1/2+O(\mu)$ . Then, in the limit $\mu\to 0$ , and the asymptotics for $I(g,\mu)$ (as a function of $\mu$ ) can be found to be the following after a grueling calculation

\displaystyle I=1-\delta I

(149)

where

$\displaystyle\delta I$	$\displaystyle=-\frac{1}{2}\left\{\mu(1+\tilde{g})\ln[\mu(1+\tilde{g})]+\mu\tilde{g}\ln(\mu\tilde{g})\right\}$
	$\displaystyle\quad+\left[\frac{1}{2}+\ln(2)\right](1+2\tilde{g})\mu$
	$\displaystyle\quad-\frac{\tilde{g}(\tilde{g}+1)}{2}\mu^{2}\ln\left[\mu\sqrt{(4+4\tilde{g})4\tilde{g}}\right]+O(\mu^{2})$	(150)

E.2 Solving for approximate wavefunctions

In the classically allowed region within $Q>0$ , standard WKB approximation gives the position basis wavefunction $\psi(Q)$ in terms of the classical momentum $P(Q,g)$ and the classical velocity $\dot{Q}(Q,g)$ as follows,

\displaystyle\psi(Q)=\frac{C}{\sqrt{\dot{Q}(Q,g;\mu)}}\cos\left(\frac{1}{\lambda}\int_{Q_{\max}}^{Q}P(Q^{\prime},g;\mu)dQ^{\prime}+\frac{\pi}{4}\right).

(151)

This wavefunction matches the boundary condition at $Q=Q_{\max}$ .

Now we obtain a wavefunction that matches the boundary condition at $Q=Q_{\min}\approx 0$ . To do so, we use a quadratic approximation for the Hamiltonian $\hat{g}$ ,

\displaystyle\hat{g}\approx\frac{1}{2}(1-\mu)\hat{P}^{2}-\frac{1}{2}(1+\mu)\hat{Q}^{2}.

(152)

This approximation captures the energy saddle point in phase space. It is justified because for a state whose energy is close to the barrier top, both the position and the momentum are small in the region near $Q=0$ . From Eq. (152), we can identify an effective mass is $(1-\mu)^{-1}$ and an effective spring constant $(1+\mu)$ . Together they give rise to the effective frequency $\omega_{b}$ , a length scale $Q_{b}$ and an energy scale $g_{b}$ ,

	$\displaystyle\omega_{b}$	$\displaystyle=\sqrt{\frac{1+\mu}{(1-\mu)^{-1}}}=\sqrt{1-\mu^{2}},$		(153)
	$\displaystyle g_{b}$	$\displaystyle=\lambda\omega_{b}=\lambda\sqrt{1-\mu^{2}}.$		(154)

Since $g_{n}$ sets a natural energy scale for motion near the barrier top, we define a new energy variable $a$ by taking the ratio of $g$ and $g_{n}$ ,

\displaystyle a=-\frac{g}{g_{n}}=-\frac{g}{\lambda\sqrt{1-\mu^{2}}}.

(155)

The time-independent Shrödinger equation corresponding to Eq. (152) is

\displaystyle-\frac{\lambda^{2}}{2}(1-\mu)\frac{\partial^{2}\psi}{\partial Q^{2}}-\frac{1}{2}(1+\mu)Q^{2}\,\psi=g\,\psi

(156)

where $\psi(Q)$ is the position basis wavefunction. Changing into the new energy variable $a$ and a new position variable as

\displaystyle Q

\displaystyle=\left(\frac{\lambda^{2}(1-\mu)}{4(1+\mu)}\right)^{\frac{1}{4}}\xi

(157)

reduces the differential equation for $\psi$ to the following standard form,

\displaystyle\frac{\partial^{2}\psi}{\partial\xi^{2}}+\left(\frac{1}{4}\xi^{2}-a\right)\psi=0,

(158)

for which the general solutions are

\displaystyle\psi=c_{1}W(a,\xi)+c_{2}W(a,-\xi)

(159)

where $W$ is a parabolic cylinder function (see section 12.14 of the Digital Library of Mathematical Functions, DLMF olver_nist_2024.) The wavefunction can be made to have even or odd parity like so,

\displaystyle\psi_{\pm}\propto W(a,\xi)\pm W(a,-\xi).

(160)

E.3 Matching the two wavefunctions to produce a quantization condition

To produce a quantization condition, we need to match the parabolic cylinder function solution in Eq. (160) to the WKB wavefunction in Eq. (151). We first apply the following asymptotics of the parabolic cylinder function $W$ for large $\xi$ (see DLMF Eq. 12.14.17 and 12.14.18),

	$\displaystyle W(a,\xi)\sim\sqrt{\frac{2k}{\xi}}\cos\omega,$		(161)
	$\displaystyle W(a,-\xi)\sim\sqrt{\frac{2}{k\xi}}\sin\omega,$		(162)

where

$\displaystyle\omega$	$\displaystyle=\frac{1}{4}\xi^{2}-a\ln\xi+\frac{1}{4}\pi$
	$\displaystyle\quad+\frac{1}{2}\arg\Gamma\left(\frac{1}{2}+ia\right)+O(\xi^{-2}),$	(163)
$\displaystyle\frac{1}{k}$	$\displaystyle=\sqrt{1+e^{2\pi a}}+e^{\pi a}.$	(164)

Substituting the asymptotics into the parabolic cylinder function solution Eq. (160) gives the following wavefunction in the large $\xi$ limit,

$\displaystyle\psi_{\pm}$	$\displaystyle=\frac{C}{\sqrt{\xi}}\cos\left[\frac{1}{4}\xi^{2}-a\ln\xi+\frac{1}{4}\pi\right.$
	$\displaystyle\quad\quad+\frac{1}{2}\arg\Gamma\left(\frac{1}{2}+ia\right)$
	$\displaystyle\quad\quad\left.\mp\arctan(\sqrt{1+e^{2\pi a}}+e^{\pi a})+O(\xi^{-2})\right].$	(165)

Next we rewrite the amplitude and the phase of the wavefunction above in terms of the classical action and velocity. Using the classical Hamiltonian obtained by removing the hats from Eq. (152),

\displaystyle g\approx\frac{1}{2}(1-\mu)P^{2}-\frac{1}{2}(1+\mu)Q^{2},

(166)

we can find the following expressions for the momentum and velocity, valid near $Q=0$ ,

$\displaystyle P(Q,g;\mu)$	$\displaystyle=\sqrt{\frac{2}{1-\mu}\left[g+\frac{1}{2}(1+\mu)Q^{2}\right]},$	(167)
$\displaystyle\dot{Q}(Q,g;\mu)$	$\displaystyle=\frac{\partial g}{\partial P}=(1-\mu)P$
	$\displaystyle=\sqrt{2(1-\mu)\left[g+\frac{1}{2}(1+\mu)Q^{2}\right]}.$	(168)

If we change into the $(\xi,a)$ variables, we obtain the following expression for the momentum,

	$\displaystyle P(Q,g;\mu)$	$\displaystyle=\lambda\left(\frac{4(1+\mu)}{\lambda^{2}(1-\mu)}\right)^{1/4}\sqrt{\frac{1}{4}\xi^{2}-a},$		(169)
	$\displaystyle\dot{Q}(Q,g;\mu)$	$\displaystyle\propto\sqrt{\frac{1}{4}\xi^{2}-a}.$		(170)

Substituting the momentum into the following (normalized) action integral, integrating, and finding the asymptotics at $\xi\to\infty$ produces

	$\displaystyle\frac{1}{\lambda}\int_{Q_{\min}}^{Q}P(Q^{\prime},g;\mu)\,\mathrm{d}Q^{\prime}$
	$\displaystyle=\int_{\xi_{\min}}^{\xi}\sqrt{\frac{1}{4}(\xi^{\prime})^{2}-a}\,\mathrm{d}\xi^{\prime}$		(171)
	$\displaystyle=\frac{1}{4}\xi^{2}-a\ln\xi+\frac{1}{2}a\ln\|a\|-\frac{1}{2}a+O(\xi^{-2}),$		(172)

where $Q_{\min}$ and $\xi_{\min}$ denote the minimum position reached by the classical orbit if we restrict to the $Q>0$ half plane ( $Q_{\min}=\xi_{\min}=0$ for $g>0$ ). Taking the $\xi\to\infty$ asymptotics of the velocity Eq. (170) gives

\displaystyle\dot{Q}

\displaystyle\propto\sqrt{\frac{1}{4}\xi^{2}-a}\sim\frac{1}{2}\xi.

(173)

Using the asymptotics of the action and the velocity above, we can rewrite the solution for the wavefunction in Eq. (165) as follows,

$\displaystyle\psi_{\pm}$	$\displaystyle=\frac{C^{\prime}}{\sqrt{\dot{Q}(Q,g;\mu)}}\cos\left[\frac{1}{\lambda}\int_{Q_{\min}}^{Q}P(Q^{\prime},g;\mu)\ \mathrm{d}Q^{\prime}\right.$
	$\displaystyle\quad\quad-\frac{1}{2}a\ln\|a\|+\frac{1}{2}a+\frac{1}{4}\pi$
	$\displaystyle\quad\quad+\frac{1}{2}\arg\Gamma\left(\frac{1}{2}+ia\right)$
	$\displaystyle\quad\quad\left.\mp\arctan(\sqrt{1+e^{2\pi a}}+e^{\pi a})+O(\xi^{-2})\right].$	(174)

This is a solution of the wavefunction that satisfies the boundary condition near $Q=0$ .

We now have two approximate wavefunctions Eqs. (151) and (174), each satisfying the boundary condition at one position extremum of the classical orbit. Now we match the two to obtain a quantization condition for the energy $g$ . Since the the arguments of the cosine function must differ by an integer multiple of $\pi$ for the solutions to match, we have the following quantization condition on the action integral (omitting higher order terms),

	$\displaystyle\ \frac{1}{\lambda}\int_{Q_{\min}}^{Q_{\max}}P(Q^{\prime},g;\mu)\ \mathrm{d}Q^{\prime}$
	$\displaystyle=\pi n+\frac{1}{2}a\ln\|a\|-\frac{1}{2}a-\frac{1}{2}\arg\Gamma\left(\frac{1}{2}+ia\right)$
	$\displaystyle\quad\quad\pm\arctan(\sqrt{1+e^{2\pi a}}+e^{\pi a}).$		(175)

This condition can be rewritten as

	$\displaystyle\frac{1}{\lambda}S(g,\mu)$	$\displaystyle=2\pi n+a\ln\|a\|-a-\mathrm{Im}\left[\ln\Gamma\left(\frac{1}{2}+ia\right)\right]$
		$\displaystyle\quad\pm\left(\pi-\arctan(e^{-\pi a})\right).$		(176)

Adding the term $(\pi\mp\pi)$ to the right-hand side allows us to shift the index $n$ by 1 for all odd parity states without changing physical content of the quantization condition,

	$\displaystyle\frac{1}{\lambda}S(g,\mu)$	$\displaystyle=2\pi\left(n+\frac{1}{2}\right)+a\ln\|a\|-a$
		$\displaystyle\quad-\mathrm{Im}\left[\ln\Gamma\left(\frac{1}{2}+ia\right)\right]\mp\arctan(e^{-\pi a}).$		(177)

Upper sign for even states, lower sign for odd states. In the limit of $a\gg 1$ (states deep in the well), the quantization condition reduces to the standard one,

\displaystyle\frac{1}{\lambda}S(g,\mu)

\displaystyle=2\pi\left(n+\frac{1}{2}\right)\mp e^{-\pi a},

(178)

where $n=0$ denotes the ground state in each parity sector. The simple form of the tunnel exponent $\pi a$ is a side effect of using a quadratic approximation for the energy barrier. We can fix this by simply go back to (177) and replace the tunnel exponent by the actual integral

\displaystyle\pi a\to\frac{1}{\lambda}I(g,\mu).

(179)

Therefore, the final quantization condition for even and odd parity states are

	$\displaystyle\ \frac{1}{\lambda}S(g,\mu)$
	$\displaystyle=2\pi\left(n+\frac{1}{2}\right)+a\ln\|a\|-a-\mathrm{Im}\left[\ln\Gamma\left(\frac{1}{2}+ia\right)\right]$
	$\displaystyle\quad\quad\mp\arctan\left\{\exp\left[-\frac{1}{\lambda}I(g,\mu)\right]\right\},$		(180)

where $S(g,\mu)$ is defined in Eq. (134), $I(g,\mu)$ in Eq. (144), and $a$ in Eq. (155). This quantization condition allows a smooth interpolation between the quantization condition for states inside and outside the double well. For states deep in the well and away from oscillatory behavior in the classically forbidden region the quantization condition is approximately

\displaystyle\ \frac{1}{\lambda}S(g,\mu)=2\pi\left(n+\frac{1}{2}\right)\mp\exp\left[-\frac{1}{\lambda}I(g,\mu)\right],

(181)

which gives a tunnel splitting in terms of the energy gap between manifolds to be

\displaystyle\frac{\delta}{E_{\text{gap}}}=\frac{1}{\pi}e^{-\frac{I}{\lambda}}

(182)

[Plot showing validity of this quantization condition].

E.4 Solving the quantization condition for energy

We are interested in the energies of states near the top of the barrier $g=0$ ( $a=0$ ). Here it is more convenient to work with the new energy variable $a$ rather than $g$ . From the asymptotics of $S(g,\mu)$ and $I(g,\mu)$ in Eq. (136), Eq. (145), and the definition of $a$ in Eq. (155), we can rewrite the asymptotics into the following form,

$\displaystyle\frac{1}{\lambda}S(g,\mu)$	$\displaystyle=\frac{1}{\lambda}S(0,\mu)+\left[\ln\left(\frac{\lambda\|a\|}{4(1-\mu^{2})^{3/2}}\right)-1\right]a$
	$\displaystyle\quad-\frac{\lambda}{4(1-\mu^{2})^{3/2}}\biggl{[}-3-14\mu^{2}$
	$\displaystyle\quad\quad-(2+4\mu^{2})\ln\left(\frac{\lambda\|a\|}{4(1-\mu^{2})^{3/2}}\right)\biggr{]}a^{2}$
	$\displaystyle\quad+O(\|a\|^{3}\ln\|a\|),$	(183)
$\displaystyle\frac{1}{\lambda}I(g,\mu)$	$\displaystyle=\pi a\left(1+\frac{1+2\mu^{2}}{2(1-\mu^{2})^{3}/2}\lambda a\right.$
	$\displaystyle\quad\left.+\frac{3+8\mu^{2}(3+\mu^{2})}{4(1-\mu^{2})^{3}}\lambda^{2}a^{2}+O(a^{3})\right).$	(184)

Plug the asymptotics into the quantization condition Eq. (180) and cancelling common terms gives an equation for $a$

	$\displaystyle\frac{1}{\lambda}S(0,\mu)+a\ln\left(\frac{\lambda}{4(1-\mu^{2})^{3/2}}\right)$
	$\displaystyle\quad-\frac{\lambda}{4(1-\mu^{2})^{3/2}}\biggl{[}-3-14\mu^{2}$
	$\displaystyle\quad\quad-(2+4\mu^{2})\ln\left(\frac{\lambda\|a\|}{4(1-\mu^{2})^{3/2}}\right)\biggr{]}a^{2}$
	$\displaystyle=2\pi\left(n+\frac{1}{2}\right)-\mathrm{Im}\left[\ln\Gamma\left(\frac{1}{2}+ia\right)\right]$
	$\displaystyle\quad\quad\mp\arctan\left(e^{-I(g,\mu)/\lambda}\right)$		(185)

Let’s only keep terms of order $a$ . Then.

	$\displaystyle\frac{1}{\lambda}S(0,\mu)+a\ln\left(\frac{\lambda}{4(1-\mu^{2})^{3/2}}\right)$
	$\displaystyle=2\pi\left(n+\frac{1}{2}\right)-\mathrm{Im}\left[\ln\Gamma\left(\frac{1}{2}+ia\right)\right]$
	$\displaystyle\quad\quad\mp\arctan\left(e^{-\pi a}\right)$		(186)

Lets define $2\pi[n_{0}(\lambda,\mu)+1/2]=\frac{1}{\lambda}S(0,\mu)$ . Then we need to solve

	$\displaystyle a\ln\left(\frac{\lambda}{4(1-\mu^{2})^{3/2}}\right)$
	$\displaystyle=2\pi\left(n-n_{0}\right)-\mathrm{Im}\left[\ln\Gamma\left(\frac{1}{2}+ia\right)\right]\mp\arctan\left(e^{-\pi a}\right)$		(187)
	$\displaystyle=2\pi\left(n-n_{0}\right)+h_{\pm}(a)$		(188)

where

\displaystyle h_{\pm}(a)=-\mathrm{Im}\left[\ln\Gamma\left(\frac{1}{2}+ia\right)\right]\mp\arctan\left(e^{-\pi a}\right).

(189)

which satisfies $h_{\pm}(a=0)=\mp\frac{\pi}{4}$ . We may use this condition to solve for $a_{\pm}(n,\mu)$ , the (rescaled) energy as a function of the principle quantum number $n$ , the parity, and the rescaled detuning parameter $\mu$ . Zeroth order approximation is

	$\displaystyle a_{\pm}^{(0)}\ln\left(\frac{\lambda}{4(1-\mu^{2})^{3/2}}\right)$	$\displaystyle=2\pi\left(n-n_{0}\right)\mp\frac{\pi}{4}$		(190)
	$\displaystyle a_{\pm}^{(0)}$	$\displaystyle=\frac{2\pi\left(n-n_{0}\right)\mp\frac{\pi}{4}}{\ln\left(\frac{\lambda}{4(1-\mu^{2})^{3/2}}\right)}$		(191)

Plug back in to $h_{\pm}(a)$ to get first order solution

	$\displaystyle a_{\pm}^{(1)}$	$\displaystyle\approx\frac{2\pi\left(n-n_{0}\right)+h_{\pm}(a_{\pm}^{(0)})}{\ln\left(\frac{\lambda}{4(1-\mu^{2})^{3/2}}\right)}$		(192)
		$\displaystyle=\frac{2\pi\left(n-n_{0}\right)-\mathrm{Im}\left[\ln\Gamma\left(\frac{1}{2}+ia_{\pm}^{(0)}\right)\right]\mp\arctan\left(e^{-\pi a_{\pm}^{(0)}}\right)}{\ln\left(\frac{\lambda}{4(1-\mu^{2})^{3/2}}\right)}$		(193)

To find the tunnel splitting, we use the formula

\displaystyle\Delta a=\left.\frac{\Delta n}{-\mathrm{d}n_{-}(a)/\mathrm{d}a}\right|_{a_{-}(n)}

(194)

where

\displaystyle\Delta n=n_{+}(a)-n_{-}(a)

(195)

and

\displaystyle 0=2\pi(n_{+}-n_{-})-2\arctan(e^{-\pi a})

(196)

therefore $\Delta n=\frac{1}{\pi}\arctan(e^{-\pi a})$ and $\mathrm{d}n_{-}(a)/\mathrm{d}a$ is defined by taking the derivative of the condition

	$\displaystyle\ln\left(\frac{\lambda}{4(1-\mu^{2})^{3/2}}\right)$
	$\displaystyle=2\pi\frac{\mathrm{d}n_{\pm}(a)}{\mathrm{d}a}-\mathrm{Re}\left[\Psi\left(\frac{1}{2}+ia\right)\right]\mp\left[-\frac{\pi}{2}\mathrm{sech}(\pi a)\right]$		(197)

$\Psi(z)$ is the digamma function. So

\displaystyle\frac{\mathrm{d}n_{-}(a)}{\mathrm{d}a}=\frac{1}{2\pi}\left\{\ln\left(\frac{\lambda}{4(1-\mu^{2})^{3/2}}\right)+\mathrm{Re}\left[\Psi\left(\frac{1}{2}+ia\right)\right]+\frac{\pi}{2}\mathrm{sech}(\pi a)\right\}

(198)

	$\displaystyle\Delta a$	$\displaystyle=\frac{\frac{1}{\pi}\arctan(e^{-\pi a})}{-\frac{1}{2\pi}\left\{\ln\left(\frac{\lambda}{4(1-\mu^{2})^{3/2}}\right)+\mathrm{Re}\left[\Psi\left(\frac{1}{2}+ia\right)\right]+\frac{\pi}{2}\mathrm{sech}(\pi a)\right\}}$		(199)
		$\displaystyle=\frac{2\arctan(e^{-\pi a})}{\ln\left(\frac{4(1-\mu^{2})^{3/2}}{\lambda}\right)-\mathrm{Re}\left[\Psi\left(\frac{1}{2}+ia\right)\right]-\frac{\pi}{2}\mathrm{sech}(\pi a)}$		(200)

Convert back to $\Delta g$ then $\Delta E$ to get splitting $\delta_{n}$ ! Only works when $\delta_{n}$ isn’t too small.

Appendix F Small tunnel splitting limit

If $\delta_{n}$ is really small, use

\displaystyle\frac{\delta_{n}}{E_{\text{gap, n}}}=\frac{1}{\pi}e^{-I/\lambda}

(201)

Plug in asymptotics for $I$ and $E_{\text{gap, n}}$ near $g\to-\frac{1}{4}(1-\mu)^{2}$ . The result is

\displaystyle\delta_{n}=\frac{1}{\pi}4K\alpha^{2}\left[\frac{e^{1+2\tilde{g}}}{(1+\tilde{g})^{1+\tilde{g}}\tilde{g}^{\tilde{g}}}\right]4^{1+2g}(\alpha^{2})^{1+2\tilde{g}}e^{-2\alpha^{2}}

(202)

Set $\tilde{g}=n-1/2$ :

\displaystyle\delta_{n}=\frac{4^{2n+1}e^{2n}}{\pi\left(n+\frac{1}{2}\right)^{n+1/2}\left(n-\frac{1}{2}\right)^{n-1/2}}K\alpha^{4n+2}e^{-2\alpha^{2}}

(203)

The pre-factor is $115.9$ for $n=1$ and $980.3$ for $n=2$ .

	$\displaystyle\quad S(g,\mu)$
	$\displaystyle=S(0,\mu)+\frac{1}{\sqrt{1-\mu^{2}}}\left[1-\ln\left(\frac{\|g\|}{4(1-\mu^{2})^{2}}\right)\right]g$
	$\displaystyle\quad+\frac{3+14\mu^{2}+(2+4\mu^{2})\ln\left(\frac{\|g\|}{4(1-\mu^{2})^{2}}\right)}{4(1-\mu^{2})^{5/2}}g^{2}$
	$\displaystyle\quad+O(\|g\|^{3}\ln\|g\|),$		(136)

Unraveling the switching dynamics in the double well of a Kerr-cat qubit

Abstract

I Introduction

II Background

III Solving for the spontaneous switching rate

III.1 Projection of the density operator onto two-level manifolds

III.2 Effective equations of motion of the density operator components

III.3 Perturbative solution of the effective equations of motion

III.3.1 Exponential ansatz in the long time limit

III.3.2 Solving for the quantum Zeno dynamics via the adiabatic elimination of ⟨Yn^⟩\braket{\hat{Y_{n}}}

III.3.3 Solving for the quasi-equilibrium population distribution via the adiabatic approximation on ⟨X^n>0⟩\braket{\hat{X}_{n>0}}

III.3.4 Final solution and summary of the derivation

IV Location of the steps in the staircase

V Application of the semi-analytical formula: activation mechanism within the well

VI Comparison to previous work

VII Conclusions

VIII Acknowledgments

Appendix A Numerical solution of the spontaneous switching rate

Appendix B Validity of approximations used in the derivation

Projection onto two-level manifolds

Adiabatic approximation

Perturbative elimination of Y^n\hat{Y}_{n} to zeroth order

Perturbative solution of X^n>0\hat{X}_{n>0}

Appendix C Matrix element estimates

Appendix D Rewriting the solution

Appendix E WKB approximation of the tunnel splitting

E.1 Characterizing classical orbits in phase space

E.2 Solving for approximate wavefunctions

E.3 Matching the two wavefunctions to produce a quantization condition

E.4 Solving the quantization condition for energy

Appendix F Small tunnel splitting limit

III.3.2 Solving for the quantum Zeno dynamics via the adiabatic elimination of $\braket{\hat{Y_{n}}}$

III.3.3 Solving for the quasi-equilibrium population distribution via the adiabatic approximation on $\braket{\hat{X}_{n>0}}$

Perturbative elimination of $\hat{Y}_{n}$ to zeroth order

Perturbative solution of $\hat{X}_{n>0}$