Dissipative Phase Transition in Systems with Two-Photon Drive and Nonlinear Dissipation near the Critical Point

Abstract

In this paper, we examine dissipative phase transition (DPT) near the critical point for a system with two-photon driving and nonlinear dissipations. The proposed mean-field theory, which explicitly takes into account quantum fluctuations, allowed us to describe properly the evolutionary dynamics of the system and to demonstrate new effects in its steady-state. We show that the presence of quantum fluctuations leads to a power-law dependence of the anomalous average at the phase transition point, with which the critical exponent is associated. Also, we investigate the effect of the quantum fluctuations on the critical point renormalization and demonstrate the existence of a two-photon pump “threshold”. It is noteworthy that the obtained results are in a good agreement with the numerical simulations.

1. Introduction

Nowadays, dissipative phase transition (DPT) in open quantum systems is among the most rapidly developing fields of quantum optics [1]. DPT can be observed when a direct manipulation of the interaction constants, external driving, or dissipation rates of the system lead to an abrupt and nonanalytical change of the system observables [2,3]. Recent publications have reported on the observation of dissipative critical phenomena and nonequilibrium quantum states in superconducting circuits [4,5], cavity quantum electrodynamics systems [6,7,8], optomechanical resonators [9], semiconductor microcavities [10], and atomic systems [11,12,13,14]. The experimental implementation of highly controllable open nonequilibrium photonic systems became possible due to nonlinear reservoir engineering. It enables the realization of the optical cavities with an engineered two-photon drive and nonlinear dissipation [15,16,17,18,19,20]. In such systems, the presence of dissipation does not destroy but stabilizes the quantum state, also known as the Schrödinger cat state [21,22,23]. Based on this state, it is possible to prepare a dynamically protected qubit [24] for further applications in quantum information processing [25,26,27,28,29,30].

The above-mentioned quantum systems are commonly studied with numerical methods. The most notable of these are integration of a master equation on a truncated Fock basis [19] and diagonalization of the Liouvillian superoperator [2,31,32]. Additionally, the complex-P-representation [33], Monte Carlo [34], quantum trajectory [35], and the quantum-absorber methods [36] are intensively utilized to investigate the properties of the nonequilibrium stationary state in such problems. However, numerics typically cannot allow one to uncover physical phenomena underlying the evolution of the system. This is one of the reasons why the DPT in open quantum systems is not well understood to date. Within the limit of a large number of photons, a qualitative picture of the occurring phenomena can be obtained using the so-called semiclassical approximation [37,38]. Nevertheless, this approach completely neglects the quantum fluctuations (QF), which can play a significant role even in the case of a large average number of photons in the system [10]. Therefore, the study of mean-field models, including QF, becomes an important and promising direction for further research, since it can provide simple and convenient analogies with the physics of equilibrium phase transitions. Besides, such a treatment can offer a correct approach for constructing theories beyond the mean field description.

In this paper, we demonstrate the mean field treatment of a DPT in systems with two-photon driven and nonlinear dissipation, which includes QF. Here, we use the formalism of Keldysh Green’s functions [39,40], since it explicitly takes into account the effects of QF. Recently it has been shown that the nonequilibrium Keldysh technique is a promising way to study nonequilibrium open quantum systems [41]. It also provides a theoretical framework for the systematic treatment of DPT [41,42,43,44,45]. The application of this approach enables us to construct the self-consistent equations of motion similar to the Gorkov equations in the mean-field approximation. We use the resulting equations to calculate the dynamics of the system observables and to demonstrate the new effects appearing in the steady-state. In analogy with Landau theory, we show that the presence of QF leads to a broadening of the DPT near the critical point. This effect is caused by an “external” QF field. Its existence leads to the power-law dependence of the anomalous average vs the pump rate at the phase transition point with which the critical exponent is associated. In addition, the QF effect causes the critical point renormalization. Both results are in a good agreement with our numerical simulations.

2. Nonlinear Cavity Including Two-Photon Processes

