Schr\"odinger cats and steady states in subharmonic generation with Kerr nonlinearities

TL;DR

This paper investigates the properties of steady states in nonlinear quantum circuits, revealing that true Schr"odinger cats cannot persist under realistic loss conditions, and characterizes the complex mixed states formed.

Contribution

It provides a detailed analysis of the steady states in Kerr nonlinear systems, highlighting the limitations of Schr"odinger cat states and describing the nature of the resulting mixed states.

Findings

True Schr"odinger cats cannot survive with photon loss.

A 'cat-like' steady state exists only at very high nonlinearity.

Steady states are complex mixed states with reduced purity.

Abstract

We discuss general properties of the equilibrium state of parametric down-conversion in superconducting quantum circuits with detunings and Kerr anharmonicities, in the strongly nonlinear regime. By comparing moments of the steady state and those of a Schr\"odinger cat, we show that true Schr\"odinger cats cannot survive in the steady state if there is any single-photon loss. A delta-function 'cat-like' steady-state distribution can be formed, but this only exists in the limit of an extremely large nonlinearity. The steady state is a mixed state, which is more complex than a mixture or linear combination of delta-functions, and whose purity is reduced by driving. We expect this general behaviour to occur in other driven, dissipative quantum subharmonic non-equilibrium open systems.

Equations

$$H_I$$

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

$${\mathcal{L}}_k^{(j)}\left[\rho \right]$$

$$\frac{H_A}{\hslash }={\Delta }_1{a}^{†}a+i\left[\frac{\mathcal{E}}{2}{a}^{†2}-h.c.\right]+\frac{{\chi }_e}{2}{a}^{†2}{a}^2,$$

$$\mathcal{E}=\frac{\kappa }{{\gamma }_2}{\mathcal{E}}_2,{\chi }_e=\chi -\frac{{\Delta }_2}{2}{\left|\frac{\kappa }{{\gamma }_2}\right|}^2.$$

$$\frac{\partial }{\partial t}{\rho }_1$$

$$+\frac{{\gamma }_e^{(2)}}{2}(2a^2{\rho }_1{a}^{†2}-a^{†2}{a}^2{\rho }_1-{\rho }_1{a}^{†2}{a}^2),$$

$${\gamma }_e^{(2)}={\gamma }_1^{(2)}+\frac{{\gamma }_2^{(1)}}{2}{\left|\frac{\kappa }{{\gamma }_2}\right|}^2.$$

$$\frac{H}{\hslash }$$

$$+2\Re ({ϵ}_p{e}^{-i{\omega }_{p}t}+{ϵ}_d{e}^{-i{\omega }_{d}t})(a_r^{†}+a_r),$$

$$\phi$$

$$U=\exp \left[it\left({\omega }_q{n}_q+{\omega }_d{n}_r+\frac{{\omega }_p+{\omega }_d}{2}{n}_s\right)-{\tilde{\xi }}_p{a}_r^{†}+{\tilde{\xi }}_p^{\ast }{a}_r\right],$$

$$\tilde{H}/\hslash$$

$$-\frac{E_J}{\hslash }\left(\cos (\tilde{\phi })+{\tilde{\phi }}^2/2\right),$$

$$\tilde{\phi }$$

$${\tilde{a}}_q$$

$$\tilde{H}\approx H_{\mathrm{shift}}+H_{\mathrm{Kerr}}+H_2,$$

$$\frac{H_{\text{ shift }}}{\hslash }$$

$$+\left({\omega }_r-{\omega }_d-{\delta }_r-2{\chi }_{rr}{\left|{\xi }_p\right|}^2\right)n_r$$

$$+\left({\omega }_s-\frac{{\omega }_p+{\omega }_d}{2}-{\delta }_s-{\chi }_{rs}{\left|{\xi }_p\right|}^2\right)n_s,$$

$$\frac{H_{\text{ Kerr }}}{\hslash }$$

$$-{\chi }_{qs}{n}_q{n}_s-{\chi }_{rs}{n}_r{n}_s,$$

$$\frac{H_2}{\hslash }$$

$$\frac{H_{sr}}{\hslash }$$

$$+g_2^{\ast }{a}_s^2{a}_r^{†}+g_2{\left(a_s^{†}\right)}^2{a}_r+{ϵ}_d{a}_r^{†}+{ϵ}_d^{\ast }{a}_r$$

$$-{\chi }_{rs}{n}_r{n}_s-\sum _{m=r,s}\frac{{\chi }_{mm}}{2}{a}_m^{†}{a}_m^2,$$

$$H_s={ϵ}_2^{\ast }{a}_s^2+{ϵ}_2{\left(a_s^{†}\right)}^2-\frac{{\chi }_{ss}}{2}{a}_s^{{†}^2}{a}_s^2.$$

$$\frac{d}{dt}{\rho }_s=-i\left[H_s,{\rho }_s\right]+\frac{{\kappa }_2}{2}\mathcal{L}\left[a_s^2\right]{\rho }_s+\frac{{\kappa }_s}{2}\mathcal{L}\left[a_s\right]{\rho }_s,$$

