Photon pair condensation by engineered dissipation

TL;DR

This paper demonstrates how engineered photon-pair dissipation can stabilize a unique quantum condensate with long-range pair correlations, using circuit-QED systems to realize such states.

Contribution

It introduces a method to stabilize photon-pair condensates with phase-nematic order through engineered dissipation, including a scheme to select a specific parity sector.

Findings

Stabilization of a photon pair condensate with long-range order

Implementation scheme in circuit-QED architecture

Control of parity sectors via Glauber dynamics

Abstract

Dissipation can usually induce detrimental decoherence in a quantum system. However, engineered dissipation can be used to prepare and stabilize coherent quantum many-body states. Here, we show that by engineering dissipators containing photon pair operators, one can stabilize an exotic dark state, which is a condensate of photon pairs with a phase-nematic order. In this system, the usual superfluid order parameter, i.e. single-photon correlation, is absent, while the photon pair correlation exhibits long-range order. Although the dark state is not unique due to multiple parity sectors, we devise an additional type of dissipators to stabilize the dark state in a particular parity sector via a diffusive annihilation process which obeys Glauber dynamics in an Ising model. Furthermore, we propose an implementation of these photon-pair dissipators in circuit-QED architecture.

Equations

$$\frac{d}{dt}\rho =-i[H,\rho ]+\mathcal{L}\rho ,$$

$$l_j=(a_j^{†2}+a_{j+1}^{†2})(a_j^2-a_{j+1}^2),$$

$$⟨a_i^{†2}{a}_j^2⟩\text{xLongrightarrow}|i-j|\to \infty ⟨a_i^{†2}⟩⟨a_j^2⟩={\bar{n}}^2⟨e^{i2({\varphi }_i-{\varphi }_j)}⟩={\bar{n}}^2,$$

$$⟨a_i^{†}{a}_j⟩=\bar{n}⟨e^{i({\varphi }_i-{\varphi }_j)}⟩=\bar{n}⟨e^{i({\varphi }_i-{\varphi }_j)+\pi }⟩=-⟨a_i^{†}{a}_j⟩=0,$$

$$|D_{2n}⟩\propto A^{†n}|0⟩,$$

$$c_j=\Gamma (a_{j+1}^{†}{a}_j+a_{j-1}^{†}{a}_j)\left(\frac{1-P_j}{2}\right),$$

$$H^{\prime }=-\frac{\chi }{2}\sum _{\lambda =L,R}{a}_{\lambda }^{†2}{a}_{\lambda }^2+g_1{T}_{-}|e⟩⟨g|+g_2{T}_{+}^{†}|f⟩⟨e|+h.c.$$

$$H_{\text{eff}}=-\frac{\chi }{2}\sum _{\lambda =L,R}{a}_{\lambda }^{†2}{a}_{\lambda }^2+\frac{1}{\delta }\left[\begin{array}{cc}|g_2{|}^2{T}_{+}^{†}{T}_{+}& g_1{g}_2{T}_{+}^{†}{T}_{-}\\ g_1^{\ast }{g}_2^{\ast }{T}_{-}^{†}{T}_{+}& |g_1{|}^2{T}_{-}^{†}{T}_{-}\end{array}\right]$$

$$|D_{2}m⟩=A^{†m}|0⟩,$$

$$A^{†}=\sum _{i=1}^L{a}_i^{†2}\frac{1}{n_i+1},$$

$$[a_j^2,A^{†}]$$

$$=\sum _{n_j=0}^{\infty }(n_j+2)|n_j⟩⟨n_j|-\sum _{n_j=0}^{\infty }(n_j+2)|n_j+2⟩⟨n_j+2|$$

$$=\sum _{n_j=0}^{\infty }(n_j+2)|n_j⟩⟨n_j|-\sum _{n_j=2}^{\infty }{n}_j|n_j⟩⟨n_j|$$

$$=2I+|1_j⟩⟨1_j|=2I+P_j^1,$$

$$[l_j,A^{†}]=[a_j^2-a_{j+1}^2,A^{†}]=P_j^1-P_{j+1}^1.$$

$$P_j^1{A}^{†}=|1_j⟩⟨1_j|\sum _{i=1}^L\sum _{n_i=0}^{\infty }\sqrt{\frac{n_i+2}{n_i+1}}|n_i+2⟩⟨n_i|=0.$$

$$l_j{A}^{†}|0⟩=(A^{†}{l}_j+P_j^1-P_{j+1}^1)|0⟩=0.$$

$$l_j(A^{†}{)}^{k+1}|0⟩=(A^{†}{l}_j+P_j^1-P_{j+1}^1)A^k|0⟩=0.$$

$$|D_2^{\prime }m⟩=\prod _{i=1}^{2n_d}{a}_{d_i}^{{}^{\prime }†}(A^{†}{)}^{n-n_d}|0⟩,$$

$$a_{d_i}^{{}^{\prime }†}=a_{d_i}^{†}\frac{1}{n_{d_i}+1}.$$

$$[l_i,a_j^{\prime }]$$

$$=\sum _{n_j=1}^{\infty }\sqrt{n_j}|n_j-1⟩⟨n_j|-\sum _{n^{\prime }=0}^{\infty }\sqrt{n_j^{\prime }+2}|n_j^{\prime }+1⟩⟨n_j^{\prime }+2|$$

$$=|0_j⟩⟨1_j|.$$

$$l_j|D_2^{\prime }m⟩$$

$$=\prod _{i=1}^{2n_d}[a_{d_i}^{{}^{\prime }†}{l}_j+{\delta }_{j,d_i}|0_j⟩⟨1_j|$$

$$-{\delta }_{j+1,d_i}|0_{j+1}⟩⟨1_{j+1}|](A^{†}{)}^{n-n_d}|0⟩$$

$$=0.$$

$$\frac{d}{dt}P(\{\sigma \},t)=-\sum _j{w}_j(\{{\sigma }_{j-1},{\sigma }_j,{\sigma }_{j+1}\})P(\{\sigma \},t)$$