First, we briefly revisit the system described in [24] and experimentally realized in [46]. Its schematic is depicted in Figure 1a, which shows two superconducting microwave cavities. The first resonator has a low Q-factor associated with the presence of significant single-photon losses at κ rate. The second resonator has a high Q-factor, so it can be used to store and protect quantum information [23,27]. A transmission line with embedded Josephson junction provides a nonlinear interaction between two cavities, which we denote as readout and storage. In addition, pump and drive microwave tones are applied to the readout cavity at frequencies ω_p and ω_r. The Hamiltonian of the system has the following form (ℏ = 1):

$$H=H_0+H_{\mathrm{drive}}+H_{\mathrm{int}},$$

(1)

where H₀ is the Hamiltonian of two linear cavities, H_drive represents resonant coherent drive of the readout cavity, as shown in Figure 1a and H_int describes generation of the two identical storage photons from the readout and pump photons, as shown in Figure 1b, and the corresponding backward process, which is shown in Figure 1c. It is important to note that non-resonant coherent drive is embedded into H_int [46]. Thus, H₀ = ω_sa^†a + ω_rb^†b where ω_s and ω_r are the storage and readout frequencies, respectively, and a/a^† and b/b^† are the annihilation/creation operators corresponding to the fundamental cavity modes. The Hamiltonian of a resonant coherent drive with amplitude ε_r and frequency ω_r has the form: H_drive = ε_r⋅b^†⋅exp[-iω_rt] + ε_r^*⋅b⋅exp[iω_rt]. The last term in the Equation (1) is determined by the following expression:

$$H_{\mathrm{int}}=-\frac{1}{2}\mu (\xi e^{-i{\omega }_{p}t}b{a}^{†}{}^2+{\xi }^{*}{e}^{i{\omega }_{p}t}{b}^{†}{a}^2),$$

(2)

where μ is a nonlinear coupling constant emerging from the presence of a Josephson junction [47]. In Equation (2) pump photons are considered in the classical field approximation with effective amplitude ξ = − i⋅ε_p/[κ/2 + i(ω_r − ω_p)], where ε_p and ω_p are amplitude and frequency of the external non-resonant coherent pump [46].

Since dissipative processes play a crucial role in the considered system, its behavior must be described by Lindblad master equation [48]:

$${\partial }_t\rho =-i\left[H,\rho \right]+\kappa D[b](\rho ),$$

(3)

where ρ is a system density matrix, κ is a single-photon loss rate, D[b] is a Liouvillian, which defines as D[b](ρ) = bρb^† – ½(b^†bρ + ρb^†b). Notably, Liouvillian acts on the density matrix ρ from both sides and the perator b is called Lindblad operator or a quantum jump operator.

Further, we assume that the characteristic time scale of the single-photon dissipation (1/κ) is much smaller than all other time scales of the system. Consequently, we can eliminate degrees of freedom associated with the readout cavity [46]. Furthermore, it is convenient to use a unitary transformation U = exp[− i(ω_{p +} ω_r)a^†a/2], which makes the effective Hamiltonian time-independent. As a result, the following description of the reduced density matrix ρ_s = Tr[ρ]_r can be obtained:

$${\partial }_t{\rho }_s=-i\left[H_{\mathrm{eff}},{\rho }_s\right]+2\gamma D[a^2]({\rho }_s),$$

(4)

where the effective Hamiltonian is given by:

$$H_{\mathrm{eff}}=\Delta a^{†}a+i\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\left(g{a}^{†}{}^2-g^{*}{a}^2\right),$$

(5)

where g = 2ξ⋅μ⋅ε_r/κ is a two-photon pump rate, γ = |ξ⋅μ|²/2κ is a two-photon dissipation rate, and Δ = (2ω_s – ω_p – ω_r)/2 is a frequency detuning. The most interesting is the regime where ω_p ≈ 2ω_s − ω_r and hence Δ << ω_p, ω_s, ω_r. Further, we suppose that the two-photon pump rate g is a real value. This can be easily achieved by tuning the complex phase of the non-resonant coherent pump amplitude.

It is worth mentioning that g is proportional to the product of the non-resonant pump and resonant drive amplitudes ξ and ε_r. However, the absorption rate is γ ∝ ξ ² and does not depend on the drive amplitude ε_r. As a result, a two-photon pump and nonlinear dissipation are effectively implemented in the storage cavity due to the presence of a linear dissipation and coherent pump in the readout cavity, as shown in Figure 1d.

3. Mean-Field Theory