$${\hat{\rho }}_1=\iint d\alpha d{\alpha }^{+}P\left(\alpha ,{\alpha }^{+}\right)\frac{|\alpha ⟩⟨{\alpha }^{+\ast }|}{⟨{\alpha }^{+\ast }|\alpha ⟩},$$

$$\frac{\partial P}{\partial t}=\left\{\frac{\partial }{\partial \alpha }\left[\gamma \alpha -\mathcal{E}\left(\alpha \right){\alpha }^{+}\right]+\frac{1}{2}\frac{{\partial }^2}{\partial {\alpha }^2}\mathcal{E}\left(\alpha \right)+h.c.\right\}P,$$

$$\text{normalcolor}\frac{\partial P\left(\vec{\beta }\right)}{\partial \tau }$$

$$P_1\left(\vec{\beta }\right)=N\exp \left[-\Phi \left(\vec{\beta }\right)\right],$$

$$\frac{\left(1-{\beta }^2\right)}{2n}\frac{\partial \Phi }{\partial \beta }$$

$$\frac{\left(1-{\beta }^{+2}\right)}{2n}\frac{\partial \Phi }{\partial {\beta }^{+}}$$

$$\Phi \left(\vec{\beta }\right)=-n\left[{\beta }^{+}\beta +\tilde{c}\ln (1-{\beta }^2)+h.c.\right],$$

$$\text{normalcolor}{P}_S\left(\vec{\beta }\right)=N{\left[(1-{\beta }^2{)}^{\tilde{c}}(1-{\beta }^{+2}{)}^{{\tilde{c}}^{\ast }}\exp (2{\beta }^{+}\beta )\right]}^n.$$

$$g^{(2)}(0)=\frac{⟨a^{†}{a}^{†}aa⟩}{⟨a^{†}a{⟩}^2},$$

Full text

I Introduction

The Schrödinger cat is a famous thought experiment (Schrödinger, 1935), where a cat is placed in a quantum superposition of two macroscopically distinct states, either alive or dead. It opens the fundamental question of whether quantum theory holds true in the macroscopic world (Haroche, 2013; Wineland, 2013; Arndt and Hornberger, 2014). Macroscopic superpositions have been experimentally realized in atoms (Monroe et al., 1996; Leibfried et al., 2005; Kovachy et al., 2015; Omran et al., 2019) and photons (Ourjoumtsev et al., 2007; Afek et al., 2010; Kirchmair et al., 2013), and have been proposed in quantum computation (Mirrahimi et al., 2014), quantum teleportation (van Enk and Hirota, 2001), quantum metrology (Joo et al., 2011) and quantum key distribution (Simon et al., 2014). One of the most common recent strategies for Schrödinger cats (Reid and Yurke, 1992) is via non-equilibrium subharmonic generation (Drummond et al., 1980, 1981) leading to discrete time symmetry-breaking or time crystals (Zhang et al., 2017), and this approach is analyzed in greater detail here.

The steady state of above-threshold subharmonic generation is known  for parametric down-conversion without anharmonicities (Drummond et al., 1980; Wolinsky and Carmichael, 1988). In this case transient Schrödinger cats are possible (Reid and Yurke, 1992; Krippner et al., 1994; Munro and Reid, 1995). Quantum subharmonic generation with Kerr anharmonicities was recently achieved in superconducting circuits (Leghtas et al., 2015), and large cat states were observed. In this experiment, the physics of the quantum steady state is different from previous studies (Sun et al., 2019). This exact solution for the steady state demonstrates how dissipation restores broken time symmetry, with potential applications to solving combinatorial optimization problems (McMahon et al., 2016).

Quantum optical and quantum circuit physics are similar, except that quantum circuits operate at microwave instead of optical frequencies. General driven quantum subharmonic generation with damping and weak nonlinearities was studied in a previous paper (Sun et al., 2019), where non-equilibrium quantum tunneling (Drummond and Kinsler, 1989) occurs. Here we focus on the cat-like properties of the steady states in the case of strong combined parametric and Kerr nonlinearities, as found in superconducting quantum circuits.

We analytically calculate the exact steady state in subharmonic generation with strong parametric and Kerr nonlinearities. This exactly soluble model has a very rich structure, while displaying the expected physics of more complex devices. We use the resulting exact correlation function to show that neither simple mixtures of coherent states nor Schrödinger cat states can occur in the steady state. This is confirmed by a numerical steady-state calculation in the number state basis.

We expect this physical result to occur in other parametric experiments with a similar dissipative, non-equilibrium behavior. A steady-state mixture of coherent states (Wolinsky and Carmichael, 1988) is achievable as a limiting case of extremely strong nonlinearities, but it is still a mixed state. This is consistent with the superconducting experiment (Leghtas et al., 2015) where an approximate Schrödinger cat was observed in a transient regime. The steady-state in the zero loss case can show a macroscopic superposition, although it is not unique, due to conserved number parity.

The outline of this paper is as follows. In Section (II) we explain our model definitions and notation, with a comparison to Josephson junction superconducting circuit theory. In Section (III) we obtain the exact steady-state solution, and explain the integration contour for the complex P-representation manifold. Section (IV) gives the diagonalization method as another alternative. In Section (V), moments are calculated both exactly and in approximations using coherent or incoherent combinations of coherent state delta-functions, for comparison purposes. Finally, Section (VI) summarizes our results.