$$+\sum _j{w}_j(\{{\sigma }_{j-1},-{\sigma }_j,{\sigma }_{j+1}\})P(\{{\sigma }_1,...,-{\sigma }_j,{\sigma }_{j+1},...\},t).$$

$$w_j(\{\uparrow \uparrow \downarrow \})=w_j(\{\downarrow \downarrow \uparrow \})=w_j(\{\downarrow \uparrow \uparrow \})=w_j(\{\uparrow \downarrow \downarrow \})=\frac{1}{2}\Gamma .$$

$$w_j(\{\uparrow \downarrow \uparrow \})=w_j(\{\downarrow \uparrow \downarrow \})=\Gamma .$$

$$w_j(\{\uparrow \uparrow \uparrow \})=w_j(\{\downarrow \downarrow \downarrow \})=\Gamma h.$$

$$w_j({\sigma }_j)=\frac{1}{2}\Gamma \left[1-\frac{1}{2}\gamma {\sigma }_j({\sigma }_{j+1}+{\sigma }_{j-1})\right],$$

$$H=-J\sum _j{\sigma }_j{\sigma }_{j+1}.$$

$$\frac{p_j(-{\sigma }_j)}{p_j({\sigma }_j)}=$$

Full text

---

With the rapid development of quantum optical technology and quantum information platforms such as cavity/circuit quantum electrodynamics (QED) Blais et al. (2007); Schoelkopf and Girvin (2008); Houck et al. (2012) and Rydberg polaritons Saffman et al. (2010), it is now possible to investigate strongly-correlated many-body physics of photons Houck et al. (2012); Hartmann et al. (2008); Carusotto and Ciuti (2013); Noh and Angelakis (2016); Cárdenas-López et al. (2019). While photons can have strong interactions in these platforms, they do not naturally thermalize, and one has to synthesize thermalization and a chemical potential to obtain many-body ground states Hafezi et al. (2015); Kapit et al. (2014); Ma et al. (2019); Schiró et al. (2012); Wang et al. (2019, 2018). Remarkably, dissipation induced by the environment, which is usually regarded as a noise source leading to decoherence of the states, can actually become a useful resource. If harnessed properly, dissipation can be used to autonomously prepare and stabilize an exotic many-body pure state as the steady/dark state of a system Plenio et al. (1999); Diehl et al. (2008); Kraus et al. (2008a); Verstraete et al. (2009); Weimer et al. (2010); Cho et al. (2011); Kastoryano et al. (2011); Diehl et al. (2011); Barreiro et al. (2012); Reiter et al. (2013); Rao and Mølmer (2014); Budich et al. (2015); Bardyn et al. (2012); Stannigel et al. (2014). In the context of analogue quantum simulation, some well-known examples of dissipative engineering schemes include the autonomous preparation and stabilization of the Bose-Einstein condensate (BEC) state Diehl et al. (2008), the Majorana-fermion state Diehl et al. (2011) and the Chern insulator state Budich et al. (2015), all of which can be thought as the ground states of non-interacting Hamiltonians. In a digital quantum simulation scheme Barends et al. (2015, 2016); Salathé et al. (2015), a class of jump operators has been realized Kraus et al. (2008b); Weimer et al. (2010); Barreiro et al. (2011), where the steady states correspond to the ground states of a specific class of interacting Hamiltonians.

Meanwhile, there have been significant experimental achievements in engineering analogue dissipators with higher-order photon jumps in small circuit-QED systems Murch et al. (2012); Shankar et al. (2014); Mirrahimi et al. (2014); Leghtas et al. (2015); Mundhada et al. (2017); Albert et al. (2018); Touzard et al. (2018). While these efforts have been motivated by autonomous error correction for a single- or two-sites system, it is interesting to investigate engineering many-body states, using these tools. Specifically, one can ask whether a strongly-correlated pure many-body state can be stabilized with an engineered analogue dissipator? In this Letter, we answer this question by proposing a type of two-photon jump operator which can dissipatively prepare and stabilize an exotic strongly correlated photon-pair condensate exhibiting phase-nematic order. Furthermore, we propose an analogue experimental realization with circuit-QED systems. The added benefit of this approach is that one does not require effective thermalization or generation of a chemical potential in a photonic system.

In order to illustrate the key idea, we start with a canonical example of Ref. Diehl et al. (2008). For an open quantum system with Markovian environment, the system dynamics can be described by the Lindblad master equation:

[TABLE]

where $H$ is the Hamiltonian of the system and the Liouvillian $\mathcal{L}\rho=\sum_{j}\kappa_{j}(2l_{j}\rho l^{\dagger}_{j}-l^{\dagger}_{j}l_{j}\rho-\rho l^{\dagger}_{j}l_{j})$ describes the dissipation associated with the jump operator $l_{j}$ with decay rate $\kappa_{j}$. Consider the dynamics of bosonic particles on 1-dimensional lattice with $H=0$ and a number-conserving jump operator of the form l_{j}$$=$$(a^{\dagger}_{j}+a^{\dagger}_{j+1})(a_{j}-a_{j+1}), connecting nearest neighbors in the lattice, where $a^{\dagger}_{j}$ is the boson creation operator on site $j$. These jump operators stabilize a BEC with fixed number of particles, i.e., a pure dark state |\,D\,\rangle=(a_{k=0}^{\dagger})^{N_{tot}}|\,0\,\rangle\propto\big{(}\sum_{j}a^{\dagger}_{j}\big{)}^{N_{tot}}|\,0\,\rangle, where $N_{tot}$ is the total number of bosons. A simple way to understand these jump operators is the mean-field picture in which $a_{j}\rightarrow\sqrt{\bar{n}}e^{i\phi_{j}}$, where $\phi_{j}$ represents the compact U(1) phase variable (mod $2\pi$).