From now on, we will consider the behavior of the following two system observables:

$$G^K(t,t)=2n(t)+1,F^K(t,t)=2\psi (t),$$

(6)

where n(t) = Tr[a^†aρ_s(t)] is an average number of photons in the storage cavity or normal average, ψ(t) = Tr[a²ρ_s(t)] is a two-particle order parameter or an anomalous average of the system, G^K(t,t) and F^K(t,t) are the normal and anomalous simultaneous Keldysh Green’s functions. The unity term in the normal Keldysh Green’s function takes into account the presence of the QF [40,41]. As a result, one can obtain equations of motion for simultaneous Keldysh Green’s functions from the Equation (4):

$$\begin{array}{l}{\partial }_t{G}^K(t,t)=g{F}^K{(t,t)}^{*}+g^{*}{F}^K(t,t)-8\gamma \langle a^{†}{}^2{a}^2\rangle ,\\ ({\partial }_t+i2\Delta )F^K(t,t)=2g{G}^K(t,t)-4\gamma \langle a^2(2a^{†}a+1)\rangle ,\end{array}$$

(7)

where angle brackets denote the quantum mechanical averaging over the density matrix of the system. As can be seen from the Equation (7), the normal and anomalous Keldysh Green’s functions are expressed through the expectation values of higher-order operators, due to the nonlinear dissipation arising in the system. We use the mean-field approximation to avoid considering higher-order equations of motion. Within this approximation, the expectation values of the operator products are replaced by the products of their expectation values [49]. However, there are several ways to implement this decoupling [50]. The two most common choices are as follows:

$$\langle a^{†}{}^2{a}^2\rangle \approx {\psi }^{*}\psi ,$$

(8)

$$\langle a^{†}{}^2{a}^2\rangle \approx {\psi }^{*}\psi +2n^2,$$

(9)

where in (8) decoupling was carried out using only the so-called “Cooper” or “pairing” channel [51]. In (9) an additional “density” channel appears [51]. Here, we assume (8) to be valid for the considered system. Justification of this choice will be given in the following text. As a result of (8), the expectation values of the higher-order operators in the Equation (7) can be expressed through Keldysh Green’s functions as follows:

$$\begin{array}{c}\langle a^{†}{}^2{a}^2\rangle \approx \frac{1}{4}{F}^K{(t,t)}^{*}{F}^K(t,t),\\ \langle a^2(2a^{†}a+1)\rangle \approx \frac{1}{2}{F}^K(t,t)G^K(t,t),\end{array}$$

(10)

and the mean-field self-consistent equations of motion can be obtained:

$$\begin{array}{c}{\partial }_t{G}^K(t,t)=g_{\mathrm{eff}}(t)F^K{(t,t)}^{*}+g_{\mathrm{eff}}{(t)}^{*}{F}^K(t,t),\\ ({\partial }_t+i2\Delta )F^K(t,t)=2g_{\mathrm{eff}}(t)G^K(t,t),\end{array}$$

(11)

where g_eff(t) = g – 2γ⋅ψ(t) is a renormalized two-photon pump rate. We will also assume that the initial state is a vacuum, which gives the following initial conditions: G^K(0,0) = 1 and F^K(0,0) = 0. The unity term in G^K emerges from the QF, as mentioned above. This enables the parametric generation of light even if there are no photons in the initial state [52]. It should also be noted that the resulting Equation (11) have a structure similar to the Gorkov equations, which play a central role in the theory of superconductivity. Another point worth mentioning is that the (11) has the following integral of motion: G^K(t,t)² − |F^K(t,t)|². Thus, it is possible to express the normal Green’s function through the anomalous one, using vacuum initial conditions: G^K(t,t) = (1 + |F^K(t,t)|²)^1/2. As a result, one can obtain the equation of motion solely for the anomalous average:

$${\partial }_t\psi (t)=-i2\Delta \psi (t)+\sqrt{1+4{\left|\psi (t)\right|}^2}\left[g-2\gamma \psi (t)\right],$$

(12)

and relate to the average number of photons:

$$n(t)=\frac{1}{2}\left(\sqrt{1+4{\left|\psi (t)\right|}^2}-1\right).$$

(13)