II Combined nonlinearity model

Firstly, we summarize the system properties and theoretical techniques used previously (Drummond et al., 1981; Bartolo et al., 2016; Sun et al., 2019), with both Kerr and parametric nonlinearities. We then treat the detailed properties of the strongly coupled case. A schematic figure of the experimental system is shown in Fig. 1. The annihilation and creation operators of the $k$-th mode in two coupled resonant cavities are $a_{k},a_{k}^{\dagger}$ at frequencies $\omega_{k}$. The frequencies have been set as $\omega_{2}\simeq 2\omega_{1}$, so the system can be externally driven simultaneously at fundamental and subharmonic frequencies, with $2\omega_{0}$ and $\omega_{0}$, although we include detunings as well.

II.1 Hamiltonian

We assume a doubly resonant nonlinear cavity with a non-interacting Hamiltonian in the rotating frame of $H_{0}=\hbar\sum\Delta_{k}a_{k}^{\dagger}a_{k}$, where $\Delta_{k}=\omega_{k}-k\omega_{0}\ll\omega_{0}$ for input lasers’ frequencies of $\omega_{0}$ and $2\omega_{0}$. Then we will reduce the driving on the subharmonic mode to zero, thus only the fundamental mode is driven as in the experiment (Leghtas et al., 2015). The interaction Hamiltonian is assumed to be given by

[TABLE]

Here $\mathcal{E}_{2}$ is the envelope amplitude of the driving for the mode $a_{2}$, while $\kappa$, $\chi$ are the parametric and Kerr nonlinearities (Drummond and Walls, 1980) respectively. Kerr nonlinearities are only included for the mode $a_{1}$.

In addition, we include single-photon and two-photon losses in this open system. Defining $H=H_{0}+H_{I}$, the master equation for the density matrix $\rho$ is

[TABLE]

Here $\gamma_{k}^{(j)}$ are the relaxation rates for $j$-photon losses in the $k$-th mode, with no two-photon losses in mode $k=2$ for simplicity. The dissipative terms are

[TABLE]

where $\hat{O}=\hat{a}_{k}^{j}$. The corresponding thermal noises are set to zero. This allows us to study the steady-state properties in the low-temperature limit, in order to understand this exactly soluble case of maximal quantum coherence.

II.2 Effective Hamiltonian and master equation

We suppose the second harmonic mode is strongly damped, as in the recent Yale experiments (Leghtas et al., 2015), giving complex single-photon loss terms defined as $\gamma_{k}=\gamma_{k}^{(1)}+i\Delta_{k}$, with single-photon losses $\gamma_{k}^{(1)}$ and detunings $\Delta_{k}$ in the $k-$th mode. An adiabatic Hamiltonian is obtained for $a\equiv a_{1}$ as:

[TABLE]

The effective driving field $\epsilon$ and nonlinearity $\chi_{e}$ are:

[TABLE]

The master equation of the reduced density matrix $\rho_{1}=\text{Tr}_{2}(\rho)$ is then obtained as

[TABLE]

with an effective two-photon loss $\gamma_{e}^{(2)}$ , where

[TABLE]

Here we have taken the detuning $\Delta_{2}$ into account. Hence the expression of the effective parameters are slightly different from those in the previous work (Sun et al., 2019), while the master equation (6) takes the same general form.

II.3 Josephson model

In this subsection, we clarify the relations between the superconducting Josephson junction experiment (Leghtas et al., 2015) and our work. In the supplemental material of the experiment (Leghtas et al., 2015), the derivation of the system Hamiltonian, similar to ours (1), has been provided in detail. Here we will make a brief comparison, so that we can connect the parameters in our Hamiltonian (1) to those in the experiment (Leghtas et al., 2015).

In the experiment (Leghtas et al., 2015), two superconducting microwave oscillators were coupled through a Josephson junction. These oscillators are the fundamental modes of two superconducting cavities. One is a high Q cavity termed “the storage”, where the steady states formed. The other is a low Q cavity termed “the readout”, to evacuate entropy from the storage cavity. The system Hamiltonian of the qubit, the readout and storage modes reads

[TABLE]

Here $a_{m}$ is the annihilation operator for the qubit $m=q$, the readout mode $m=r$ and storage mode $m=s$, respectively, and $n_{m}=a_{m}^{\dagger}a_{m}$ is the corresponding number operator. $E_{J}$ is the Josephson energy, and $\varphi$ is the phase across the junction, which can be decomposed as the linear combination of the phase across each mode, with $\varphi_{m}$ denoting the contribution of mode $m$ to the zero point fluctuations of $\varphi$. The system is irradiated by the drive and pump inputs with complex amplitudes $\epsilon_{d}$ , $\epsilon_{p}$ and frequencies $\omega_{d}$, $\omega_{p}$, respectively.

In order to eliminate the system frequencies and the pump amplitude, we make use of the rotating frame of

[TABLE]