The dark-state condition $l_{j}|\,D\,\rangle=0$ gives rise to the mean-field solution: $\phi_{j+1}-\phi_{j}=0$ (mod $2\pi$), suggesting a phase locking between neighboring sites. In this case, the dark state has long-range order, i.e., $\langle a^{\dagger}_{i}a_{j}\rangle$ $\xLongrightarrow{|i-j|\rightarrow\infty}$ $\langle a^{\dagger}_{i}\rangle\langle a_{j}\rangle$ $=\bar{n}$, where we have introduced mean-field order parameter $\langle a_{j}\rangle$ $\approx$ $\sqrt{\bar{n}}\langle e^{i\phi}\rangle$ and $\phi$ is the uniform phase, after the spontaneous breaking of a $U(1)$ symmetry. While this order parameter is fragile in 1D, these ideas can be generalized to higher dimensions where the long-range order can become robust.

Pair jump operators.— In this work, we propose a quartic jump operator connecting site $j$ and $j+1$ of the form

[TABLE]

in a 1D lattice, as shown in  Fig. 1(a). This can be generalized to 2D and 3D by assigning a jump on each link of the lattice. Before deriving the exact form of the wave function, we consider a mean-field solution in which we take $a_{j}\rightarrow\sqrt{\bar{n}}e^{i\phi_{j}}$. The dark-state condition $l_{j}|\,D\,\rangle=0$ hence gives rise to the mean-field solution: $2(\phi_{j+1}-\phi_{j})=0$ (mod $2\pi$), suggesting a locking of twice the phase variables between neighboring sites. This leads to the mean-field order parameter $\langle a_{j}^{2}\rangle=\bar{n}\langle e^{2i\phi_{j}}\rangle$ rather than $\langle a_{j}\rangle$, as in the previous case. In fact, we now have $\langle a_{j}\rangle=0$. For correlation functions, we get

[TABLE]

for $\ i\neq j$. While $(2\phi_{j})$ exhibits long-range order, $\phi_{j}$ does not since $\phi_{j}$ can flip by $\pi$ and still satisfy the dark-state condition [c.f. Fig. 1(b,c)]. This photon pair condensate exhibits phase-nematic order. A similar state has been studied in the context of Josephson junction arrays Douçot and Vidal (2002),the symmetry breaking phase of a photon pair hopping Hamiltonian Lebreuilly et al. (2019) and fragmented many-body state in the ultra-cold atomic system Mueller et al. (2006); Bader and Fischer (2009); Fischer and Kang (2015). The oriented rods without an arrow head in Fig. 1(c) represent the local order parameter $\langle e^{i2\phi_{j}}\rangle$ for such a state, which does not differentiate the $\pi$-phase flip of $\phi_{j}$ and corresponds to the spontaneous breaking of a $U(1)/\mathbb{Z}_{2}$ symmetry.

Exact solutions.— For the system where the Hamiltonian $H=0$ and the jump operator is described in Eq. (2).The steady state density matrix is given by $\rho_{ss}=|\,D\,\rangle\langle\,D\,|$ where $|\,D\,\rangle$ is annihilated by the jump operators in Eq. (2) satisfying $l_{j}|\,D\,\rangle=0.$

We find that the dark state $|\,D_{2n}\,\rangle$ can be described as a condensate of $n$ two-photon bound states:

[TABLE]

where $A^{\dagger}=\sum_{j}a_{j}^{{\dagger}2}(n_{j}+1)^{-1}$ 222 See Supplemental Material for more details on the derivation dark state exact solution. is the creation operator of quasi-particles related to photon pair bound state and $n_{j}=a^{\dagger}_{j}a_{j}$ is the on-site number operator. Note that the extra normalization factor $(n_{j}+1)^{-1}$ in the definition creation operator $A^{\dagger}$ only affects the relative weights of different photon pair spatial configurations, but not the essence of the pair condensation.

One can easily see that the single-particle correlator $\langle a^{\dagger}_{i}a_{j}\rangle$ (for $i\neq j$) vanishes because $a_{i}|\,D_{2n}\,\rangle$ and $a_{j}|\,D_{2n}\,\rangle$ have zero overlap since the photon occupation on site $i$ and $j$ becomes odd respectively. On the other hand, the pair correlation $\langle a^{{\dagger}2}_{i}a_{j}^{2}\rangle$ is flat since a_{i}^{2}|\,D_{2n}\,\rangle$$=$$a_{j}^{2}|\,D_{2n}\,\rangle\propto|\,D_{2(n-1)}\,\rangle due to the fact that taking a pair out of the condensate at any site results in the same condensate with $n-1$ pairs of photons. This is just a manifestation of the definition of a pair condensate.

We numerically simulate the time-dependent master equation Eq. (1) for an open 1D chain via quantum trajectory method with time-evolving block decimation (TEBD) algorithm Gambetta et al. (2008); Vidal (2004), with results shown in Fig. 2. We start with a product Fock state $|\,2,0,2,0,...\,\rangle$, and the jump operator drives the system to the steady (dark) state. We see from Fig. 2(a) that the single-particle correlation function $\langle a^{{\dagger}}_{L/4}a_{L/4+10}\rangle$ remains zero at all times, while the pair correlation function $\langle a^{{\dagger}2}_{L/4}a_{L/4+10}^{2}\rangle$ grows rapidly with an exponential saturation until reaching the steady state. The whole time evolution resembles a cooling process. The cooling time is independent of the system size, as seen in the plot where the total number of sites is varied as $L=16,24,32$. The exponential saturation behavior and the cooling time is manifest on a logarithmic scale, see Fig. 2(b).

We also plot the pair correlators as a function of the distance between two sites, i.e., $\langle a^{{\dagger}2}_{L/4}a_{L/4+j}^{2}\rangle$ versus time $t$, as shown in Fig. 2(c). We see that the state has almost flat correlation when reaching the dark state, consistent with the prediction from the analytical solution shown above. Before reaching the dark state, the correlator is not flat and decays with distance. This is due to the fact that correlation between more distant sites needs more time to be built up. Fig. 2(d) shows the equilibrium time $T_{eq}$ as function of distance. The equilibrium time $T_{{\rm eq}}$ is defined as the time it takes for the correlator $\langle a^{{\dagger}2}_{L/4}a^{2}_{L/4+j}\rangle$ to reach $80\%$ of its steady state value. The spreading of the correlation function follows the Lieb-Robinson light cone beheavior. In addition, we have observed that the introduction of a Kerr non-linearity in the form $H=Ua^{{\dagger}2}_{j}a_{j}^{2}$ in the system Hamiltonian leads to an exponential decay of correlator (in 1D) as a function of the distance $j-1$. The decay becomes faster increasing $U$.