We use the resulting equations of motion of the proposed mean-field theory to calculate the time evolution of the anomalous average ψ(t) and the average number of photons n(t). A comparison of the time evolution obtained from the Equations (12) and (13) and the numerical simulation of the master equation [53,54] is shown in Figure 2a,b. We perform numerical simulations using integration of a master equation on a truncated Fock basis. For this basis, equations of motion for the density operator matrix elements ρ_m,k ≡ <m|ρ_s|k>, where |k> is a Fock state, will have the following form:

$$\begin{array}{l}{\partial }_t{\rho }_m{}_{,k}=-i\Delta (m-k){\rho }_m{}_{,k}+\frac{g}{2}\left[\sqrt{m(m-1)}{\rho }_{m-2}{}_{,k}-\sqrt{(k+1)(k+2)}{\rho }_m{}_{,k+2}\right]-\\ -\frac{g^{*}}{2}\left[\sqrt{(m+1)(m+2)}{\rho }_{m+2}{}_{,k}-\sqrt{k(k-1)}{\rho }_m{}_{,k-2}\right]-\gamma [\left(k(k-1)+m(m-1)\right){\rho }_m{}_{,k}-\\ -2\sqrt{(k+1)(k+2)(m+1)(m+2)}{\rho }_{m+2}{}_{,k+2}].\end{array}$$

(14)

Since a vacuum is assumed as the initial quantum state, the initial condition for the matrix elements was chosen as ρ_m,k(t = 0) = δ_0,kδ_m,k. Here δ_m,k is the Kronecker delta function, which leads to ρ_0,0 = 1 and ρ_m,k = 0 for all other elements. As long as (14) is a system of linear differential equations, it can be solved numerically using the fourth-order Runge-Kutta method. Photon-number cutoff is chosen in accordance with the the value of two-photon pump rate and is justified by checking the low probabilities for large photon numbers. To make sure that the cutoff error is negligible, we increased the cutoff number and checked the convergence of the solution.

Figure 2 shows that the mean-field theory qualitatively describes well the evolution of calculated system observables. As pointed above, the influence of QF is tremendous in a near-zero-time region by analogy with the parametric generation of light from a vacuum state. As can be seen from Figure 2, the semiclassical solution is exact zero, since generation cannot start without QF. At a longer time scale, the presence of QF significantly affects the dynamics of the system observables, which leads to the formation of a stationary state.

For the case of zero frequency detuning, the following mean-field stationary solutions can be found from Equation (12) for the anomalous average and the average number of photons in the storage cavity:

$$\psi =g/2\gamma ,n=\frac{1}{2}\left(\sqrt{1+4{\psi }^2}-1\right).$$

(15)

In addition, using the exact stationary solution for the density matrix for Δ = 0 [55], one can find ψ and n in the explicit form:

$$\psi =g/2\gamma ,n=\psi \tanh \left(\psi \right)$$

(16)

A comparison of the stationary solutions calculated by the mean-field theory (15) and the exact stationary density matrix (16) demonstrates very good agreement of the anomalous averages, as shown in Figure 2a. However, it can be seen that the average number of photons differs between the mean-field and exact solutions (Figure 2b). Nevertheless, both solutions have the same asymptotic behavior in the limit of a small (g << γ) and large (g >> γ) two-photon pump rate.

4. Properties of the Steady State

For the case of non-zero frequency detuning, we derive the following biquadratic equation on the modulus of the anomalous average:

$${\Delta }^2{\left|\psi \right|}^2=\frac{1}{4}(1+4{\left|\psi \right|}^2)(g^2-4{\gamma }^2{\left|\psi \right|}^2),$$

(17)

where the unity term corresponds to the presence of QF. Dropping QF in Equation (17), one gets a semiclassical solution [33]:

$$\left|\psi \right|\approx \theta (g^2-{\Delta }^2)\sqrt{g^2-{\Delta }^2}/2\gamma ,$$

(18)

where θ(x) is a Heaviside step function. Equation (18) yields non-zero anomalous average ψ when the frequency detuning is below than the two-photon pump rate (Δ < g), and zero ψ value after passing the critical point Δ = g (green curve in Figure 3b). For open quantum systems, this phenomenon is also known as DPT [3,33,56,57], and was studied for a considered system earlier [2,33]. Here, we can observe DPT in the “thermodynamic limit”, when the average number of photons n or the anomalous average ψ tends to infinity, and one can usually ignore QF [57]. For our system, this limit can be realized in the case of a large two-photon pump rate (g >> γ).