with $\tilde{\xi}_{p}\approx\xi_{p}e^{-i\omega_{p}t}$ and $\xi_{p}\approx-i\epsilon_{p}/\left(\frac{\kappa_{r}}{2}+i\left(\omega_{r}-\omega_{p}\right)\right)$. Thus the Hamiltonian takes the form of

[TABLE]

If we expand the term $\cos(\tilde{\varphi})$ up to the fourth order, and only keep non-rotating terms, the Josephson Hamiltonian then reads,

[TABLE]

with

[TABLE]

Here, the Hamiltonian $H_{\text{ Kerr }}$ corresponds to self-Kerr and cross-Kerr coupling terms, with $\chi_{mm}=\frac{E_{J}}{\hbar}\varphi_{m}^{4}/2$ and $\chi_{mm^{\prime}}=\frac{E_{J}}{\hbar}\varphi_{m}^{2}\varphi_{m^{\prime}}^{2}$. In the Hamiltonian $H_{2}$, the first two terms are nonlinear couplings between the storage and readout modes with $g_{2}=\chi_{sr}\xi_{p}^{*}/2$, which lead to the subharmonic generation. The other terms corresponds to the weak coherent drive $\epsilon_{d}$ on the readout mode.

II.4 Josephson parameters

In this paper we will focus on the evolution of the storage and readout modes. As given in the supplemental material of Ref.  (Leghtas et al., 2015), the Hamiltonian for the reduced system is

[TABLE]

where $\Delta_{d}=\omega_{r}-\omega_{d}-\delta_{r}-2\chi_{rr}\left|\xi_{p}\right|^{2}$ and $\Delta_{p}=-\Delta_{d}+2\left(\omega_{s}-\frac{\omega_{p}+\omega_{d}}{2}-\delta_{s}-\chi_{rs}\left|\xi_{p}\right|^{2}\right)$. In order to include the losses and the quantum noises, master equations have been analyzed in Ref. (Leghtas et al., 2015), where the single photon damping $\sqrt{\kappa_{r}}a_{r}$ and $\sqrt{\kappa_{s}}a_{s}$ have been considered.

In our notations, we set $a_{1}=a_{s}$ and $a_{2}=a_{r}$. Hence, this is similar to our initial Hamiltonian (1), with $\Delta_{1}=(\Delta_{p}+\Delta_{d})/2$, $\Delta_{2}=\Delta_{d}$, $\mathcal{E}_{2}=-i\epsilon_{d}$, $\kappa=2g_{2}$, $\chi=-\chi_{ss}$, $\gamma_{1}^{(1)}=\kappa_{s}/2$ and $\gamma_{2}^{(1)}=\kappa_{r}/2$. In our initial Hamiltonian (1), we have omitted the cross-Kerr term $\chi_{rs}$ and the self-Kerr term on the second harmonic mode $\chi_{rr}$ for simplicity.

In fact, the same approximation was used to derive the adiabatic Hamiltonian in Ref. (Leghtas et al., 2015) as well. In their supplemental material, they have shown that the effect of the cross-Kerr term is negligibly small and thus can be ignored. Since our main results are obtained under the adiabatic approximation, these omissions are valid in our situation.

With the detunings and $\chi_{rr}$ omitted, the adiabatic approximation can be applied in the region where $g_{2}/\kappa_{r},\epsilon_{d}/\kappa_{r},\chi_{rs}/\kappa_{r}\sim\delta$ and $\chi_{ss}/\kappa_{r},\kappa_{s}/\kappa_{r}\sim\delta^{2}$ with the small dimensionless parameter $\delta\ll 1$. By neglecting terms of order $\delta$ and higher, the adiabatic Hamiltonian has been derived in the supplemental material of Ref. (Leghtas et al., 2015), which reads

[TABLE]

The corresponding master equation takes the form:

[TABLE]

with $\kappa_{2}=4\left|g_{2}\right|^{2}/\kappa_{r}$ and $\epsilon_{2}=-2ig_{2}\epsilon_{d}/\kappa_{r}$. Compared with our adiabatic Hamiltonian (4) and master equation (6), we find the parameter mappings for this experiment to be: $\mathcal{E}=2\epsilon_{2}$, $\chi_{e}=-\chi_{ss}$ and $\gamma_{e}^{(2)}=\kappa_{2}$ with $\gamma_{1}^{(2)}=0$ and $\Delta_{1}=\Delta_{2}=0$.

III Exact steady-state solution

This master equation has an exact analytic solution for the steady-state, including damping, driving and detunings together with all the nonlinear couplings. We note that this is neither an energy eigenstate nor a thermal state, but rather a unique non-equilibrium solution to the steady-state.

III.1 Complex P-representation

To obtain the exact solution, we introduce a generalized P-representation (Drummond and Gardiner, 1980) transformation of the single-mode density matrix. If we expand the reduced quantum density matrix in terms of coherent state projection operators and a complex P-distribution $P\left(\alpha,\alpha^{+},t\right)$ , one then obtains

[TABLE]

where $\left|\alpha\right\rangle$ is a coherent state and $d\alpha d\alpha^{+}$ is a surface integral measure over a closed surface, so that boundary terms will vanish on integration by parts. The adiabatic Hamiltonian results in a single-mode Fokker-Planck equation for $P$,