Parity sectors.—The above analytical and numerical analyses only consider a simplified situation where the initial condition has all even number of photons. We note that even for fixed total photon number, the dark-state subspace has extensive degeneracies $2^{L-1}$. ($L$ is the number of sites), labeled by the local parity $P_{j}=(-1)^{n_{j}}$ on each site. The exact wave function we wrote down above in Eq. (4) is only the exact wave function for the sector where the parities of all sites are all even, i.e., $P_{j}=1$ for all $j$, which we call a “pure pair condensate”. On top of that, there are odd-parity “defects”, which are created in pairs from the pure pair condensate.

The wave function of a particular defect configuration can be given by $|\,D^{\prime}_{n_{d}}\,\rangle\propto\prod_{i=1}^{2n_{d}}a^{\prime{\dagger}}_{d_{i}}|\,D_{2(n-n_{d})}\,\rangle$, where $n_{d}$ denotes the number of pairs of odd-parity defects and their positions are labeled by $d_{i}$. Several solutions of the parity-sector problem are discussed as follows.

To begin with, we note that the different parity sectors are not coupled together via the jump operators. Similar to what has been considered in the numerical simulation in Fig. 2, one can start with an initial product state in the all-even sector (easy to prepare experimentally with pulses in the presence of onsite nonlinearity). In this case, the jump operator will only drive the system to the dark state $|\,D_{2n}\,\rangle$ in Eq. (4), in the absence of unwanted noise. We note that noise is always present in experimental systems which either brings the state to different sectors. Therefore, one only expects to prepare the targeting dark state with the jump operators before the unwanted decoherence dominates. If one aims to stabilize the dark state, extra measurement or stabilizing schemes are needed as discussed below.

A more general solution is: by imposing measurement and feedback operation on the parity of each site, i.e., $P_{j}=(-1)^{n_{j}}$, it is possible to keep projecting the many-body state to a particular parity sector. A jump operator describing such measurement and feedback operation to stabilize parity-all-even sector is

[TABLE]

where $\Gamma$ is the hopping rate. Note that this jump operator applies a hopping term connecting that particular site to its nearest neighbors conditioned by the the parity on the particular site being odd. It causes the odd-parity defect to take a random walk and eventually annihilate with another parity defect, as illustrated in Fig. 3(a). This diffusive defect annihilation process resembles a chemical reaction described by the formula: $\mathbf{df}+\mathbf{df}\rightarrow 0$, where $\mathbf{df}$ stands for a single defect. Therefore, no defect will exist in the steady state if the total photon number is even ($2n$), and so the steady state becomes a pure pair condensate. We call such a process “healing”.

Since the parity measurement at time $t+dt$ will post-select the direction (left or right) to which the defect has hopped at time $t$, the defect dynamics can be exactly mapped to a classical stochastic dynamics of the diffusive annihialation problem. For a 1D chain, the dynamics of defect density exhibits power law decay: $m(t)$ $\sim$ =$(\Gamma t)^{-1/2}$ Lee (1994); Racz (1985a)

The classical Monte Carlo simulation (with 100 sites) of the time evolution of the average defect density $m(t)$ quantitatively agrees with the quantum TEBD simulation (with 30 sites) for a 1D chain, as shown in Fig. 3(b). The former is plotted in a log-log scale in Fig. 3(c), where the associted fitting curves confirms an asymptotic power-law $m(t)$ $\sim(\Gamma t)^{-1/2}$. We have checked in both types of simulations that the time evolution of $m(t)$ is almost independent of system size.

In the presence of additional noise, such as incoherent single-photon hopping described by the jump operator l^{\prime}_{j}$$=$$a^{\dagger}_{j}a_{j+1}, there exists a finite defect production rate $\Gamma h$. The defect production can be balanced by the diffusive defect annihilation process with hopping rate $\Gamma$, and the defect density approaches a residual steady-state density $m_{s}$ with a characteristic relaxation time $\tau$. We have the following scaling (when $h\ll 1$): $m_{s}\sim h^{1/\delta}$ and $\tau\sim\Gamma^{-1}h^{-\Delta}$ Racz (1985a). Note that even a large hopping rate $\Gamma$ cannot reduce the steady-state density $m_{s}$, but only the healing time $\tau$. A generic scaling law $\Delta\xi=1/\delta$ should be satisfied Racz (1985b). For a 1D chain, we have $\delta=2$ and $\Delta=1$ 333See Supplemental Material which includes Ref. Racz (1985a); Glauber (1963); Felderhof and Suzuki (1971); Racz (1985b); Torney and McConnell (1983); Lee (1994).

In Fig. 3(d), we use TEBD to calculate the steady-state pair correlators $\langle a^{{\dagger}2}_{L/4}a_{L/4+j}^{2}\rangle$ for system sizes $L=16,24,32$, in the following situations: i. ideal case (no noise: $h=0$); ii. dirty case (in the presence of single photon hopping noise: $h>0$, without healing: $\Gamma=0$); iii. healed case (with noise and healing: $h,\Gamma>0$). We see that in the presence of noise which proliferate the parity defects, the steady state ends up with a mixed state and the correlators decays exponentially, while the healing process significantly slows down the decay.

Due to the complexity of the parity operator, one typically can only realize such a jump operator with an active parity measurement in circuit-QED setup through either continuous Cohen et al. (2017) or discretized repeated Ofek et al. (2016); Hu et al. (2018) measurement schemes, instead of using continuous autonomous stabilization. Nevertheless, in the situation that we first impose hard-core condition for occupation more than three photons, i.e., $a^{{\dagger}3}_{j}=0$, the parity condition can be simply converted to occupation condition, and we can effectively re-express the jump operator in Eq. (5) as c^{\prime}_{j}$$=$$(a^{\dagger}_{j}a_{j+1}$$+$$\text{H.c.})n_{j}(n_{j}-2), up to a constant factor of 2. This jump operator can potentially be implemented continuously and hence autonomously stabilize the targeting pure photon pair condensate. Similarly, one can either actively or autonomously monitor and stabilize the total photon number in the system.