We will consider the situation when the average number of photons n and the anomalous average ψ are the finite quantities and it is absolutely necessary to take into account the effects of QF. Within the framework of Keldysh formalism, this leads to the following modification of the semiclassical solution (18):

$$\left|\psi \right|={\left(\sqrt{g^2{\gamma }^2+\frac{1}{4}{\left(g^2-{\Delta }^2-{\gamma }^2\right)}^2}+\frac{1}{2}\left(g^2-{\Delta }^2-{\gamma }^2\right)\right)}^{1/2}/2\gamma .$$

(19)

As can be seen from Figure 3a–c, the presence of QF makes the DPT wider near the critical point. This fact is well described by the proposed mean-field theory. To obtain a simple physical explanation for this broadening, we shall demonstrate an analogy with Landau theory of phase transitions [58]. For that, it is necessary to rewrite the Equation (17), sorting the contributions in descending order of powers of ψ:

$$4{\gamma }^2{\left|\psi \right|}^4+({\Delta }^2-g^2+{\gamma }^2){\left|\psi \right|}^2=\frac{1}{4}{g}^2,$$

(20)

$$4B\cdot {\eta }^3+2A\cdot \eta =h,$$

(21)

where η is an order parameter, h is an external field, A and B are series expansion coefficients of the thermodynamic potential Φ = Aη² + Bη⁴ − hη. Minimization of the Φ yields the expression in (21). Comparing (20) and (21), we can conclude that the presence of QF leads to the analog of the external field h. The presence of such an external field in Landau theory breaks the symmetry of the system. As a result, the difference between the two phases disappears, as well as the discrete phase transition point [58]. Therefore, phase transition itself is broadening, which is observed in our system.

Also, it is well known that at the phase transition point defined by a condition A = 0, the response to the external field is nonlinear and determined by the power-law η∝(h/B)^1/δ [58]. Here δ is one of the critical exponents for the nonzero external field and it is equal to δ = 3 in Landau theory. Thus, a similar phenomenon should be observed in our system as well. To confirm this statement, we examine the behavior of the anomalous average ψ as a function of the pump rate g at the critical point Δ = g. Here we are interested in the regime g >> γ. Figure 3d demonstrates a comparison between the predictions of the mean-field theory and the exact numerical simulation, which is governed by the following power-law |ψ|∝(g/γ)^1/δ and the corresponding critical exponent: δ = 1.47 for the numerical solution and δ = 2 for the mean-field theory.

It is important to note that the experimental verification of the predicted power-law is very feasible within the framework of the discussed setup (Figure 1). The frequency detuning Δ can be controlled by changing the non-resonant pump frequency ω_p, and the two-photon pump rate g by increasing the amplitude ε_p of the external non-resonant coherent pump wave.

Second non-trivial result associated with QF in the proposed mean-field theory is the renormalization of the critical point. Setting the coefficient before |ψ|² in (20) to zero yields the following expression for the case of g ⩾ γ:

$${\Delta }_0(g)=\sqrt{g^2-{\gamma }^2}.$$

(22)

The dependence of the critical point on a pump rate g obtained from the mean-field theory (22) is shown in Figure 4a. As one can see, the mean-field theory predicts the existence of a threshold pump rate, which is equal to g_th = γ.

To understand the physical meaning of the critical point renormalization, it is necessary to consider the case of a small pump rate (g << γ). In this limit, the anomalous average ψ is quite small and one can neglect the quartic part in the Equation (20). Thus, the response of the system to the fluctuation field is determined only by a quadratic contribution. The crucial fact here is that the coefficient before |ψ|² is a positive quantity for arbitrary frequency detuning Δ and g < γ. In Landau theory, such behavior corresponds to the symmetric phase of the system (A > 0). Further increasing the pump rate g would decrease the corresponding coefficient (Δ² − g² + γ²). As result, it becomes zero at the threshold point g_th = γ and Δ = 0, which violates the linear response and makes the quartic coefficient dominant. This transformation of the system’s response correlates precisely with the emergence of the discussed critical behavior.