[TABLE]

where we define $\gamma\equiv\gamma_{1}=\gamma_{1}^{(1)}+i\Delta_{1}$. We also introduce an effective complex nonlinear decay of $g=\gamma_{e}^{(2)}+i\chi_{e}$, and a function $\mathcal{E}\left(\alpha\right)=\mathcal{E}-g\alpha^{2}$. The notation $h.c.$ indicates hermitian conjugate terms obtained by the replacement of $\alpha\rightarrow\alpha^{+}$, and the conjugation of all complex parameters. As in our previous work (Sun et al., 2019), we introduce dimensionless parameters: $\epsilon=\mathcal{E}/g$, $n=\left|\epsilon\right|$, $c=\gamma/\left(gn\right)$, $\tau=\mathcal{E}t$, $\beta=\alpha/\sqrt{\epsilon}$, and $e^{i\theta}=g/|g|=n/\epsilon$, so that the Fokker-Planck equation can be simplified to the form:

[TABLE]

With this transformation, time is scaled relative to the two-photon driving rate. Here $c$ is a complex dimensionless single-photon loss and detuning, and $n$ is the photon number at which saturation of the mode occupation occurs due to the nonlinear losses.

The steady-state solution of the scaled Fokker-Planck equation (18) can be derived via the potential method (Graham and Haken, 1971a, b; Risken, 1972; Seybold and Risken, 1974)

[TABLE]

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

[TABLE]

These equations (20) are obtained by inserting the form (19) into the Fokker-Planck equation (18) and requiring that $\partial P_{1}/\partial\tau=0$ in the steady-state.

By solving the differential equations (20) directly, the exact steady-state solution with quantum noise can be expressed via the potential:

[TABLE]

with $\tilde{c}=c-1/n$. Thus, the steady-state probability distribution is

[TABLE]

This is the exact zero-temperature steady-state solution for the density matrix. Written in this way, we can see how it scales with the effective driving field $n$ occurring in the exponent. Apart from $n$, all the parameters here can have complex values, which is necessary when treating the situations in recent quantum circuit experiments (Leghtas et al., 2015).

In the case where the power $\tilde{c}$ has a negative real part, if a real planar complex manifold is chosen, one obtains singular peaks at the boundaries where $\left|\beta\right|,\left|\beta^{+}\right|=\pm 1$, as shown in Fig (2). This would give boundary terms on partial integration, causing errors. Integration over phase-space distributions requires vanishing boundary terms. Instead, one must choose a curved topological structure with cuts on the complex integration manifold. This leads to branch points, rather than local potential minima. This is why there is no quantum tunneling, although transient Schrödinger cats can be formed in this type of experiment (Leghtas et al., 2015).

As a result, this physical situation requires a completely different phase-space manifold to that investigated in the previous work (Sun et al., 2019), where the real part of $\tilde{c}$ is positive. In that case, there is quantum tunneling between local potential minima on a finite, bounded manifold. To define the distribution for strong coupling, one must choose complex integration contours which are closed, continuous (Drummond and Gardiner, 1980; Drummond et al., 1981; Bartolo et al., 2016), and without boundaries. This is obtained by inserting cuts at the branch-points for $\beta=\pm 1$ and $\beta^{+}=\pm 1$, combined with complex Pochhammer contours. This method is used to represent the beta and hypergeometric special functions (Pochhammer, 1890; MacRobert, 1938; Bateman, 1955). One way to visualize this is to imagine the contours drawn on both sides of two sheets of paper, one for $\beta$ and one for $\beta^{+}$.

III.2 Moments and Correlations

The second-order correlation function of the single-mode intra-cavity field is defined as

[TABLE]

where the $k$-th moment can be calculated with P-representation integrals as

[TABLE]

It is well known that nonclassical effects like photon anti-bunching will occur if $g^{(2)}(0)<1$ and classical bunching takes place if $g^{(2)}(0)>1$. Thus, $g^{(2)}(0)$ is often used to distinguish classical from non-classical behavior  (Scully and Zubairy, 1997).

The exact solution for the moments (Drummond et al., 1981) is obtained by expanding the term $e^{2n\beta^{+}\beta}=\sum_{m}\left(2n\right)^{m}\beta^{m}\beta^{+m}/m!$ in Eq. (22). In this way, we obtain the form of moment after normalization and integration over the complex manifold, as:

[TABLE]

Here ${}_{2}F_{1}$ is the hypergeometric function, and $N^{\prime}$ is the normalization factor,

[TABLE]

The case of real $\tilde{c}$ has been investigated in Ref (Wolinsky and Carmichael, 1988; Krippner et al., 1994), where there was no anharmonic nonlinearity, and a real manifold was used. It was suggested that the steady-state distribution approaches a set of $\delta$ functions in strong coupling limits. The case without single-photon loss and anharmonic nonlinearity has also been studied in Ref. (Gilles et al., 1994), where one always has $\tilde{c}=-1/n$. In this case, steady-state Schrödinger cats can be achieved with initial Fock states. Other work studying this potential in different parameter regimes was used to benchmark our numerical results, given below (Bartolo et al., 2016).