Finally, we can also passively stabilize the parity sector via energetic constraint in the case of hard-core condition ${a^{\dagger}}^{3}_{j}=0$. This is achieved by assign the following energy penalty term $\delta H=-V_{0}n_{j}(n_{j}-2)$ in the Hamiltonian $H$ Bretheau et al. (2015). Therefore, the configuration with single-photon occupation (odd parity) on any site is projected out of the low-energy sector.

Experimental realization with Circuit QED.—We illustrate the experimental scheme with a two-site jump operator l$$=$$(a_{L}^{{\dagger}2}+a_{R}^{{\dagger}2})(a_{L}^{2}-a_{R}^{2}). The generalization to a 1D chain is straightforward. Consider a system consisting of the two high Q cavities $a_{L}$ and $a_{R}$, and an anharmonic oscillator. The anharmonic oscillator is modeled by a three-level system ($|\,g\,\rangle,|\,e\,\rangle$ and $|\,f\,\rangle$) and is coupled to a cavity $a_{\lambda}$ and a Josephson junction mode $J_{\lambda}$, where $\lambda=L,R$, at both sides respectively. Both junction modes are driven by a two-tone drive $\Omega_{\lambda}(t)=\Omega_{\lambda,r}e^{i\omega_{r}t}+\Omega_{\lambda,b}e^{i\omega_{b}t}$ as shown in Fig. 4.

We engineer a two-photon jump operator via four-wave mixing induced from the junction modes $J_{L}$ and $J_{R}$. The drive $\omega_{r}$ ($\omega_{b}$) is used to introduce an exchange of two photons of cavity mode $a^{2}_{\lambda}$ ($a^{{\dagger}2}_{\lambda}$) with the excitation $g\rightarrow e$ ($e\rightarrow f$). The four-wave mixing interaction of the pump $\omega_{r}$ ($\omega_{b}$) is proportional to $\sum_{\lambda}\Omega_{\lambda,r}a^{2}_{\lambda}|\,e\,\rangle\langle\,g\,|$ ($\sum_{\lambda}\Omega_{\lambda,b}a^{{\dagger}2}_{\lambda}|\,f\,\rangle\langle\,e\,|$). The minus sign in the jump operator can be engineered by introducing a $\pi$ phase shift between $\Omega_{L,r}$ and $\Omega_{R,r}$. The effective Hamiltonian is of the form

[TABLE]

where $T_{\pm}=a^{2}_{L}\pm a^{2}_{R}$, $\chi$ is the Kerr nonlinearity induced from the junction modes and $g_{1}$ ($g_{2}$) is proportional to the Rabi frequency $\Omega_{r}$ ($\Omega_{b}$) as shown in Fig. 4a.

To obtain the jump operator in Eq. (2), we combine the two-photon loss process ($T_{-}$) and two-photon creation process ($T_{+}^{\dagger}$) in equation above. We can achieve this by detuning the two four-wave mixing processes by $\delta$ as shown in Fig. 4b so that only a cascade of two such processes is possible Albert et al. (2018). The detunning $\delta\gg g_{1},g_{2}$ allows a two-photon exchange processes via a Raman transition. The effective Hamiltonian becomes

[TABLE]

where the 2 $\times$ 2 matrix acts on the anharmonic oscillator basis $\{|\,f\,\rangle,|\,g\,\rangle\}$. In Eq. (6), the term $T^{\dagger}_{+}T_{-}|\,f\,\rangle\langle\,g\,|$ gives the desired two-photon process coupled to $g\leftrightarrow f$ transition. Assuming the decay rate($\kappa$) of process $f\rightarrow g$ is much greater than all of the other coupling constant in Eq. (6). The system can be described $H_{s}=-\frac{\chi}{2}\sum_{\lambda=L,R}a_{\lambda}^{\dagger 2}a_{\lambda}^{2}+\frac{|g_{1}|^{2}}{\delta}(a_{L}^{\dagger 2}-a_{R}^{\dagger 2})(a_{L}^{2}-a_{R}^{2})$ and the jump operator $l_{s}=\sqrt{\frac{2}{\kappa}}\frac{g_{1}g_{2}}{\delta}(a_{L}^{\dagger 2}+a_{R}^{\dagger 2})(a_{L}^{2}-a_{R}^{2})$. The jump operator gives rise to Eq. (1) for the array case. The self-Kerr and the cross-Kerr nonlinearity terms in Hamiltonian $H_{s}$ can be eliminated by adding an extra pair of Josephson junction and a two-level system 444 See Supplemental Material for more details..

Conclusion and outlook.— We have discovered a photon pair jump operator which can dissipatively prepare and stabilize an exotic two-photon pair-condensate with phase-nematic order, with a circuit-QED implementation. We have further proposed a conditional hopping operator to stabilize the dark state in a particular parity sector. Such a scheme can also be realized with Rydberg polaritons or ion-trap systems. An interesting future direction would be using such higher-order dissipators for autonomous quantum error correction using bosonic codes.

Appendix A Exact Solution of the Dark State

A.1 Pure Two-Photon Condensate

In this section, we are going to show that the analytical expression of the dark state of the jump operators $l_{j}=(a_{j}^{{\dagger}2}+a_{j+1}^{{\dagger}2})(a_{j}^{2}-a_{j+1}^{2})$, where $j$ is integer and $1\leq j\leq L$, is

[TABLE]

where

[TABLE]

and the number of photon is $2m$ .

We start by introducing two useful identities. (1) $[l_{j},A^{\dagger}]=P^{1}_{j}-P^{1}_{j+1}$. (2) $P^{1}_{j}A^{\dagger}=0$, where $l_{j}=a_{j}-a_{j+1}$ is the annihilation part of the jump operator $l_{j}$ and $P^{1}_{j}$ is the projector that project a state to one photon Fock state at site j, or $P^{1}_{j}=|\,1_{j}\,\rangle\langle\,1_{j}\,|$.