One can observe such a transition from the analysis of the second derivative ∂²|ψ|/∂Δ² in the region Δ ≈ 0, which is shown in Figure 4b. It can be seen that its behavior changes from the linear asymptotic ∝g to the inverse proportion ∝ 1/g after passing the threshold point g_th = γ. Furthermore, considering the dynamics of the system’s susceptibility χ = ∂|ψ|/∂g as a function of the pump rate g and frequency detuning Δ, we also observe the existence of a threshold pump g_th. Its behavior in the framework of the mean-field theory is shown in Figure 4c. One can see that for g < γ the maximum susceptibility is at the point Δ_max = 0. However, after passing the threshold point g_th = γ one can observe a shift of the maximum, which indicates the emergence of the critical behavior. The numerical simulations of the system’s susceptibility shown in Figure 4d also predict the existence of g_th ≈ 1.952γ, which is in qualitative agreement with the proposed Keldysh formalism.

Finally, justification and constraints of our mean-field treatment should be discussed. It is well known that the DPT in the considered system is associated with transformation from the Schrödinger cat to the squeezed-vacuum state [33] as shown in Figure 3b. For both states, the anomalous average ψ ≠ 0, but they differ in the decoupling of channels, as shown in Equations (8) and (9). As it turns out, the proposed decoupling is rigorous only for a large number of photons region, which corresponds to the cat state. However, in the region of a small number of photons, corresponding to the squeezed state, the anomalous average |ψ| << 1. From Equation (13) follows that: |ψ| >> n. Consequently, in the region of a small number of photons, the “Cooper” channel dominates over the “density” channel and the proposed decoupling is valid. Also, we have shown above that and near the critical point, the proposed description allows us to achieve a good qualitative agreement with the exact results of numerical calculations. We believe that the observed mismatch can be explained by the presence of a fluctuation region, where the QF significantly affects the behavior of the system. Thus, for further exploration, it seems to be necessary to apply more precise methods, such as the Keldysh functional integral [41], the 2PI effective action [59], and the renormalization group [60,61].

5. Conclusions and Outlook

In this paper, we develop a mean-field theory for a system with two-photon driving and dissipation which explicitly takes quantum fluctuations into account. Consideration of quantum fluctuations allows us to describe properly the evolution dynamics of the system and to demonstrate the new effects in the steady-state. We show that the dissipative phase transition broadening near the critical point is naturally conditioned by the fluctuation field. Its presence leads to a power-law dependence of the anomalous average at the phase transition point with which the critical exponent is associated. This counterintuitive effect cannot be found in the semiclassical approximation. The reason lies in the crucial role of the quantum fluctuations, which significantly affect the behavior of a given system at the critical point and should be taken into account even in the case of a large average number of photons. Also, we investigate the effect of the quantum fluctuations on the critical point renormalization and demonstrate the existence of a two-photon pump “threshold”. It is noteworthy that the obtained results are in a good agreement with numerical simulations. An important point here is that intuitively free energy analysis is out of consideration for dissipative phase transitions. However, obtained counter-intuitive results are largely based on analogies with Landau theory. We believe that further investigations will clarify this close relationship between purely non-equilibrium formalism and traditional equilibrium techniques. Such a relationship is most likely the result of Boltzmann-Gibbs-like stationary distribution created by the quantum Langevin equation, which describes the additional QF. Thus, the rigorous justification of the proposed equilibrium description based on the thermodynamic potential is worth extensive study.

The results presented in this paper can be applied to the development of the new nonequilibrium quantum states with controllable properties [62,63]. They may also be beneficial for applications in quantum information processing [64] and quantum metrology [65], as we believe that the quantum photonic state and systems spectral characteristics in the critical point have a non-conventional properties due to the observed nontrivial scaling behavior. Investigated quantum systems, in combination with terahertz to gigahertz light convertors [66], can be included into fully quantum photonic networks. This strategy opens new avenues for utilization of important and promising terahertz light sources such as quantum cascade lasers. These compact and powerful sources of light can be used for optical drive of the studied finite-component system consisting of a several bosonic modes interacting through a Josephson junction. In addition, quantum cascade lasers can be very convenient as a source of broadband continuous single-mode tuning [67] due to the experimental need to change the frequency detuning discussed above.