IV Numerical diagonalization

As a comparison and independent check of these exact results, we have also solved the master equation Eq. (6) numerically by expanding the density operator in a number state basis. The steady state of the system corresponds to the eigenstate of the Liouvillian operator with zero eigenvalue. This steady state density operator is then used to compute the statistical moments of interest. In this approach, which is valid for small photon number, we numerically diagonalize the Liouville operator of the master equation, with a photon number cutoff. This allows us to compare the analytical and numerical approaches. We find that there is excellent agreement between the two methods.

IV.1 Number state basis

In order to verify our analytic results, we applied these numerical number state methods to the same case. We expand the density operator $\rho$ in the number-state basis, where its matrix elements $\rho_{kl}$ are defined as

[TABLE]

Then the master equation (6) takes the form:

[TABLE]

Here the Einstein summation convention has been used on identical indices. And $T_{ij}^{kl}$ is a four-dimensional transition matrix, which describes the transition from the state $\rho_{kl}$ to the state $\rho_{ij}$. It can be written as

[TABLE]

with

[TABLE]

The system can be characterized by the eigenvectors of the transition matrix $T_{ij}^{kl}$. The steady state of the system corresponds to the eigenvector with zero eigenvalue (Kinsler and Drummond, 1991).

IV.2 Transition matrix elements

Within the numerical calculation, we must use a photon number cutoff $N$ to make the transition matrix finite, $0\leq i,j,k,l\leq N$. This approximation is valid if the high-photon-number states play negligible roles in determining the system’s evolution. We check that the cut-off is set to a high enough value by repeating the calculation with a higher cutoff and checking that no change occurs.

Hence, the four-dimensional matrix $T_{ij}^{kl}$ can be reduced to a two-dimensional one $T_{\bar{\alpha}}^{\bar{\beta}}$ with this truncation, so that

[TABLE]

with

[TABLE]

where

[TABLE]

Here $\delta_{\bar{\alpha}}^{\bar{\beta}}$ is a Kronecker delta, and $\bar{\alpha}$, $\bar{\beta}$ are in the range of $[1,(N+1)^{2}]$.

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

[TABLE]

where the coefficients $A_{k}$ define the initial state. We order the indices $k$ by the size of the real part of the eigenvalues, $\Re(\epsilon_{k})\geq\Re(\epsilon_{k+1})$. Therefore, $\epsilon_{0}$ is the stable eigenvalue with $\epsilon_{0}=0$, and $\rho_{\bar{\alpha}}^{(0)}$ corresponds to the stable state.

With the numerical expansion method, the stable state $\rho_{\bar{\alpha}}^{(0)}$ can be obtained by solving the eigenvalue problem of the transition matrix $T_{\bar{\alpha}}^{\bar{\beta}}$. Then the average photon numbers $\left\langle a^{\dagger}a\right\rangle$ and the second order correlation functions $g^{(2)}(0)$ can be obtained directly by

[TABLE]

where $\rho^{(0)}$ is the stable-state matrix reshaped from the stable-state vector $\rho_{\bar{\alpha}}^{(0)}$. The numerical results are shown below in Fig. 3 and 4 below with green dots. They agree with the analytic results very well. This confirms the validity of our analytic calculations.

V Moments and Schrödinger Cat

comparisons

We will use these exact analytic and approximate numerical results to check the validity of approximate delta-function steady-state distributions which we introduce below (LABEL:distribution). These correspond to the physical assumptions that one has either a quantum superposition or a quantum mixture of two coherent states with opposite signs. As we show below, neither assumption is correct in the steady-state of this driven, non-equilibrium quantum system.

V.1 Experimental parameter values

For numerical evaluations of the steady-state moments, we obtain the parameters of the recent experiment (Leghtas et al., 2015), using the results of Section (II). In our notation, we obtain that for these recent quantum circuit experiments, $\gamma/2\pi=3.98$kHz, $g/2\pi=(7.96-4i)$ kHz and $\mathcal{E}=(-19.2-0.07i)$ kHz. Thus, we have $\tilde{c}=-0.279+0.093i$ and $\epsilon=-1.92-0.97i$. Since the real part of $\tilde{c}$ is negative, there will be singularities occurring at $\beta=\pm 1$ or $\beta^{+}=\pm 1$.

From now on, we will treat the strong coupling regime, which corresponds to the parameter region of $\Re(\tilde{c})<0$. Using the definitions of $\tilde{c}$ and $g$, we have

[TABLE]

Considering that $n>0$, it follows that $\Re(\tilde{c})<0$ is equivalent to $\gamma_{e}^{(2)}(\gamma_{1}^{(1)}-\gamma_{e}^{(2)})+\chi_{e}(\Delta_{1}-\chi_{e})<0$. This is satisfied if there is either a weak single-photon damping $\gamma_{1}^{(1)}$ or strong nonlinear couplings $\chi_{e}$, $\gamma_{e}^{(2)}$. It is easily checked, provided there are no detunings, that $\Re(\tilde{c})\geq-1/n$ and the limit $\tilde{c}\to-1/n$ occurs if $\gamma_{1}^{(1)}\ll\gamma_{e}^{(2)}$ or $\gamma_{1}^{(1)}\ll\chi_{e}$.