Proof of identity (1) : The commutator $[a^{2}_{j},A^{\dagger}]$ is

[TABLE]

where $I$ is the identity operator.

Therefore, we have the commutator

[TABLE]

Proof of identity (2) : Since the $A^{\dagger}$ operator always creates more than 1 photons, $A^{\dagger}|\,s\,\rangle$ is always orthogonal to $|\,1_{j}\,\rangle$ for any arbitrary state $|\,s\,\rangle$ and hence $P^{1}_{j}A^{\dagger}|\,s\,\rangle$ is always zero. Or equivalently,

[TABLE]

With these two useful identities, we can proceed to proof that Eq. ((7)) is a dark state under jump operator $l_{j}$.

Proof: When $m=1$, we have

[TABLE]

Here we use the identity (1).

Assume that for an integer $k$, $l_{j}(A^{\dagger})^{k}|\,0\,\rangle=0$.

Then for $m=k+1$, we have

[TABLE]

The first term in the right hand side of above equation is zero because we assume $l_{j}(A^{\dagger})^{k}|\,0\,\rangle=0$. The second and the third terms are zero because of the identity (2).

By mathematical induction, equation $l_{j}(A^{\dagger})^{m}|\,0\,\rangle=0$ holds for all integer $m$ greater than zero. And therefore, $|\,D_{2}m\,\rangle$ is the dark state under the dissipation $l_{j}$.

A.2 With Parity Defects

The dark state in the Eq (7) is not unique in the absence of the healing process. One can always create localized parity defects as described in the main text. The dark state with parity defects of a particular defect configuration can be written as

[TABLE]

where $2n_{d}$ is the number of defect, $d_{i}$ is the position of the ith defect and

[TABLE]

To show that this is also a dark state, let’s consider the commutator $[l_{i},a^{\prime}_{j}]$.

If $i\neq j$ and $i\neq j-1$, then these two operators commute with each others. When $i=j$, we have

[TABLE]

Also, when $j=i+1$, the commutator becomes $[l_{i},a^{\prime}_{j}]=-|\,0_{j}\,\rangle\langle\,1_{j}\,|$.

With this commutator, we can see that

[TABLE]

The first term in the last equation above is zero because of the identity (1) and the second term and the third term vanishes since the state $(A^{\dagger})^{n-n_{d}}|\,0\,\rangle$ is always orthogonal to $|\,1_{j}\,\rangle$ by parity symmetry.

Appendix B Effective Classical Theory and Simulation for the Healing Process and the Steady States with Parity Defects

In this section, we discuss the effective classical theory of the diffusive annihilation (healing) process of parity defects mentioned in the main text.

We start the discussion by focusing on the 1D situation, while the generalization to higher dimension is straightforward and will be discussed later. We label the defect number on a 1D periodic lattice with $L$ sites by \{\mathbf{m}\}$$\equiv$$\{m_{1},...,m_{j},m_{j+1},...,m_{L}\}, where $j$ labels the sites. For later convenience, we map the diffusive annihilation process onto a periodic classical Ising spin chain with $L$ sites as in Ref. Racz (1985a). The spin configuration is labeled by \{\sigma\}$$\equiv$$\{\sigma_{1},...,\sigma_{j},\sigma_{j+1},...,\sigma_{L}\}, consisting of the stochastic spin variables $\sigma_{j}=\pm 1$ ($1$ stands for $\uparrow$ and $-1$ stands for $\downarrow$). The correspondence between the spin and defect configuration is through the mapping $m_{j}=(1-\sigma_{j}\sigma_{j+1})/2$. The presence of a defect at $j$ (i.e., $m_{j}=1$) is equivalent to the presence of a domain wall between opposite spins at site $j$ and $j+1$ on the spin chain, i.e., $\uparrow_{j}\downarrow_{j+1}$ or $\downarrow_{j}\uparrow_{j+1}$. One can thus think the defect is living on the bond of the effective spin chain. The spin flip on the left/right of the domain wall corresponds to the moving of defect (domain wall) towards left/right, as illustrated by the processes $\uparrow\uparrow\downarrow\downarrow\Rightarrow\uparrow\downarrow\downarrow\downarrow$ and $\uparrow\uparrow\downarrow\downarrow\Rightarrow\uparrow\uparrow\uparrow\downarrow$ respectively. The spin flip in the middle of two domain walls correspond to the annihilation of two defects (domain walls), as illustrated by the processes $\downarrow\uparrow\downarrow\Rightarrow\downarrow\downarrow\downarrow$ and $\uparrow\downarrow\uparrow\Rightarrow\uparrow\uparrow\uparrow$. The inverse of such annihilation processes gives rise to the pair production processes of defects (domain walls).

The probability distribution functions of the effective spin configuration are represented as $P(\{\sigma\},t)$. There are $2^{L}$ such distribution functions in total, and they satisfy the following classical master equation first introduced by Glauber Glauber (1963):

[TABLE]

Here, the first term on the right hand side describes the transition from the current configuration $\{\sigma_{j}\}$ at time $t$ to the new configuration $\{\sigma_{1},...,-\sigma_{j},\sigma_{j+1},...\}$ with the $j^{\text{th}}$ spin being flipped, while the second term describe the transition from the configuration $\{\sigma_{1},...,-\sigma_{j},\sigma_{j+1},...\}$ to the current configuration $\{\sigma_{j}\}$ by flipping the $j$th spin. Therefore, all the stochastic dynamics are captured by a single spin flip. The flipping rate ($w_{j}$) of the $j$th spin is determined by the local spin configuration involving the sites $j$ and $j\pm 1$. First, we require the defect to hop either to the left or right with hopping rate $\frac{1}{2}\Gamma$, leading to the following conditions:

[TABLE]

Next, we require neighboring defects to annihilate with each other with rate $\Gamma$, which is given by the defect hopping rate and the fact that double occupation of the defects on the same site is equivalent to zero defect. This leads to

[TABLE]

Finally, we require the pair production rate of defects to be $h$, leading to the condition

[TABLE]

Therefore, our classical stochastic model for the diffusive annihilation process on the 1D lattice is clearly exactly by Eq. (B) to Eq. (20).

Although we will consider analytical analysis of the master equation later, we can first simulate this effective dynamics numerically with the classical Monte Carlo simulation. Starting from a random spin configuration $\{\sigma\}$, we randomly update the spin configuration with the three types of random processes (defect hopping, pair annihilation, and pair production) with the corresponding rates given by Eq. (19) to Eq. (21). By averaging multiple stochastic histories, we can evaluate the physical variables, such as the average defect density $m$ as shown in Fig. 3 in the main text (where we have chosen $h=0$, $L=100$, several values of $\Gamma$ and averaged over $1000$ histories). Through fitting, we get the asymptotic power-law decay of the defect numbers $m(t,h=0)$ $\sim(\Gamma t)^{-1/2}$.

Now we consider the analytical analysis, following the treatment in Ref. Racz (1985a) and Glauber (1963). We first propose an approximated form of the transition rates

[TABLE]

where we choose $\gamma=(1-h)/(1+h)$. With this formula, we get the following conditions for the rates: (1) $w_{j}(\{\uparrow\uparrow\downarrow\})=w_{j}(\{\downarrow\downarrow\uparrow\})=w_{j}(\{\downarrow\uparrow\uparrow\})=w_{j}(\{\uparrow\downarrow\downarrow\})=\frac{1}{2}\Gamma$, (2) $w_{j}(\{\uparrow\downarrow\uparrow\})=w_{j}(\{\downarrow\uparrow\downarrow\})=\frac{1}{2}\Gamma(1+\gamma)=\Gamma/(1+h)$, and (3) $w_{j}(\{\uparrow\uparrow\uparrow\})=w_{j}(\{\downarrow\downarrow\downarrow\})=\frac{1}{2}\Gamma(1-\gamma)=\Gamma h/(1+h)$. Note that condition (1) here is always the same as Eq. (19) in our exact classical stochastic model. On the other hand, in the $h=0$ ($\gamma=1$) case, conditions (2) and (3) are also the same as Eq. (20) and  (21) in the exact classical stochastic model. For small $h$, condition (2) and (3) remain a good approximation to the exact conditions Eq. (20) and  (21). One can think of this approximation as introducing additional effective short-range repulsion between the defects which slightly decreases the defect annihilation rate. This short-range repulsion is expected to not change the scaling property in the $h\rightarrow 0$ limit.

The reason to choose this approximated form in Eq. (22) is due to its equivalence to a kinetic Ising model as first pointed out by Glauber Glauber (1963). The Hamiltonian of the Ising model is

[TABLE]

When reaching equilibrium at temperature $T$, the ratio between the reduced probability distributions of the $j$th spin is determined by the Boltzman distribution, i.e.,

[TABLE]

On the other hand, the approximated stochastic model defined by Eq. (B) and Eq. (22) will approach a steady-state distribution with the ratio

[TABLE]

By comparing Eq. (20) with Eq. (19), we see that they are exactly the same with the correspondence

[TABLE]

This suggests that the approximated stochastic model is equivalent to a kinetic Ising model, whose steady-state distribution is the same as the equilibrium distribution at the effective temperature $T$ according to Eq. (26).

With the above equivalence, we can get steady state distribution for our exact stochastic model in the small $h$ limit via the equilibrium distribution. Due to the mapping m_{j}$$=$$(1$$-$$\sigma_{j}\sigma_{j+1})/2, we can first evaluate the equilibrium distribution of the correlation function in the Ising model, which has the exact solution

[TABLE]

Therefore, in the small $h$ limit, we get the average defect density at steady state as

[TABLE]

Since the correlation $\sigma_{j}\sigma_{j+1}$ is proportional to energy, its relaxation time $\tau$ can be determined by the long-time decay of homogenous energy perturbation as studied in Ref. Felderhof and Suzuki (1971), leading to the result

[TABLE]

As we have mentioned in the main text, in general we have the following three asymptotic scaling laws for the diffusive annihilation process: m(t,h=0)$$\sim$$(\Gamma t)^{-\xi}, $m_{s}\sim h^{1/\delta}$ and $\tau\sim\Gamma^{-1}h^{-\Delta}$. In addition, there is a scaling relation $\Delta\xi=1/\delta$, derived by Ref. Racz (1985b). From the property of equilibrium and kinetic Ising models, we have just obtained $\delta=2$ and $\Delta=1$. Using the scaling relation, we can derive that $\xi=1/2$, which has been confirmed by our classical Monte Carlo simulation discussed above and shown in Fig. 3 in the main text. We also note that exact analytical derivation of $\xi$ can be found in Ref. Torney and McConnell (1983).

We also note that in the presence of noise inducing pair production of defects ($h>0$), the steady state of the system is a mixed state of different defect configurations, of which the density matrix in the small-$h$ limit can be described by the equilibrium density matrix of the effective Ising model at the corresponding effective temperature $T$ according to Eq. (26). The detailed pair condensate description of each defect configuration is encoded by the wavefunctions in Eq. (LABEL:defectwavefunction) in the main text.

So far, we have focused in the context of 1D chain, while we note such a stochastic model and the corresponding master equations can be straightforwardly generalized to higher dimension, as discussed in Ref. Lee (1994). We note for dimension equal or larger than the critical dimension, i.e., $d\geq d_{c}=2$ Lee (1994), a simple mean-field description for the diffusive annihilation process gives the correct scaling m(t,h$$=$$0)$$\sim$$(\Gamma t)^{-1} as we have discussed in the main text.

Appendix C Derivation for the Circuit-QED Realization

The experimental setup is shown in Fig. 5. The two-level system and the two junction modes $J_{3}$ and $J_{4}$ are added to cancel the Kerr nonlinearity in the cavity $a_{L}$ and $a_{R}$.

We start by writing the main results of this derivation. The effective dynamics of the system is

[TABLE]