Considering that nonlinear losses are always weak, the relation $\gamma_{1}^{(1)}\ll\gamma_{e}^{(2)}$ can occur with large $\kappa$ refer to Eq. (7). Thus the limit $\tilde{c}\to-1/n$ occurs either with large nonlinearities $\kappa$ or $\chi$.

V.2 Delta-function approximations

To understand the physics more clearly, we note that in the limit of $\tilde{c}\rightarrow-1/n$, the exact solution is a product of simple poles with opposite contour integration directions. These can be integrated using Cauchy’s theorem, and correspond to a delta-function solution, so the ratio of the probabilities at the singularities is

[TABLE]

If we assume this is also true approximately for $\tilde{c}\neq-1/n$, we obtain a real distribution (Wolinsky and Carmichael, 1988) in the form of

[TABLE]

We now contrast this with an idealized, even cat state $\left|\psi\right\rangle_{cat}\varpropto\left[\left|\sqrt{\epsilon}\right\rangle+\left|-\sqrt{\epsilon}\right\rangle\right]$, where the P-representation takes the form after normalization

[TABLE]

The factor is $e^{-2n}$ ($e^{2n}$), rather than $e^{-4n}$ ($e^{4n}$) in Eq. (LABEL:distribution), so even if the steady state does evolve to a delta-function distribution (LABEL:distribution), it will be a mixed state instead of a true cat state.

In this case, the density matrix can be derived to have the following form,

[TABLE]

Here

[TABLE]

The purity of this limiting form can then be obtained as

[TABLE]

which is a monotonic decreasing function of $n$ since

[TABLE]

for $n>0$. Thus, the driving will weaken the purity of the steady state since $n$ is proportional to the driving $\mathcal{E}_{2}$.

It is obvious that we will have $p\to 1$ in the limit of $n\to 0$. Thus the delta-function distribution tends to be a true Schrödinger cat state in this limit. However, since $|\epsilon|=n\to 0$, the steady state will actually reduce to a vacuum state. This is natural that a non-driven system can be expressed as a vacuum state. In the opposite limit of $n\to\infty$, the delta-function steady-state distribution (LABEL:distribution) will reduce to the mixed state $\rho_{mix}$ since $p\to 0$. Therefore, a pure Schrödinger cat state is unreachable in the steady state of the system, even using an approximate delta-function solution.

The parity $\hat{\mathcal{P}}=(-1)^{a^{\dagger}a}$ can also be studied directly with the complex P-distribution (LABEL:distribution). In the P-representation, the parity operator is equivalent to the average of $\mathcal{P}=\exp(-2n\beta^{+}\beta)$. In the steady state of the delta-function approximation, we have $\mathcal{P}_{ss}=\text{sech}(2n)$. This means that $\mathcal{P}_{ss}=1$ in the case of $n=0$, and $\mathcal{P}_{ss}=0$ in the limit of $n\to\infty$. It is consistent with the density matrix (40) which is a vacuum state when $n=0$ and a mixed state when $n\to\infty$. Parity is not conserved because of the finite single-photon loss.

V.3 Steady-state distributions

The exact steady-state distributions (22) with different parameters are shown in Fig. 2, plotted on a finite manifold. We see that delta-function distribution will be obtained approximately with large $|\Re(\tilde{c})|$ and small $|\Im(\tilde{c})|$, and reduced to classical mixture of coherent states with large $n$. However, these graphs also demonstrate that the probability does not vanish at the boundaries, which means that with $\Re\left(\tilde{c}\right)<0$ on this bounded manifold, the potential solution when restricted to this planar manifold is no longer a solution to the original master equation, since boundary terms from integration by parts are non-vanishing.

An inspection of Fig 2 shows that when assuming a real, bounded manifold, the distribution is neither a true delta function, nor does it vanish at the boundaries, which is the reason why the exact complex contour manifold is essential when there are poles.

V.4 Moment comparisons

As a result, the true steady states are clearly neither mixtures of delta functions nor Schrödinger cats. This difference can be quantified by using the steady-state distribution (LABEL:distribution), to compare moments. The approximate $k$-th moment is obtained directly with the definition (24) as,

[TABLE]

Similarly, the moment can be written down directly with the cat state (LABEL:distribution-cat) as:

[TABLE]

We have compared the average steady-state photon number $\left\langle a^{\dagger}a\right\rangle$ and the second order correlation function $g^{(2)}(0)$ changing with $c$ in Fig. 3. The results of Fig. 3, show that the delta-function distribution (LABEL:distribution) is only attainable when $\tilde{c}\to-1/n$, which is valid when $\gamma_{1}^{(1)}\ll\gamma_{e}^{(2)}$ or $\gamma_{1}^{(1)}\ll\chi_{e}$, if there are no detunings. Mathematically, it is obtained by reaching the steady state first and then taking the limit $\gamma_{1}^{(1)}\to 0$, which is different from the magenta circles where we take $\gamma_{1}^{(1)}=0$ exactly and then get the steady states assuming some particular parity (Gilles et al., 1994). Number parity is conserved only if $\gamma_{1}^{(1)}=0$, and non-conserved if $\gamma_{1}^{(1)}\neq 0$. Thus the ordering of the limit is important, which leads to the gap between the red line with $\tilde{c}\to-1/n$ (a mixed state) and the magenta circles (a pure cat state) in Fig. 3. In addition, the delta-function distribution can also be obtainable in the region of extremely strong nonlinearity as the limit $\tilde{c}\to-1/n$ suggests, which is more practical than the case $\gamma_{1}^{(1)}=0$.

In Fig. 3 the results of the delta-function distributions never agree with those of the cat states. This is consistent with the discussion above that the steady state of the system is always a mixed state (40) instead of a pure cat state. Although there are crosses for the exact results of the steady state and those of the pure cat state, they are always at different $\tilde{c}$ for $\left\langle a^{\dagger}a\right\rangle$ and $g^{(2)}(0)$. The exact steady state is therefore different from both the cat state and a mixture of delta-functions. Hence we can’t generate a pure steady-state cat state, unless the system has no single-photon losses.

We have stated that in the limit of small $n$, the delta-function distribution (LABEL:distribution) tends to an approximate Schrödinger cat. Now we show how $\left\langle a^{\dagger}a\right\rangle$ and $g^{(2)}(0)$ change with $n$ in Fig. 4. It is natural that the average photon number $\left\langle a^{\dagger}a\right\rangle$ increases with large driving $\mathcal{E}_{2}\propto n$ as shown in Fig. 4(a). It also shows that in the region of small $n$, their photon numbers agree with each other, but $g^{(2)}(0)$ has a different behavior.

This means that even with $n\to 0$, the delta-function steady-state distribution (LABEL:distribution) is still different from the distribution of a Schrödinger cat. We also show in Fig. 4 that in the limit of $\tilde{c}\to-1/n$, the exact steady state will approach the delta-function steady-state distribution, although as before, this is not a cat state.

It is directly checked with Eqs. (44) and (45) that the second-order correlation functions are

[TABLE]

Thus in the limit of $n\to 0$, we have $g_{cat}^{(2)}(0)/g_{lim}^{(2)}(0)\to 4$ with $g_{cat}^{(2)}(0)\to\infty$ and $g_{lim}^{(2)}(0)\to\infty$. This tendency can be found in the Fig. 4. In addition, we will also have $g_{cat}^{(2)}(0)>g_{lim}^{(2)}(0)>1$ over the full range of $n$. This means that their probability distributions are both super-Poissonian (Scully and Zubairy, 1997). From all the discussions above, we demonstrate that the delta-function steady-state distribution (LABEL:distribution) is different from the Schrödinger cat state, even if $n\to 0$.

Pure steady-state cats can occur in systems without single-photon loss and anharmonic nonlinearity (Gilles et al., 1994). If we neglect the single-photon loss in our system from the beginning, the steady-state solution is obtained from solving $\partial\rho_{1}/\partial t=0$ in Eq. (6). We expand the density operator in the coherent state basis as $\rho_{1}\left(t=\infty\right)=\iint c_{\alpha,\alpha^{\prime}}|\alpha\rangle\langle\alpha^{\prime}|\,d^{2}\alpha d^{2}\alpha^{\prime}$. Substituting into Eq. (6) with $\gamma_{1}^{(1)}=0$, for arbitrary $c_{\alpha,\alpha^{\prime}}$ we have

[TABLE]

Thus the steady-state density matrix with no single-photon damping takes the form,

[TABLE]

where the coefficients $c_{\alpha,\alpha^{\prime}}$ are determined by the initial states. This is consistent with earlier work (Gilles et al., 1994), which however had no Kerr anharmonic term. In the P-representation, the distribution reads in this undamped case,

[TABLE]

which is also a delta-function distribution. The possible pure state solutions are coherent states and cat states. Since the parity is conserved without single-photon loss according to the master equation (6), Schrödinger cats can be achieved if the initial states are eigenstates of the parity, such as Fock states. These steady-state Schrödinger cats with $\gamma_{1}^{(1)}=0$ and initial vacuum states have been graphed in Fig. 3 (a) and (b), where a gap between them and the results for the limit $\gamma_{1}^{(1)}\to 0$, which is a mixture, can be observed.

VI Summary

We have studied the steady states of quantum subharmonic generation with strong nonlinearity, which has been experimental achieved (Leghtas et al., 2015). By comparing the correlation functions, we conclude that true Schrödinger cats cannot survive in the steady state unless there is no single-photon loss. With single-photon loss included, the steady state for subharmonic generation will reduce to a delta-function steady-state distribution (LABEL:distribution) only if there is an extremely strong nonlinearity. More generally, the exact solution is always more complex than any type of delta-function, whether a pure or mixed state. To obtain this exact behavior, the correct integration manifold is a Pochhammer contour which samples both sheets of a double Riemann sheet contour. Intriguingly, this reflects some of the character of the transient macroscopic superposition that occurs on the path to the steady-state.

Acknowledgements.

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