The coupling constant $g_{1}$, $g_{2}$ and $g_{3}$ will be defined later. All the non-linear terms in the Hamiltonian can be canceled with each others by choosing $\chi=\frac{|g_{1}|^{2}}{2\delta_{1}}=\frac{|g_{3}|^{2}}{2\delta_{2}}$. We can achieve this by tuning the driving frequency $\omega_{r}$,$\omega_{b}$ and $\omega_{T}$. Therefore we can obtain a pure dissipative process that generate photon pair condensate as its steady state.

The Hamiltonian of the two cavities, anharmonic oscillator and the Josephson junction modes.

[TABLE]

where $\omega_{c}$ is the resonance frequency of the cavity, $\omega_{j}$ is the energy of anharmonic oscillator which is modeled by the three-level system ($|\,g\,\rangle$, $|\,e\,\rangle$ and $|\,f\,\rangle$), $\omega_{TLS}$ is the resonance frequency of the two-level system and $\omega_{J}$ is the resonance frequency of the Josephson junction modes and we assume that these four Josephson junction modes have the same frequency for simplicity.

The phases cross the $i$th junction modes are

[TABLE]

where $\varphi_{a}$, $\varphi_{an}$ and $\varphi_{J}$ are the amplitude participation ratios of the respective modes in the junction and $\Sigma$ is the annihilation operator of the anharmonic oscillator and we only consider the lowest three levels. Therefore, $\Sigma=|\,1\,\rangle\langle\,0\,|+\sqrt{2}|\,2\,\rangle\langle\,1\,|$. The driving fields are defined as

[TABLE]

C.1 Displaced Frame

In this section ,we follow the derivation in [1] to move the contribution of the drivings to the phase degree of freedom, we perform the following transformations:

• •

1. Go into the rotating frame of tone $\omega_{r}$. The junction modes becomes $J_{1(2)}\rightarrow J_{1(2)}e^{-i\omega_{r}t}$.

• •

2. Displace the two junction modes by $J_{1(2)}\rightarrow J_{1(2)}-\frac{\Omega_{r,1(2)}}{\omega_{J}-\omega_{r}}$

• •

3. Move back to the original frame. $J_{1(2)}\rightarrow J_{1(2)}e^{i\omega_{r}t}$.

After the transformation, the contribution of driving $\omega_{r}$ is completely absorbed by the phase $\varphi^{\prime}_{i}(t)=\varphi_{\lambda}a_{\lambda}+\varphi_{an}\Sigma_{an}+\varphi_{J}(J_{i}-\frac{\Omega_{\lambda,1}e^{i\omega_{d1}t}}{\omega_{c}-\omega_{d1}})+h.c.$ for i = 1,2. . We can perform the same transformation for driving $\omega_{b}$ and $\omega_{T}$ and expand the cosine term to the fourth order. The resulting Hamiltonian becomes

[TABLE]

where

[TABLE]

and

[TABLE]

C.2 Expansion of $[\tilde{\varphi}_{\lambda}(t)]^{4}$

We now expand the $[\tilde{\varphi}_{\lambda}(t)]^{4}$ in order to obtain the desired interaction in the non-rotating frame. After expansion, we perform normal-ordering to the expansion terms, and take all Stark shifts and Lamb shifts into account $\omega\rightarrow\tilde{\omega}$.

By tuning the driving

[TABLE]

and performing rotating wave approximation, we have

[TABLE]

where

[TABLE]

where the Kerr-nonlinearities are $\chi=E_{J}\varphi_{c}^{4}$, $\chi_{a,an}=E_{J}\varphi_{c}^{2}\varphi_{an}^{2}$ and $\chi_{a,T}=E_{J}\varphi_{c}^{2}\varphi_{T}^{2}$ .

C.3 Large Detuning Limit

Next, when $\delta_{1},\delta_{2}\gg g_{1},g_{2},\chi$, we can adiabatically eliminate the level $|\,f\,\rangle$ and $|\,1\,\rangle$ by performing the Schrieffer-Wolff transformation. The unitary operator is

[TABLE]

By keeping the terms up to $O(\frac{g_{1(2)}^{2}}{\delta})$ and neglecting the terms proportional to $|\,e\,\rangle\langle\,e\,|$ and $|\,1\,\rangle\langle\,1\,|$, we have

[TABLE]

where

[TABLE]

C.4 Large $\kappa_{f}$ Limit

Then we take the dissipation process of the anharmonic oscillator $|\,f\,\rangle\rightarrow|\,g\,\rangle$ and $|\,1\,\rangle\rightarrow|\,0\,\rangle$ into account. The dynamics of the system can be described by the master equation

[TABLE]

where $\kappa_{f}$ is the decay rate from state $|\,f\,\rangle$ to $|\,g\,\rangle$ . Assume that the decay rate $\kappa_{f}\gg\frac{|g_{1}|^{2}}{\delta_{1}},\frac{|g_{2}|^{2}}{\delta_{1}},\frac{|g_{1}g_{2}|}{\delta_{1}},\frac{|g_{3}|^{2}}{\delta_{2}}$, we can derive the effective dynamics of the two cavities.

Consider an operator $\hat{O}$ that only depends on the degrees of freedom of cavities. The equation of motion of $\hat{O}$ can be written as

[TABLE]

where

[TABLE]

Also, the dynamics of $|\,g\,\rangle\langle\,f\,|$ is

[TABLE]

Since $\kappa_{f}$ is much greater than all the other parameters in the effective Hamiltonian, we can assume that $\langle|\,g\,\rangle\langle\,g\,|\rangle\approx 1$, $\langle|\,f\,\rangle\langle\,f\,|\rangle\ll 1$ and $\langle|\,f\,\rangle\langle\,g\,|\rangle\approx O(\frac{g^{2}}{\kappa_{f}\delta})$. In this regime, the operator $|\,g\,\rangle\langle\,f\,|$ reaches its stationary state in the time scale that is much smaller than the dynamics of the cavity operator $\hat{O}$. Thus, we can replace the operator $|\,g\,\rangle\langle\,f\,|$ by its stationary state

[TABLE]

Therefore, Eq.(42) becomes

[TABLE]

The equation of motion can be converted to the master equation of the cavity mode as described in Eq.(30).