Introduction
Quantum sensing constitutes a prominent research frontier at the interface of quantum science and precision metrology. Its goal is to convert coherence, squeezing, entanglement, and quantum manipulation into high-precision estimations of weak physical quantities under the constraints of quantum noise, measurement feedback, and resource limitations. In recent years, quantum sensing has not only advanced fundamental physical tests but also provided key technical support for the development of quantum communication, quantum computing, and low-power quantum devices [4 , 5 , 6 , 7 , 8 , 9 ] [1, 2, 3]. Among various physical carriers, resonators have become an important platform for precise measurement and weak signal detection due to their adjustable frequency, high quality factor, and good compatibility with superconducting circuits, optical mechanical systems, and cavity quantum electrodynamics platforms [9 , 10 ] . Therefore, how to measure a weak driving signal that is much smaller than the thermal noise or detection noise level is one of the core issues in resonator quantum sensing research.
Realistic quantum systems inevitably couple to thermal reservoirs, measurement channels and control lines, giving rise to ubiquitous dissipation. Most studies adopt linear damping and Markov approximation, using the Lindblad master equation to describe the energy exchange of the system with the environment [11 ] . However, many practical resonators exhibit amplitude-dependent damping: the system gains energy at small amplitudes and is suppressed at large amplitudes due to nonlinear mechanisms, eventually forming stable self-sustained oscillations. This behavior is typically described by the vdP model, whose classical form includes the competition between negative linear damping and cubic nonlinear damping [12 ] . When approaching the self-sustained oscillation threshold or the low amplitude limit, the response of the classical vdP oscillator can become very steep, and weak drive can even lead to sensitivity divergence. This nonlinear amplification mechanism is believed to be related to the high-sensitivity active amplification of biological mechanical receptors such as cochlear hair bundle cells, enabling biological systems to detect extremely weak mechanical stimuli in strong background noise [13 , 14 ] . After quantizing the vdP oscillator, zero-point fluctuations, transitions to the heat reservoir, quantum measurement feedback, and finite excitation numbers introduce noise that is absent in the classical model. The steady state of the quantum vdP oscillator is no longer completely described by a single deterministic trajectory; its amplitude diffusion, phase diffusion, and synchronization behaviors need to be studied using quantum master equations and phase space distributions, etc. Existing theoretical studies have analyzed the locking and synchronization thresholds of quantum vdP oscillators using quantities such as quantum phase space distributions and two-oscillator correlations [15 , 16 ] . These results indicate that the zero-drive sensitivity divergence observed in the classical model cannot be directly extrapolated to the quantum regime: quantum noise changes the response scale near the threshold and creates a new trade-off between ultra-high sensitivity and measurement stability. When the investigation extends from a single oscillator to two coupled oscillators, information transfer pathways and noise backflow emerge as critical factors governing sensing performance. Non-reciprocal coupling can establish directional information flow between coupled resonators or cavity systems, allowing signals to be transmitted from one oscillator to another while suppressing disturbances and measurement feedback in the opposite direction. This provides the possibility of constructing directional weak signal sensors and unidirectional quantum transducers [17 , 18 ] . Meanwhile, studies of two-oscillator quantum synchronization [19 , 20 , 21 , 22 ] focus on whether two self-sustaining oscillators can establish a stable phase relationship and how this phase relationship is affected by quantum noise and dissipation. These two lines of research respectively indicate that non-reciprocal coupling facilitates directional information transmission and that the quantum vdP oscillator serves as an important platform for studying quantum synchronization. One key question arising from this is: When two quantum vdP oscillators are connected through non-reciprocal coupling, can the non-reciprocal properties improve the quantum sensing precision of weak driving signals?
In this article, we employ two vdP oscillators to implement quantum sensing measurements. A weak coherent drive acts on one vdP resonator. Through coherent and dissipative coupling, the information of the drive is transferred to the other resonator. The amplitude information of the drive is then extracted via homodyne measurement. When the incoherent drive and single-particle loss are both smaller than the two-particle loss, the measurement precision is limited and will not diverge. Complete non-reciprocal coupling will result in a measurement precision that is twice as high as that obtained with reciprocal coupling. When the incoherent drive and dissipation, as well as the coherent coupling strength are all equal to 0, the dark mode forms under zero coherent drive. This is beneficial for the measurement of extremely weak drive, ultimately achieving the Heisenberg limit of measurement precision in the time domain. Under strong incoherent drive, the measurement precision increases as the two-particle loss rate decreases, and the factor by which the non-reciprocal coupling exceeds the reciprocal coupling also increases as the two-particle loss rate decreases. At the critical incoherent driving condition, only when the weak driving force is near zero, can the measurement uncertainty obtained with non-reciprocal coupling tend to 0 as the two-particle loss rate approaches 0. At this point, compared with reciprocal coupling, the improvement in measurement precision obtained through non-reciprocal coupling will also increase as the two-particle loss rate decreases.
The rest of this article is organized as follows. In section II, the combination of coherent and dissipative couplings realizes the nonreciprocal coupling in the cascaded vdP oscillators. In section III, we obtain finite measurement precision by employing weak incoherent drive for the quantum vdP resonator. When only dissipation coupling is present, a dark mode appears at zero drive in section IV. In section V, measurement precision is obtained under strong incoherent driving conditions. In section VI, the sensing precision for weak drive at critical incoherent drive is investigated. We make a conclusion in section VII.
Nonreciprocal coupling in vdP oscillators
We consider two resonant quantum vdP oscillators coupled together via coherent and dissipative couplings, whose corresponding dynamical equations can be written in the form of a Lindblad master equation under the Markov approximation. In the rotating frame, the density matrix is described by [18 , 23 ]
ρ ˙ = − i [ H d + H c , ρ ] + γ 1 + D a † [ ρ ] + γ 1 − D a [ ρ ] + γ 2 D a 2 [ ρ ] + Γ 1 + D b † [ ρ ] + Γ 1 − D b [ ρ ] + Γ 2 D b 2 [ ρ ] + G D a + b [ ρ ] (1) \dot{\rho}=-i[H_d+H_c,\rho]+\gamma_1^+\mathcal{D}_{a^\dagger}[\rho]+\gamma_1^-\mathcal{D}_{a}[\rho]+\gamma_2\mathcal{D}_{a^2}[\rho]+\Gamma_1^+\mathcal{D}_{b^\dagger}[\rho]+\Gamma_1^-\mathcal{D}_{b}[\rho]+\Gamma_2\mathcal{D}_{b^2}[\rho]+\mathcal{G}\mathcal{D}_{a+b}[\rho] \tag*{(1)} ρ ˙ = − i [ H d + H c , ρ ] + γ 1 + D a † [ ρ ] + γ 1 − D a [ ρ ] + γ 2 D a 2 [ ρ ] + Γ 1 + D b † [ ρ ] + Γ 1 − D b [ ρ ] + Γ 2 D b 2 [ ρ ] + G D a + b [ ρ ] (1) where a a a b b b H d = i Ω ( a † − a ) H_d=i\Omega(a^\dagger-a) H d = i Ω ( a † − a ) Ω \Omega Ω D c [ ρ ] = 2 c ρ c † − c † c ρ − ρ c † c \mathcal{D}_{c}[\rho]=2c\rho c^\dagger-c^\dagger c\rho-\rho c^\dagger c D c [ ρ ] = 2 c ρ c † − c † c ρ − ρ c † c γ 1 + \gamma_1^+ γ 1 + Γ 1 + \Gamma_1^+ Γ 1 + γ 1 − \gamma_1^- γ 1 − Γ 1 − \Gamma_1^- Γ 1 − γ 2 \gamma_2 γ 2 Γ 2 \Gamma_2 Γ 2 G G G H c H_c H c
H c = i g ( a † b − a b † ) (2) H_c=ig(a^\dagger b-ab^\dagger) \tag*{(2)} H c = i g ( a † b − a b † ) (2) where g g g
According to the Lindblad master equation, the expectation values of the operators ⟨ a ⟩ \langle a\rangle ⟨ a ⟩ ⟨ b ⟩ \langle b\rangle ⟨ b ⟩
⟨ a ˙ ⟩ = ( γ 1 + − γ 1 − − G ) ⟨ a ⟩ + ( g − G ) ⟨ b ⟩ − 2 γ 2 ⟨ a † a 2 ⟩ + Ω (3) \langle\dot{a}\rangle=(\gamma_1^+-\gamma_1^--G)\langle a\rangle+(g-G)\langle b\rangle-2\gamma_2\langle a^\dagger a^2\rangle+\Omega \tag*{(3)} ⟨ a ˙ ⟩ = ( γ 1 + − γ 1 − − G ) ⟨ a ⟩ + ( g − G ) ⟨ b ⟩ − 2 γ 2 ⟨ a † a 2 ⟩ + Ω (3) ⟨ b ˙ ⟩ = ( Γ 1 + − Γ 1 − − G ) ⟨ b ⟩ − ( g + G ) ⟨ a ⟩ − 2 Γ 2 ⟨ b † b 2 ⟩ . (4) \langle\dot{b}\rangle=(\Gamma_1^+-\Gamma_1^--G)\langle b\rangle-(g+G)\langle a\rangle-2\Gamma_2\langle b^\dagger b^2\rangle. \tag*{(4)} ⟨ b ˙ ⟩ = ( Γ 1 + − Γ 1 − − G ) ⟨ b ⟩ − ( g + G ) ⟨ a ⟩ − 2 Γ 2 ⟨ b † b 2 ⟩ . (4) When both dissipative and coherent couplings are present simultaneously, the effective coupling between the two systems becomes non-reciprocal, i.e., g − G ≠ g + G g-G\ne g+G g − G = g + G G = g G=g G = g a a a b b b b b b a a a
Replacing the operator a a a b b b α \alpha α β \beta β γ 1 + − γ 1 − − G = 0 \gamma_1^+-\gamma_1^--G=0 γ 1 + − γ 1 − − G = 0 Γ 1 + − Γ 1 − − G = 0 \Gamma_1^+-\Gamma_1^--G=0 Γ 1 + − Γ 1 − − G = 0 β \beta β Ω → 0 \Omega\to0 Ω → 0
β = g Γ 2 ( Ω 2 γ 2 ) 1 / 9 , G = g ; (5) \beta=\frac{g}{\Gamma_2}\left(\frac{\Omega}{2\gamma_2}\right)^{1/9},\ G=g; \tag*{(5)} β = Γ 2 g ( 2 γ 2 Ω ) 1/9 , G = g ; (5) β = Ω G − g , ∣ G − g ∣ ≫ Ω . (6) \beta=\frac{\Omega}{G-g},\ |G-g|\gg\Omega. \tag*{(6)} β = G − g Ω , ∣ G − g ∣ ≫ Ω. (6) The susceptibility χ c = d α d Ω \chi_c=\frac{d\alpha}{d\Omega} χ c = d Ω d α
χ c = − g Γ 2 1 9 ( 1 2 γ 2 Ω 8 ) 1 / 9 , G = g ; (7) \chi_c=-\frac{g}{\Gamma_2}\frac{1}{9}\left(\frac{1}{2\gamma_2\Omega^8}\right)^{1/9},\ G=g; \tag*{(7)} χ c = − Γ 2 g 9 1 ( 2 γ 2 Ω 8 1 ) 1/9 , G = g ; (7) χ c = 1 G − g , ∣ G − g ∣ ≫ Ω . (8) \chi_c=\frac{1}{G-g},\ |G-g|\gg\Omega. \tag*{(8)} χ c = G − g 1 , ∣ G − g ∣ ≫ Ω. (8) It shows that the susceptibility is divergent at zero drive χ c ∝ Ω − 8 / 9 \chi_c\propto\Omega^{-8/9} χ c ∝ Ω − 8/9 χ c ∝ Ω − 2 / 3 \chi_c \propto\Omega^{-2/3} χ c ∝ Ω − 2/3 [24 ] . This indicates that complete non-reciprocal coupling can significantly enhance the measurement susceptibility of weak driving signals in classical vdP oscillators.
Quantum vdP Oscillators with Weak Incoherent Driving
In this section, we investigate whether the quantum vdP oscillator still exhibits divergent sensitivity at the zero coherent drive as Ω → 0 \Omega\to0 Ω → 0
When the two-particle loss rates are large enough, i.e., γ 1 + ≪ γ 2 \gamma_1^+ \ll\gamma_2 γ 1 + ≪ γ 2 Γ 1 + ≪ Γ 2 \Gamma_1^+ \ll\Gamma_2 Γ 1 + ≪ Γ 2 Ω ≪ 1 \Omega\ll1 Ω ≪ 1 [25 , 26 ] and truncating to the lowest three energy levels for each subsystem, the expectation value of b b b Figure 1 , the analytical results obtained via the non-Hermitian method and the numerical results from the master equation become increasingly consistent as Ω \Omega Ω O ( Ω ) O(\Omega) O ( Ω )
Figure 1. The comparison graph of the b b b γ 1 + = Γ 1 + = 0.25 \gamma_1^+=\Gamma_1^+=0.25 γ 1 + = Γ 1 + = 0.25 γ 1 − = Γ 1 − = 0.05 \gamma_1^-=\Gamma_1^-=0.05 γ 1 − = Γ 1 − = 0.05 γ 2 = Γ 2 = 5.0 \gamma_2=\Gamma_2=5.0 γ 2 = Γ 2 = 5.0 g = G = 0.02 g=G=0.02 g = G = 0.02
⟨ a ⟩ = β Ω Δ 1 + O ( Ω 3 ) , (9) \langle a\rangle= \frac{\beta\Omega}{\Delta_1} + O(\Omega^3), \tag*{(9)} ⟨ a ⟩ = Δ 1 β Ω + O ( Ω 3 ) , (9) ⟨ b ⟩ = − ( g + G ) Ω Δ 1 + O ( Ω 3 ) , (10) \langle b\rangle= -\frac{(g+G)\Omega}{\Delta_1} + O(\Omega^3), \tag*{(10)} ⟨ b ⟩ = − Δ 1 ( g + G ) Ω + O ( Ω 3 ) , (10) ⟨ b † b ⟩ = − ( g + G ) 2 Ω 2 g 4 + O ( Ω 3 ) , (11) \langle b^\dagger b\rangle= -\frac{(g+G)^2\Omega^2}{g^4} + O(\Omega^3), \tag*{(11)} ⟨ b † b ⟩ = − g 4 ( g + G ) 2 Ω 2 + O ( Ω 3 ) , (11) where Δ 1 = μ ν + g 2 − G 2 \Delta_1 = \mu\nu+ g^2 - G^2 Δ 1 = μν + g 2 − G 2 μ = γ 1 + + γ 1 − + G \mu= \gamma_1^+ + \gamma_1^- + G μ = γ 1 + + γ 1 − + G ν = Γ 1 + + Γ 1 − + G \nu= \Gamma_1^+ + \Gamma_1^- + G ν = Γ 1 + + Γ 1 − + G
When g 2 + ( γ 1 + + γ 1 − ) ( Γ 1 + + Γ 1 − ) ≫ g ( γ 1 + + γ 1 − + Γ 1 + + Γ 1 − ) g^2+(\gamma_1^+ + \gamma_1^-)(\Gamma_1^+ + \Gamma_1^-) \gg g(\gamma_1^+ + \gamma_1^- + \Gamma_1^+ + \Gamma_1^-) g 2 + ( γ 1 + + γ 1 − ) ( Γ 1 + + Γ 1 − ) ≫ g ( γ 1 + + γ 1 − + Γ 1 + + Γ 1 − ) χ G = g / χ G = 0 = 2 \chi_{G=g}/\chi_{G=0}=2 χ G = g / χ G = 0 = 2 M = ( b † + b ) M=(b^\dagger+b) M = ( b † + b )
Δ M = ⟨ ( b † + b ) 2 ⟩ − ⟨ ( b † + b ) ⟩ 2 = 1 + O ( Ω 2 ) . (12) \Delta M=\langle(b^\dagger+b)^2\rangle-\langle(b^\dagger+b)\rangle^2=1+O(\Omega^2). \tag*{(12)} Δ M = ⟨( b † + b ) 2 ⟩ − ⟨( b † + b ) ⟩ 2 = 1 + O ( Ω 2 ) . (12) Then, we obtain the same sensitivity factor as that of the susceptibility
δ Ω G = g / δ Ω G = 0 = 2. (13) \delta\Omega_{G=g}/\delta\Omega_{G=0}=2. \tag*{(13)} δ Ω G = g / δ Ω G = 0 = 2. (13) Steady-State Dark-Mode Approximation
When γ 1 + = γ 1 − = 0 \gamma_1^+=\gamma_1^-=0 γ 1 + = γ 1 − = 0 Γ 1 + = Γ 1 − = 0 \Gamma_1^+=\Gamma_1^-=0 Γ 1 + = Γ 1 − = 0 g = 0 g=0 g = 0 G > 0 G>0 G > 0 c = a + b 2 c=\frac{a+b}{\sqrt{2}} c = 2 a + b d = a − b 2 d=\frac{a-b}{\sqrt{2}} d = 2 a − b β ≃ − α \beta\simeq-\alpha β ≃ − α
We obtain the first-moment equations
d d t ⟨ a ⟩ = Ω − G ⟨ a + b ⟩ − 2 γ 2 ⟨ a † a 2 ⟩ , (14) \frac{d}{dt}\langle a\rangle=\Omega-G\langle a+b\rangle-2\gamma_2\langle a^\dagger a^2\rangle, \tag*{(14)} d t d ⟨ a ⟩ = Ω − G ⟨ a + b ⟩ − 2 γ 2 ⟨ a † a 2 ⟩ , (14) d d t ⟨ b ⟩ = − G ⟨ a + b ⟩ − 2 γ 2 ⟨ b † b 2 ⟩ . (15) \frac{d}{dt}\langle b\rangle=-G\langle a+b\rangle-2\gamma_2\langle b^\dagger b^2\rangle. \tag*{(15)} d t d ⟨ b ⟩ = − G ⟨ a + b ⟩ − 2 γ 2 ⟨ b † b 2 ⟩ . (15) Under the factorization approximation
⟨ a † a 2 ⟩ ≃ ∣ α ∣ 2 α , ⟨ b † b 2 ⟩ ≃ ∣ β ∣ 2 β , (16) \langle a^\dagger a^2\rangle\simeq|\alpha|^2\alpha,\qquad\langle b^\dagger b^2\rangle\simeq|\beta|^2\beta, \tag*{(16)} ⟨ a † a 2 ⟩ ≃ ∣ α ∣ 2 α , ⟨ b † b 2 ⟩ ≃ ∣ β ∣ 2 β , (16) the semiclassical equations become
α ˙ = Ω − G ( α + β ) − 2 γ 2 ∣ α ∣ 2 α , (17) \dot{\alpha}=\Omega-G(\alpha+\beta)-2\gamma_2|\alpha|^2\alpha, \tag*{(17)} α ˙ = Ω − G ( α + β ) − 2 γ 2 ∣ α ∣ 2 α , (17) β ˙ = − G ( α + β ) − 2 γ 2 ∣ β ∣ 2 β . (18) \dot{\beta}=-G(\alpha+\beta)-2\gamma_2|\beta|^2\beta. \tag*{(18)} β ˙ = − G ( α + β ) − 2 γ 2 ∣ β ∣ 2 β . (18) Assuming the steady-state amplitudes are real, we set
α = s + q , β = − s + q , (19) \alpha=s+q,\qquad\beta=-s+q, \tag*{(19)} α = s + q , β = − s + q , (19) where s s s q q q
α + β = 2 q . (20) \alpha+\beta=2q. \tag*{(20)} α + β = 2 q . (20) The steady-state equations are
Ω − 2 G q − 2 γ 2 ( s + q ) 3 = 0 , (21) \Omega-2Gq-2\gamma_2(s+q)^3=0, \tag*{(21)} Ω − 2 Gq − 2 γ 2 ( s + q ) 3 = 0 , (21) 2 G q + 2 γ 2 ( − s + q ) 3 = 0. (22) 2Gq+2\gamma_2(-s+q)^3=0. \tag*{(22)} 2 Gq + 2 γ 2 ( − s + q ) 3 = 0. (22) Subtracting the two equations gives
Ω − 4 γ 2 s 3 − 12 γ 2 s q 2 = 0. (23) \Omega-4\gamma_2s^3-12\gamma_2sq^2=0. \tag*{(23)} Ω − 4 γ 2 s 3 − 12 γ 2 s q 2 = 0. (23) Adding the two equations gives
Ω − 4 G q − 12 γ 2 s 2 q − 4 γ 2 q 3 = 0. (24) \Omega-4Gq-12\gamma_2s^2q-4\gamma_2q^3=0. \tag*{(24)} Ω − 4 Gq − 12 γ 2 s 2 q − 4 γ 2 q 3 = 0. (24) For fixed G > 0 G > 0 G > 0 Ω → 0 \Omega\to0 Ω → 0
s = ( Ω 4 γ 2 ) 1 / 3 , q = Ω 4 G + O ( Ω 5 / 3 , Ω 3 ) . (25) s = \left(\frac{\Omega}{4\gamma_2}\right)^{1/3},\quad q = \frac{\Omega}{4G} + O(\Omega^{5/3},\Omega^3). \tag*{(25)} s = ( 4 γ 2 Ω ) 1/3 , q = 4 G Ω + O ( Ω 5/3 , Ω 3 ) . (25) Therefore, the steady-state amplitudes are
⟨ a ⟩ s s = ( Ω 4 γ 2 ) 1 / 3 + Ω 4 G + ⋯ , (26) \langle a\rangle_{\mathrm{ss}} = \left(\frac{\Omega}{4\gamma_2}\right)^{1/3} + \frac{\Omega}{4G} + \cdots, \tag*{(26)} ⟨ a ⟩ ss = ( 4 γ 2 Ω ) 1/3 + 4 G Ω + ⋯ , (26) ⟨ b ⟩ s s = − ( Ω 4 γ 2 ) 1 / 3 + Ω 4 G + ⋯ . (27) \langle b\rangle_{\mathrm{ss}} = -\left(\frac{\Omega}{4\gamma_2}\right)^{1/3} + \frac{\Omega}{4G} + \cdots. \tag*{(27)} ⟨ b ⟩ ss = − ( 4 γ 2 Ω ) 1/3 + 4 G Ω + ⋯ . (27) According to the error propagation formula, the measurement uncertainty of Ω \Omega Ω M = b + b † M = b + b^{\dagger} M = b + b †
δ Ω = ⟨ M 2 ⟩ − ∣ ⟨ M ⟩ ∣ 2 ∣ ∂ Ω M ∣ = 6 ( 4 Ω 2 γ 2 ) 1 / 3 . (28) \delta\Omega= \frac{\sqrt{\langle M^2\rangle- \lvert\langle M\rangle\rvert^2}}{\lvert\partial_{\Omega}M\rvert} = 6(4\Omega^2\gamma_2)^{1/3}. \tag*{(28)} δ Ω = ∣ ∂ Ω M ∣ ⟨ M 2 ⟩ − ∣⟨ M ⟩ ∣ 2 = 6 ( 4 Ω 2 γ 2 ) 1/3 . (28) The characteristic time for reaching the steady state is determined by the eigenvalue with the smallest absolute value of its real part. (see the detail in Appendix § 6.2 ):
T c = ( γ 2 Ω 2 2 ) − 1 / 3 . (29) T_c = \left(\frac{\gamma_2\Omega^2}{2}\right)^{-1/3}. \tag*{(29)} T c = ( 2 γ 2 Ω 2 ) − 1/3 . (29) The measurement uncertainty of Ω \Omega Ω
δ Ω = 12 / T c . (30) \delta\Omega= 12/T_c. \tag*{(30)} δ Ω = 12/ T c . (30) Therefore, as Ω \Omega Ω [27 ] is achieved in the limit of vanishing Ω \Omega Ω Ω → 0 \Omega\to0 Ω → 0
Large incoherent drive
We assume the incoherent drive γ 1 + \gamma_1^+ γ 1 + γ 1 − \gamma_1^- γ 1 − γ 1 + − γ 1 − − G ≫ γ 2 \gamma_1^+ - \gamma_1^- - G \gg\gamma_2 γ 1 + − γ 1 − − G ≫ γ 2 Γ 1 + − Γ 1 − − G ≫ Γ 2 \Gamma_1^+ - \Gamma_1^- - G \gg\Gamma_2 Γ 1 + − Γ 1 − − G ≫ Γ 2 [28 , 29 , 30 , 31 ] is used in the large-occupation regime. As can be seen from Figure 2 , the result obtained by this approximation is close to the exact numerical solution when Γ 2 \Gamma_2 Γ 2 § 11 )
Figure 2. Comparison diagram of using Kramers-Moyal approximation and numerical calculation in the case of non-reciprocal coupling. Here, the dimensionless parameters are chosen as: γ 1 + = 1 \gamma_1^+ = 1 γ 1 + = 1 Γ 1 + = 0.70 \Gamma_1^+ = 0.70 Γ 1 + = 0.70 γ 1 − = 0.2 \gamma_1^- = 0.2 γ 1 − = 0.2 Γ 1 − = 0.15 \Gamma_1^- = 0.15 Γ 1 − = 0.15 γ 2 = 0.2 \gamma_2 = 0.2 γ 2 = 0.2 Ω = 0.001 \Omega= 0.001 Ω = 0.001 g = G = 0.05 g = G = 0.05 g = G = 0.05
⟨ b † b ⟩ s s = Γ 1 + − Γ 1 − − G 2 Γ 2 + O ( g 2 , G 2 , g G , Ω 2 ) , (31) \langle b^{\dagger}b\rangle_{\mathrm{ss}} = \frac{\Gamma_1^+ - \Gamma_1^- - G}{2\Gamma_2} + O(g^2,G^2,gG,\Omega^2), \tag*{(31)} ⟨ b † b ⟩ ss = 2 Γ 2 Γ 1 + − Γ 1 − − G + O ( g 2 , G 2 , g G , Ω 2 ) , (31) ⟨ b ⟩ s s = − ( g + G ) ( γ 1 + − γ 1 − − G ) ( Γ 1 + − Γ 1 − − G ) Ω γ 2 Γ 2 ( 3 γ 1 + − γ 1 − − G ) ( 3 Γ 1 + − Γ 1 − − G ) [ 1 + ( g 2 − G 2 ) ( γ 1 + − γ 1 − − G ) ( Γ 1 + − Γ 1 − − G ) γ 2 Γ 2 ( 3 γ 1 + − γ 1 − − G ) ( 3 Γ 1 + − Γ 1 − − G ) ] − 1 . (32) \langle b\rangle_{\mathrm{ss}} = -\frac{(g+G)(\gamma_1^+ - \gamma_1^- - G)(\Gamma_1^+ - \Gamma_1^- - G)\Omega}{\gamma_2\Gamma_2(3\gamma_1^+ - \gamma_1^- - G)(3\Gamma_1^+ - \Gamma_1^- - G)}
\left[
1+\frac{(g^2-G^2)(\gamma_1^+ - \gamma_1^- - G)(\Gamma_1^+ - \Gamma_1^- - G)}
{\gamma_2\Gamma_2(3\gamma_1^+ - \gamma_1^- - G)(3\Gamma_1^+ - \Gamma_1^- - G)}
\right]^{-1}. \tag*{(32)} ⟨ b ⟩ ss = − γ 2 Γ 2 ( 3 γ 1 + − γ 1 − − G ) ( 3 Γ 1 + − Γ 1 − − G ) ( g + G ) ( γ 1 + − γ 1 − − G ) ( Γ 1 + − Γ 1 − − G ) Ω [ 1 + γ 2 Γ 2 ( 3 γ 1 + − γ 1 − − G ) ( 3 Γ 1 + − Γ 1 − − G ) ( g 2 − G 2 ) ( γ 1 + − γ 1 − − G ) ( Γ 1 + − Γ 1 − − G ) ] − 1 . (32) For the non-reciprocal coupling G = g G = g G = g M = b + b † M = b + b^{\dagger} M = b + b † Ω \Omega Ω
δ Ω ∣ n r ≃ 2 γ 2 Γ 2 ( 3 γ 1 + − γ 1 − − G ) ( 3 Γ 1 + − Γ 1 − − G ) 4 g ( γ 1 + − γ 1 − − G ) Γ 1 + − Γ 1 − − G . (33) \left.\delta\Omega\right|_{\mathrm{nr}} \simeq\frac{2\gamma_2\sqrt{\Gamma_2}(3\gamma_1^+ - \gamma_1^- - G)(3\Gamma_1^+ - \Gamma_1^- - G)}
{4g(\gamma_1^+ - \gamma_1^- - G)\sqrt{\Gamma_1^+ - \Gamma_1^- - G}}. \tag*{(33)} δ Ω ∣ nr ≃ 4 g ( γ 1 + − γ 1 − − G ) Γ 1 + − Γ 1 − − G 2 γ 2 Γ 2 ( 3 γ 1 + − γ 1 − − G ) ( 3 Γ 1 + − Γ 1 − − G ) . (33) This result indicates that when two-particle loss rate Γ 2 \Gamma_2 Γ 2
For the reciprocal coupling G = 0 G = 0 G = 0 g > 0 g > 0 g > 0 Γ 2 → 0 \Gamma_2 \to0 Γ 2 → 0 Ω \Omega Ω
δ Ω r ∝ 1 Γ 2 . (34) \delta\Omega_r \propto\frac{1}{\sqrt{\Gamma_2}}. \tag*{(34)} δ Ω r ∝ Γ 2 1 . (34) When the two-particle loss rate Γ 2 \Gamma_2 Γ 2 Γ 2 \Gamma_2 Γ 2 η \eta η Ω \Omega Ω
η = δ Ω r δ Ω ∣ n r ∝ 1 Γ 2 . (35) \eta= \frac{\delta\Omega_r}{\left.\delta\Omega\right|_{\mathrm{nr}}} \propto\frac{1}{\Gamma_2}. \tag*{(35)} η = δ Ω ∣ nr δ Ω r ∝ Γ 2 1 . (35) This result indicates that in the presence of strong incoherent drive and weak two-particle loss rate, the use of non-reciprocal coupling can effectively facilitate the sensing and measurement of weakly driving signals.
Critical incoherent drive
In the previous section, we have already explored the case of strong driving. In this part, we investigate the regime where incoherent drive is near its critical value. The linear rates obey the critical condition
γ 1 + − γ 1 − = Γ 1 + − Γ 1 − = G . (36) \gamma_1^+ - \gamma_1^- = \Gamma_1^+ - \Gamma_1^- = G. \tag*{(36)} γ 1 + − γ 1 − = Γ 1 + − Γ 1 − = G . (36) The reciprocal and nonreciprocal cases are compared at fixed γ 1 + \gamma_1^+ γ 1 + Γ 1 + \Gamma_1^+ Γ 1 + γ 2 \gamma_2 γ 2 g g g Ω \Omega Ω Γ 2 \Gamma_2 Γ 2 G = g G=g G = g
0 < g ≤ min ( γ 1 + , Γ 1 + ) . (37) 0 < g \le\min(\gamma_1^+,\Gamma_1^+). \tag*{(37)} 0 < g ≤ min ( γ 1 + , Γ 1 + ) . (37) The complete non-reciprocal coupling: G=g
For the complete non-reciprocal coupling, the large-occupation, time-scale-separated noise-dominated theory gives (see the detail in Appendix D)
G = g , Ω → 0 , Γ 2 → 0 : N b ∣ Ω = 0 ∝ Γ 2 − 1 / 2 , ∣ ∂ Ω ⟨ b ⟩ ∣ ∝ Γ 2 − 1 / 2 . (38) G = g, \Omega\to0, \Gamma_2 \to0:
\qquad
N_b|_{\Omega=0} \propto\Gamma_2^{-1/2},\quad|\partial_\Omega\langle b\rangle| \propto\Gamma_2^{-1/2}. \tag*{(38)} G = g , Ω → 0 , Γ 2 → 0 : N b ∣ Ω = 0 ∝ Γ 2 − 1/2 , ∣ ∂ Ω ⟨ b ⟩ ∣ ∝ Γ 2 − 1/2 . (38) At fixed nonzero drive, with nonzero upstream amplitude and finite effective noise, the asymptotic response is instead
G = g , Ω ≠ 0 , Γ 2 → 0 : ∣ ⟨ ∂ Ω ⟨ b ⟩ ⟩ ∣ ∝ Γ 2 − 1 / 3 , N b ∝ Γ 2 − 2 / 3 . (39) G = g, \Omega\ne0, \Gamma_2 \to0:
\qquad
|\langle\partial_\Omega\langle b\rangle\rangle| \propto\Gamma_2^{-1/3},\quad N_b \propto\Gamma_2^{-2/3}. \tag*{(39)} G = g , Ω = 0 , Γ 2 → 0 : ∣ ⟨ ∂ Ω ⟨ b ⟩⟩ ∣ ∝ Γ 2 − 1/3 , N b ∝ Γ 2 − 2/3 . (39) From Figure 3 , it can be verified that the derived relationship equations obtained from the above two scenarios are reasonable within a very narrow range of the Γ 2 \Gamma_2 Γ 2 Ω \Omega Ω
Figure 3. Evolution graph of the steady-state expectation value of model b b b Γ 2 \Gamma_2 Γ 2 γ 1 + = 1 \gamma_1^+ = 1 γ 1 + = 1 Γ 1 + = 0.50 \Gamma_1^+ = 0.50 Γ 1 + = 0.50 γ 2 = 0.0001 \gamma_2 = 0.0001 γ 2 = 0.0001 g = G = 0.2 g = G = 0.2 g = G = 0.2
δ Ω c ∝ Γ 2 1 / 4 , Ω → 0 , (40) \delta\Omega_c \propto\Gamma_2^{1/4},\quad\Omega\to0, \tag*{(40)} δ Ω c ∝ Γ 2 1/4 , Ω → 0 , (40) δ Ω c ∝ Γ 2 0 , Ω > 0. (41) \delta\Omega_c \propto\Gamma_2^0,\quad\Omega> 0. \tag*{(41)} δ Ω c ∝ Γ 2 0 , Ω > 0. (41) This result indicates that only when Ω \Omega Ω Ω \Omega Ω Γ 2 \Gamma_2 Γ 2 Ω = 0 \Omega= 0 Ω = 0 Ω \Omega Ω
δ Ω c / δ Ω ∣ n r ∝ Γ 2 1 / 4 , Ω → 0 , (42) \delta\Omega_c/\delta\Omega|_{nr} \propto\Gamma_2^{1/4},\quad\Omega\to0, \tag*{(42)} δ Ω c / δ Ω ∣ n r ∝ Γ 2 1/4 , Ω → 0 , (42) δ Ω c / δ Ω ∣ n r ∝ Γ 2 1 / 2 , Ω > 0. (43) \delta\Omega_c/\delta\Omega|_{nr} \propto\Gamma_2^{1/2},\quad\Omega> 0. \tag*{(43)} δ Ω c / δ Ω ∣ n r ∝ Γ 2 1/2 , Ω > 0. (43) This indicates that with the non-reciprocal coupling, increasing the strength of incoherent driving can effectively enhance the measurement precision of weak driving intensity.
The reciprocal coupling: G=0
In the absence of dissipative coupling, the two systems interact only via coherent coupling. According to the error propagation formula, the measurement uncertainty of Ω \Omega Ω
δ Ω = Var ( X ) ∣ ∂ ⟨ b ⟩ s s ∂ Ω ∣ = ∣ g ∣ 2 n b , 0 + 1 2 , (44) \delta\Omega= \frac{\sqrt{\operatorname{Var}(X)}}{\left|\frac{\partial\langle b\rangle_{\mathrm{ss}}}{\partial\Omega}\right|} = \frac{|g|}{\sqrt{2}}\sqrt{n_{b,0}+\frac{1}{2}}, \tag*{(44)} δ Ω = ∂ Ω ∂ ⟨ b ⟩ ss Var ( X ) = 2 ∣ g ∣ n b , 0 + 2 1 , (44) where the zero-drive (Ω = 0 \Omega= 0 Ω = 0 n b , 0 = ⟨ b † b ⟩ 0 n_{b,0} = \langle b^\dagger b\rangle_0 n b , 0 = ⟨ b † b ⟩ 0 [32 ] . As can be seen from Figure 4 , the approximate analytical result obtained here is essentially close to the result obtained through precise numerical calculation. Consequently,
Figure 4. The graph showing the relationship between the zero-drive occupation n b , 0 n_{b,0} n b , 0 g g g γ 1 + = γ 1 − = 0.50 \gamma_1^{+}=\gamma_1^{-}=0.50 γ 1 + = γ 1 − = 0.50 Γ 1 + = Γ 1 − = 0.50 \Gamma_1^{+}=\Gamma_1^{-}=0.50 Γ 1 + = Γ 1 − = 0.50 γ 2 = 0.2 \gamma_2=0.2 γ 2 = 0.2 Ω = 0 \Omega=0 Ω = 0 Γ 2 = 0 \Gamma_2=0 Γ 2 = 0 G = 0 G=0 G = 0
n b , 0 G a u s s = ( 1 + 4 Γ 1 + γ 2 g 2 ) γ 1 + + Γ 1 + 4 γ 2 . (45) n_{b,0}^{\mathrm{Gauss}}=\left(1+\frac{4\Gamma_1^{+}\gamma_2}{g^2}\right)\sqrt{\frac{\gamma_1^{+}+\Gamma_1^{+}}{4\gamma_2}}. \tag*{(45)} n b , 0 Gauss = ( 1 + g 2 4 Γ 1 + γ 2 ) 4 γ 2 γ 1 + + Γ 1 + . (45) Substituting Eq. (45) into Eq. (44) yields
δ Ω G a u s s = 1 2 [ ( g 2 + 4 Γ 1 + γ 2 ) γ 1 + + Γ 1 + 4 γ 2 + g 2 2 ] 1 / 2 . (46) \delta\Omega_{\mathrm{Gauss}}=\frac{1}{\sqrt{2}}\left[(g^2+4\Gamma_1^{+}\gamma_2)\sqrt{\frac{\gamma_1^{+}+\Gamma_1^{+}}{4\gamma_2}}+\frac{g^2}{2}\right]^{1/2}. \tag*{(46)} δ Ω Gauss = 2 1 ( g 2 + 4 Γ 1 + γ 2 ) 4 γ 2 γ 1 + + Γ 1 + + 2 g 2 1/2 . (46) This result is clearly independent of Γ 2 \Gamma_2 Γ 2 Γ 2 \Gamma_2 Γ 2
Conclusions
In this article, we investigate a quantum sensing system composed of two vdP oscillators. In the case of weak incoherent dissipation and driving, the measurement precision of the weak driving strength is independent of the two-particle loss rate. The non-reciprocal coupling can improve the measurement precision by a factor of two compared with the reciprocal coupling. When the coherent coupling, incoherent driving and dissipation are all zero, a dark mode emerges at the zero-driving point, which enables the system to achieve a measurement precision reaching the Heisenberg limit. When incoherent drive is much larger than the two-particle loss rate, the non-reciprocal coupling enables the measurement uncertainty for the weak coherent drive to decrease as the two-particle loss rate is reduced. The non-reciprocal coupling has a greater effect on improving measurement precision than the reciprocal coupling, and this effect increases as the two-particle loss rate decreases.
The cascaded quantum vdP oscillator discussed in this article can be realized on many different experimental platforms, such as trapped ion platforms [32 , 33 , 34 ] , superconducting circuits [35 , 36 , 37 ] , cavity optomechanics [38 , 39 , 40 ] , etc.
Acknowledgements
This research was supported by the National Natural Science Foundation of China (Grant No. 12365001), and the Bagui Youth Top Talent Training Program.
Non-Hermitian
The non-Hermitian evolution under the no-jump approximation can be written as
i d d t ∣ ψ ⟩ = H e f f ∣ ψ ⟩ (47) i\frac{d}{dt}|\psi\rangle=H_{\mathrm{eff}}|\psi\rangle \tag*{(47)} i d t d ∣ ψ ⟩ = H eff ∣ ψ ⟩ (47) where the effective Hamiltonian is
H e f f = i Ω ( a † − a ) + i g ( a † b − a b † ) − i { γ 1 + a a † + Γ 1 + a † a + γ 2 a † 2 a 2 + Γ 1 + b b † + γ 1 − b † b + γ 2 b † 2 b 2 + G ( a † + b † ) ( a + b ) } (48) H_{\mathrm{eff}}=i\Omega(a^\dagger-a)+ig(a^\dagger b-ab^\dagger)-i\left\{\gamma_1^{+}aa^\dagger+\Gamma_1^{+}a^\dagger a+\gamma_2a^{\dagger2}a^2+\Gamma_1^{+}bb^\dagger+\gamma_1^{-}b^\dagger b+\gamma_2b^{\dagger2}b^2+G(a^\dagger+b^\dagger)(a+b)\right\} \tag*{(48)} H eff = i Ω ( a † − a ) + i g ( a † b − a b † ) − i { γ 1 + a a † + Γ 1 + a † a + γ 2 a † 2 a 2 + Γ 1 + b b † + γ 1 − b † b + γ 2 b † 2 b 2 + G ( a † + b † ) ( a + b ) } (48) Using the bosonic commutation relations a a † = a † a + 1 aa^\dagger=a^\dagger a+1 a a † = a † a + 1 b b † = b † b + 1 bb^\dagger=b^\dagger b+1 b b † = b † b + 1
H e f f ′ = i Ω ( a † − a ) + i g ( a † b − a b † ) − i { γ 1 + + γ 1 − + μ a † a + ν b † b + γ 2 a † 2 a 2 + γ 2 b † 2 b 2 + G ( a † b + b † a ) } (49) H_{\mathrm{eff}}'=i\Omega(a^\dagger-a)+ig(a^\dagger b-ab^\dagger)-i\left\{\gamma_1^{+}+\gamma_1^{-}+\mu a^\dagger a+\nu b^\dagger b+\gamma_2a^{\dagger2}a^2+\gamma_2b^{\dagger2}b^2+G(a^\dagger b+b^\dagger a)\right\} \tag*{(49)} H eff ′ = i Ω ( a † − a ) + i g ( a † b − a b † ) − i { γ 1 + + γ 1 − + μ a † a + ν b † b + γ 2 a † 2 a 2 + γ 2 b † 2 b 2 + G ( a † b + b † a ) } (49) where μ = γ 1 + + γ 1 − + G \mu=\gamma_1^{+}+\gamma_1^{-}+G μ = γ 1 + + γ 1 − + G ν = Γ 1 + + Γ 1 − + G \nu=\Gamma_1^{+}+\Gamma_1^{-}+G ν = Γ 1 + + Γ 1 − + G γ 1 + + γ 1 − \gamma_1^{+}+\gamma_1^{-} γ 1 + + γ 1 −
We look for right eigenstates
H e f f ′ ∣ ψ ⟩ = λ ∣ ψ ⟩ . (50) H_{\mathrm{eff}}'|\psi\rangle=\lambda|\psi\rangle. \tag*{(50)} H eff ′ ∣ ψ ⟩ = λ ∣ ψ ⟩ . (50) In the weak-driving limit, λ = O ( Ω 2 ) \lambda=O(\Omega^2) λ = O ( Ω 2 ) λ c n m \lambda c_{nm} λ c nm
∣ ψ ⟩ = c 00 ∣ 00 ⟩ + c 01 ∣ 01 ⟩ + c 10 ∣ 10 ⟩ + c 11 ∣ 11 ⟩ + c 20 ∣ 20 ⟩ + c 02 ∣ 02 ⟩ + c 12 ∣ 12 ⟩ + c 21 ∣ 21 ⟩ + c 22 ∣ 22 ⟩ . (51) |\psi\rangle=c_{00}|00\rangle+c_{01}|01\rangle+c_{10}|10\rangle+c_{11}|11\rangle+c_{20}|20\rangle+c_{02}|02\rangle+c_{12}|12\rangle+c_{21}|21\rangle+c_{22}|22\rangle. \tag*{(51)} ∣ ψ ⟩ = c 00 ∣00 ⟩ + c 01 ∣01 ⟩ + c 10 ∣10 ⟩ + c 11 ∣11 ⟩ + c 20 ∣20 ⟩ + c 02 ∣02 ⟩ + c 12 ∣12 ⟩ + c 21 ∣21 ⟩ + c 22 ∣22 ⟩ . (51) Here, the annihilation operators satisfy
a ∣ n , m ⟩ = n ∣ n − 1 , m ⟩ , (52) a|n,m\rangle=\sqrt{n}|n-1,m\rangle, \tag*{(52)} a ∣ n , m ⟩ = n ∣ n − 1 , m ⟩ , (52) b ∣ n , m ⟩ = m ∣ n , m − 1 ⟩ . (53) b|n,m\rangle=\sqrt{m}|n,m-1\rangle. \tag*{(53)} b ∣ n , m ⟩ = m ∣ n , m − 1 ⟩ . (53) Substituting Eq. (51) into Eq. (50) , we obtain the coefficients
c 00 = 1 + O ( Ω 2 ) , c 10 = β Ω Δ 1 + O ( Ω ) , c_{00}=1+O(\Omega^2),\quad c_{10}=\frac{\beta\Omega}{\Delta_1}+O(\Omega), c 00 = 1 + O ( Ω 2 ) , c 10 = Δ 1 β Ω + O ( Ω ) , c 01 = − ( g + G ) Ω Δ 1 + O ( Ω ) , c_{01}=-\frac{(g+G)\Omega}{\Delta_1}+O(\Omega), c 01 = − Δ 1 ( g + G ) Ω + O ( Ω ) , c 20 , c 11 , c 02 = O ( Ω 2 ) , c 21 , c 12 = O ( Ω 3 ) , c 22 = O ( Ω 4 ) , (54) c_{20},c_{11},c_{02}=O(\Omega^2),\quad c_{21},c_{12}=O(\Omega^3),\quad c_{22}=O(\Omega^4), \tag*{(54)} c 20 , c 11 , c 02 = O ( Ω 2 ) , c 21 , c 12 = O ( Ω 3 ) , c 22 = O ( Ω 4 ) , (54) where Δ 1 = μ ν + g 2 − G 2 \Delta_1=\mu\nu+g^2-G^2 Δ 1 = μν + g 2 − G 2
⟨ a ⟩ = ⟨ ψ ∣ a ∣ ψ ⟩ , ⟨ b ⟩ = ⟨ ψ ∣ b ∣ ψ ⟩ . (55) \langle a\rangle=\langle\psi|a|\psi\rangle,\qquad\langle b\rangle=\langle\psi|b|\psi\rangle. \tag*{(55)} ⟨ a ⟩ = ⟨ ψ ∣ a ∣ ψ ⟩ , ⟨ b ⟩ = ⟨ ψ ∣ b ∣ ψ ⟩ . (55) Explicitly,
⟨ a ⟩ = c 00 ∗ c 10 + 2 c 10 ∗ c 20 + c 01 ∗ c 11 + 2 c 11 ∗ c 21 + 2 c 02 ∗ c 12 + 2 c 12 ∗ c 22 . (56) \langle a\rangle=c_{00}^{*}c_{10}+\sqrt{2}c_{10}^{*}c_{20}+c_{01}^{*}c_{11}+\sqrt{2}c_{11}^{*}c_{21}
+\sqrt{2}c_{02}^{*}c_{12}+\sqrt{2}c_{12}^{*}c_{22}. \tag*{(56)} ⟨ a ⟩ = c 00 ∗ c 10 + 2 c 10 ∗ c 20 + c 01 ∗ c 11 + 2 c 11 ∗ c 21 + 2 c 02 ∗ c 12 + 2 c 12 ∗ c 22 . (56) ⟨ b ⟩ = c 00 ∗ c 01 + c 10 ∗ c 11 + c 20 ∗ c 21 + 2 c 01 ∗ c 02 + 2 c 11 ∗ c 12 + 2 c 21 ∗ c 22 . (57) \langle b\rangle=c_{00}^{*}c_{01}+c_{10}^{*}c_{11}+c_{20}^{*}c_{21}
+\sqrt{2}c_{01}^{*}c_{02}+\sqrt{2}c_{11}^{*}c_{12}+\sqrt{2}c_{21}^{*}c_{22}. \tag*{(57)} ⟨ b ⟩ = c 00 ∗ c 01 + c 10 ∗ c 11 + c 20 ∗ c 21 + 2 c 01 ∗ c 02 + 2 c 11 ∗ c 12 + 2 c 21 ∗ c 22 . (57) At lowest order O ( Ω ) O(\Omega) O ( Ω )
⟨ a ⟩ = β Ω Δ 1 + O ( Ω 3 ) (58) \langle a\rangle=\frac{\beta\Omega}{\Delta_1}+O(\Omega^3) \tag*{(58)} ⟨ a ⟩ = Δ 1 β Ω + O ( Ω 3 ) (58) and
⟨ b ⟩ = − ( g + G ) Ω Δ 1 + O ( Ω 3 ) . (59) \langle b\rangle=-\frac{(g+G)\Omega}{\Delta_1}+O(\Omega^3). \tag*{(59)} ⟨ b ⟩ = − Δ 1 ( g + G ) Ω + O ( Ω 3 ) . (59) ⟨ b † b ⟩ = − ( g + G ) 2 Ω 2 g 4 + O ( Ω 3 ) . (60) \langle b^\dagger b\rangle=-\frac{(g+G)^2\Omega^2}{g^4}+O(\Omega^3). \tag*{(60)} ⟨ b † b ⟩ = − g 4 ( g + G ) 2 Ω 2 + O ( Ω 3 ) . (60) Appendix B: characteristic time of dark mode
The quantum-Langevin equations of the operators a a a b b b
d d t a = Ω − G ( a + b ) − 2 γ 2 a † a 2 + 2 G z i n ( t ) + 2 γ 2 a † a i n ( t ) , (61) \frac{d}{dt}a=\Omega-G(a+b)-2\gamma_2a^\dagger a^2+\sqrt{2G}z_{\mathrm{in}}(t)+\sqrt{2\gamma_2}a^\dagger a_{\mathrm{in}}(t), \tag*{(61)} d t d a = Ω − G ( a + b ) − 2 γ 2 a † a 2 + 2 G z in ( t ) + 2 γ 2 a † a in ( t ) , (61) d d t b = − G ( a + b ) − 2 γ 2 b † b 2 + 2 G z i n ( t ) + 2 γ 2 b † b i n ( t ) . (62) \frac{d}{dt}b=-G(a+b)-2\gamma_2b^\dagger b^2+\sqrt{2G}z_{\mathrm{in}}(t)+\sqrt{2\gamma_2}b^\dagger b_{\mathrm{in}}(t). \tag*{(62)} d t d b = − G ( a + b ) − 2 γ 2 b † b 2 + 2 G z in ( t ) + 2 γ 2 b † b in ( t ) . (62) We express a a a b b b a = ⟨ a ⟩ + δ a a=\langle a\rangle+\delta a a = ⟨ a ⟩ + δ a b = ⟨ b ⟩ + δ b b=\langle b\rangle+\delta b b = ⟨ b ⟩ + δ b δ a \delta a δ a δ b \delta b δ b
d d t δ a = − G ( δ a + δ b ) − 2 γ 2 ⟨ a ⟩ 2 δ a † − 4 γ 2 ∣ ⟨ a ⟩ ∣ 2 δ a + A i n ( t ) , (63) \frac{d}{dt}\delta a=-G(\delta a+\delta b)-2\gamma_2\langle a\rangle^2\delta a^\dagger-4\gamma_2|\langle a\rangle|^2\delta a+A_{\mathrm{in}}(t), \tag*{(63)} d t d δ a = − G ( δ a + δ b ) − 2 γ 2 ⟨ a ⟩ 2 δ a † − 4 γ 2 ∣ ⟨ a ⟩ ∣ 2 δ a + A in ( t ) , (63) d d t δ b = − G ( δ a + δ b ) − 2 γ 2 ⟨ b ⟩ 2 δ b † − 4 γ 2 ∣ ⟨ b ⟩ ∣ 2 δ b + B i n ( t ) , (64) \frac{d}{dt}\delta b=-G(\delta a+\delta b)-2\gamma_2\langle b\rangle^2\delta b^\dagger-4\gamma_2|\langle b\rangle|^2\delta b+B_{\mathrm{in}}(t), \tag*{(64)} d t d δ b = − G ( δ a + δ b ) − 2 γ 2 ⟨ b ⟩ 2 δ b † − 4 γ 2 ∣ ⟨ b ⟩ ∣ 2 δ b + B in ( t ) , (64) where the quantum noise operators are given by
A i n ( t ) = 2 G z i n ( t ) + 2 γ 2 ⟨ a † ⟩ a i n ( t ) , (65) A_{\mathrm{in}}(t)=\sqrt{2G}z_{\mathrm{in}}(t)+\sqrt{2\gamma_2}\langle a^\dagger\rangle a_{\mathrm{in}}(t), \tag*{(65)} A in ( t ) = 2 G z in ( t ) + 2 γ 2 ⟨ a † ⟩ a in ( t ) , (65) B i n ( t ) = 2 G z i n ( t ) + 2 γ 2 ⟨ b † ⟩ b i n ( t ) . (66) B_{\mathrm{in}}(t)=\sqrt{2G}z_{\mathrm{in}}(t)+\sqrt{2\gamma_2}\langle b^\dagger\rangle b_{\mathrm{in}}(t). \tag*{(66)} B in ( t ) = 2 G z in ( t ) + 2 γ 2 ⟨ b † ⟩ b in ( t ) . (66) Let us set
r = ( δ a δ a † δ b δ b † ) , F i n ( t ) = ( A i n ( t ) A i n † ( t ) B i n ( t ) B i n † ( t ) ) , \mathbf{r}=
\begin{pmatrix}
\delta a\\
\delta a^\dagger\\
\delta b\\
\delta b^\dagger
\end{pmatrix},
\qquad
\mathbf{F}_{\mathrm{in}}(t)=
\begin{pmatrix}
A_{\mathrm{in}}(t)\\
A_{\mathrm{in}}^\dagger(t)\\
B_{\mathrm{in}}(t)\\
B_{\mathrm{in}}^\dagger(t)
\end{pmatrix}, r = δ a δ a † δ b δ b † , F in ( t ) = A in ( t ) A in † ( t ) B in ( t ) B in † ( t ) , quantum-Langevin equation is rewritten as
d d t r = M r + F i n ( t ) , (67) \frac{d}{dt}\mathbf{r}=\mathbf{M}\mathbf{r}+\mathbf{F}_{\mathrm{in}}(t), \tag*{(67)} d t d r = Mr + F in ( t ) , (67) in which, the evolution matrix:
M = ( − 4 γ 2 ∣ ⟨ a ⟩ ∣ 2 − G − 2 γ 2 ⟨ a ⟩ 2 − G 0 − 2 γ 2 ( ⟨ a ⟩ ∗ ) 2 − 4 γ 2 ∣ ⟨ a ⟩ ∣ 2 − G 0 − G − G 0 − 4 γ 2 ∣ ⟨ b ⟩ ∣ 2 − G − 2 γ 2 ⟨ b ⟩ 2 0 − G − 2 γ 2 ( ⟨ b ⟩ ∗ ) 2 − 4 γ 2 ∣ ⟨ b ⟩ ∣ 2 − G ) . (68) \mathbf{M}=
\begin{pmatrix}
-4\gamma_2|\langle a\rangle|^2-G & -2\gamma_2\langle a\rangle^2 & -G & 0\\
-2\gamma_2(\langle a\rangle^*)^2 & -4\gamma_2|\langle a\rangle|^2-G & 0 & -G\\
-G & 0 & -4\gamma_2|\langle b\rangle|^2-G & -2\gamma_2\langle b\rangle^2\\
0 & -G & -2\gamma_2(\langle b\rangle^*)^2 & -4\gamma_2|\langle b\rangle|^2-G
\end{pmatrix}. \tag*{(68)} M = − 4 γ 2 ∣ ⟨ a ⟩ ∣ 2 − G − 2 γ 2 (⟨ a ⟩ ∗ ) 2 − G 0 − 2 γ 2 ⟨ a ⟩ 2 − 4 γ 2 ∣ ⟨ a ⟩ ∣ 2 − G 0 − G − G 0 − 4 γ 2 ∣ ⟨ b ⟩ ∣ 2 − G − 2 γ 2 (⟨ b ⟩ ∗ ) 2 0 − G − 2 γ 2 ⟨ b ⟩ 2 − 4 γ 2 ∣ ⟨ b ⟩ ∣ 2 − G . (68) As Ω \Omega Ω M \mathbf{M} M
λ 1 = − ( γ 2 Ω 2 2 ) 1 / 3 + O ( Ω ) , (69) \lambda_1=-\left(\frac{\gamma_2\Omega^2}{2}\right)^{1/3}+O(\Omega), \tag*{(69)} λ 1 = − ( 2 γ 2 Ω 2 ) 1/3 + O ( Ω ) , (69) λ 2 = − 3 ( γ 2 Ω 2 2 ) 1 / 3 + O ( Ω ) , (70) \lambda_2=-3\left(\frac{\gamma_2\Omega^2}{2}\right)^{1/3}+O(\Omega), \tag*{(70)} λ 2 = − 3 ( 2 γ 2 Ω 2 ) 1/3 + O ( Ω ) , (70) λ 3 = − 2 G − ( γ 2 Ω 2 2 ) 1 / 3 + O ( Ω ) , (71) \lambda_3=-2G-\left(\frac{\gamma_2\Omega^2}{2}\right)^{1/3}+O(\Omega), \tag*{(71)} λ 3 = − 2 G − ( 2 γ 2 Ω 2 ) 1/3 + O ( Ω ) , (71) λ 4 = − 2 G − 3 ( γ 2 Ω 2 2 ) 1 / 3 + O ( Ω ) . (72) \lambda_4=-2G-3\left(\frac{\gamma_2\Omega^2}{2}\right)^{1/3}+O(\Omega). \tag*{(72)} λ 4 = − 2 G − 3 ( 2 γ 2 Ω 2 ) 1/3 + O ( Ω ) . (72) The characteristic time for reaching the steady state is determined by the eigenvalue with the smallest absolute value:
T c = 1 ∣ λ 1 ∣ = ( γ 2 Ω 2 2 ) − 1 / 3 . (73) T_c = \frac{1}{\lvert\lambda_1\rvert} = \left(\frac{\gamma_2\Omega^2}{2}\right)^{-1/3}. \tag*{(73)} T c = ∣ λ 1 ∣ 1 = ( 2 γ 2 Ω 2 ) − 1/3 . (73) Weak-coupling Kramers–Moyal analysis
In this Appendix, we formulate the weak-drive response of two quantum vdP oscillators coupled via coherent exchange and collective dissipation. We distinguish exact moment identities from the single-oscillator asymptotics and from additional closure assumptions.
We consider two bosonic modes with annihilation operators a a a b b b
ρ ˙ = [ Ω ( a † − a ) + g ( a † b − a b † ) , ρ ] + γ 1 + D [ a † ] ρ + γ 1 − D [ a ] ρ + γ 2 D [ a 2 ] ρ + Γ 1 + D [ b † ] ρ + Γ 1 − D [ b ] ρ + Γ 2 D [ b 2 ] ρ + G D [ a + b ] ρ . (74) \dot{\rho} = \left[\Omega(a^\dagger-a)+g(a^\dagger b-ab^\dagger),\rho\right]+\gamma_1^+\mathcal{D}[a^\dagger]\rho+\gamma_1^-\mathcal{D}[a]\rho+\gamma_2\mathcal{D}[a^2]\rho+\Gamma_1^+\mathcal{D}[b^\dagger]\rho+\Gamma_1^-\mathcal{D}[b]\rho+\Gamma_2\mathcal{D}[b^2]\rho+G\mathcal{D}[a+b]\rho. \tag*{(74)} ρ ˙ = [ Ω ( a † − a ) + g ( a † b − a b † ) , ρ ] + γ 1 + D [ a † ] ρ + γ 1 − D [ a ] ρ + γ 2 D [ a 2 ] ρ + Γ 1 + D [ b † ] ρ + Γ 1 − D [ b ] ρ + Γ 2 D [ b 2 ] ρ + G D [ a + b ] ρ . (74) where
D [ c ] ρ = 2 c ρ c † − c † c ρ − ρ c † c . (75) \mathcal{D}[c]\rho=2c\rho c^\dagger-c^\dagger c\rho-\rho c^\dagger c. \tag*{(75)} D [ c ] ρ = 2 c ρ c † − c † c ρ − ρ c † c . (75) We assume weak coherent drive and weak interoscillator coupling. The Kramers–Moyal approximation is used in the large-occupation regime, away from the critical point where the net linear gain vanishes.
Collective dissipator
The collective dissipator can be decomposed as
D [ a + b ] ρ = D [ a ] ρ + D [ b ] ρ + D c r [ a , b ] ρ , (76) \mathcal{D}[a+b]\rho=\mathcal{D}[a]\rho+\mathcal{D}[b]\rho+\mathcal{D}_{\mathrm{cr}}[a,b]\rho, \tag*{(76)} D [ a + b ] ρ = D [ a ] ρ + D [ b ] ρ + D cr [ a , b ] ρ , (76) D c r [ a , b ] ρ = 2 a ρ b † + 2 b ρ a † − a † b ρ − b † a ρ − ρ a † b − ρ b † a . (77) \mathcal{D}_{\mathrm{cr}}[a,b]\rho=2a\rho b^\dagger+2b\rho a^\dagger-a^\dagger b\rho-b^\dagger a\rho-\rho a^\dagger b-\rho b^\dagger a. \tag*{(77)} D cr [ a , b ] ρ = 2 a ρ b † + 2 b ρ a † − a † b ρ − b † a ρ − ρ a † b − ρ b † a . (77) We define the effective local one-particle loss rates
μ a = γ 1 − + G , μ b = Γ 1 − + G , (78) \mu_a=\gamma_1^-+G,\qquad\mu_b=\Gamma_1^-+G, \tag*{(78)} μ a = γ 1 − + G , μ b = Γ 1 − + G , (78) and the corresponding net linear gain rates
δ a = γ 1 + − μ a = γ 1 + − γ 1 − − G , (79) \delta_a=\gamma_1^+-\mu_a=\gamma_1^+-\gamma_1^--G, \tag*{(79)} δ a = γ 1 + − μ a = γ 1 + − γ 1 − − G , (79) δ b = Γ 1 + − μ b = Γ 1 + − Γ 1 − − G . (80) \delta_b=\Gamma_1^+-\mu_b=\Gamma_1^+-\Gamma_1^--G. \tag*{(80)} δ b = Γ 1 + − μ b = Γ 1 + − Γ 1 − − G . (80) The local Liouvillians are
L a ρ = γ 1 + D [ a † ] ρ + μ a D [ a ] ρ + γ 2 D [ a 2 ] ρ , (81) \mathcal{L}_a\rho=\gamma_1^+\mathcal{D}[a^\dagger]\rho+\mu_a\mathcal{D}[a]\rho+\gamma_2\mathcal{D}[a^2]\rho, \tag*{(81)} L a ρ = γ 1 + D [ a † ] ρ + μ a D [ a ] ρ + γ 2 D [ a 2 ] ρ , (81) L b ρ = Γ 1 + D [ b † ] ρ + μ b D [ b ] ρ + Γ 2 D [ b 2 ] ρ . (82) \mathcal{L}_b\rho=\Gamma_1^+\mathcal{D}[b^\dagger]\rho+\mu_b\mathcal{D}[b]\rho+\Gamma_2\mathcal{D}[b^2]\rho. \tag*{(82)} L b ρ = Γ 1 + D [ b † ] ρ + μ b D [ b ] ρ + Γ 2 D [ b 2 ] ρ . (82) The term that couples the two oscillators is
L c ρ = [ g ( a † b − a b † ) , ρ ] + G D c r [ a , b ] ρ . (83) \mathcal{L}_c\rho=[g(a^\dagger b-ab^\dagger),\rho]+G\mathcal{D}_{\mathrm{cr}}[a,b]\rho. \tag*{(83)} L c ρ = [ g ( a † b − a b † ) , ρ ] + G D cr [ a , b ] ρ . (83) Population Kramers–Moyal equation
For oscillator a a a
p n ( a ) = ⟨ n ∣ ρ a ∣ n ⟩ (84) p_n^{(a)}=\langle n\lvert\rho_a\lvert n\rangle \tag*{(84)} p n ( a ) = ⟨ n ∣ ρ a ∣ n ⟩ (84) denote the diagonal Fock-state population. The relevant transitions are
n ⟶ n + 1 , W + ( a ) ( n ) = 2 γ 1 + ( n + 1 ) , (85) n\longrightarrow n+1,\qquad W_+^{(a)}(n)=2\gamma_1^+(n+1), \tag*{(85)} n ⟶ n + 1 , W + ( a ) ( n ) = 2 γ 1 + ( n + 1 ) , (85) n ⟶ n − 1 , W − ( a ) ( n ) = 2 μ a n , (86) n\longrightarrow n-1,\qquad W_-^{(a)}(n)=2\mu_a n, \tag*{(86)} n ⟶ n − 1 , W − ( a ) ( n ) = 2 μ a n , (86) n ⟶ n − 2 , W 2 ( a ) ( n ) = 2 γ 2 n ( n − 1 ) . (87) n\longrightarrow n-2,\qquad W_2^{(a)}(n)=2\gamma_2 n(n-1). \tag*{(87)} n ⟶ n − 2 , W 2 ( a ) ( n ) = 2 γ 2 n ( n − 1 ) . (87) The exact population equation is
p ˙ n ( a ) = 2 γ 1 + [ n p n − 1 ( a ) − ( n + 1 ) p n ( a ) ] + 2 μ a [ ( n + 1 ) p n + 1 ( a ) − n p n ( a ) ] + 2 γ 2 [ ( n + 1 ) ( n + 2 ) p n + 2 ( a ) − n ( n − 1 ) p n ( a ) ] . (88) \dot{p}_n^{(a)}=2\gamma_1^+\left[np_{n-1}^{(a)}-(n+1)p_n^{(a)}\right]+2\mu_a\left[(n+1)p_{n+1}^{(a)}-np_n^{(a)}\right]+2\gamma_2\left[(n+1)(n+2)p_{n+2}^{(a)}-n(n-1)p_n^{(a)}\right]. \tag*{(88)} p ˙ n ( a ) = 2 γ 1 + [ n p n − 1 ( a ) − ( n + 1 ) p n ( a ) ] + 2 μ a [ ( n + 1 ) p n + 1 ( a ) − n p n ( a ) ] + 2 γ 2 [ ( n + 1 ) ( n + 2 ) p n + 2 ( a ) − n ( n − 1 ) p n ( a ) ] . (88) For a jump process with jump size Δ n \Delta n Δ n
A k ( n ) = ∑ Δ n ( Δ n ) k W ( n , Δ n ) . (89) A_k(n)=\sum_{\Delta n}(\Delta n)^kW(n,\Delta n). \tag*{(89)} A k ( n ) = Δ n ∑ ( Δ n ) k W ( n , Δ n ) . (89) The first two Kramers–Moyal coefficients are
A 1 ( a ) ( n ) = 2 γ 1 + ( n + 1 ) − 2 μ a n − 4 γ 2 n ( n − 1 ) , (90) A_1^{(a)}(n)=2\gamma_1^+(n+1)-2\mu_a n-4\gamma_2n(n-1), \tag*{(90)} A 1 ( a ) ( n ) = 2 γ 1 + ( n + 1 ) − 2 μ a n − 4 γ 2 n ( n − 1 ) , (90) A 2 ( a ) ( n ) = 2 γ 1 + ( n + 1 ) + 2 μ a n + 8 γ 2 n ( n − 1 ) . (91) A_2^{(a)}(n)=2\gamma_1^+(n+1)+2\mu_a n+8\gamma_2n(n-1). \tag*{(91)} A 2 ( a ) ( n ) = 2 γ 1 + ( n + 1 ) + 2 μ a n + 8 γ 2 n ( n − 1 ) . (91) The corresponding Fokker–Planck equation is
∂ p a ( n , t ) ∂ t = − ∂ ∂ n [ A 1 ( a ) ( n ) p a ( n , t ) ] + 1 2 ∂ 2 ∂ n 2 [ A 2 ( a ) ( n ) p a ( n , t ) ] . (92) \frac{\partial p_a(n,t)}{\partial t}=-\frac{\partial}{\partial n}\left[A_1^{(a)}(n)p_a(n,t)\right]+\frac{1}{2}\frac{\partial^2}{\partial n^2}\left[A_2^{(a)}(n)p_a(n,t)\right]. \tag*{(92)} ∂ t ∂ p a ( n , t ) = − ∂ n ∂ [ A 1 ( a ) ( n ) p a ( n , t ) ] + 2 1 ∂ n 2 ∂ 2 [ A 2 ( a ) ( n ) p a ( n , t ) ] . (92) For oscillator b b b
γ 1 + ⟶ Γ 1 + , μ a ⟶ μ b , γ 2 ⟶ Γ 2 . \gamma_1^+\longrightarrow\Gamma_1^+,\qquad\mu_a\longrightarrow\mu_b,\qquad\gamma_2\longrightarrow\Gamma_2. γ 1 + ⟶ Γ 1 + , μ a ⟶ μ b , γ 2 ⟶ Γ 2 . Off-critical limit-cycle occupation
We first assume
δ a > 0 , δ b > 0 , (93) \delta_a>0,\qquad\delta_b>0, \tag*{(93)} δ a > 0 , δ b > 0 , (93) and sufficiently large occupations.
For n ≫ 1 n\gg1 n ≫ 1 a a a
A 1 ( a ) ( n ) ≃ 2 δ a n − 4 γ 2 n 2 . (94) A_1^{(a)}(n)\simeq2\delta_a n-4\gamma_2n^2. \tag*{(94)} A 1 ( a ) ( n ) ≃ 2 δ a n − 4 γ 2 n 2 . (94) The nonzero stable drift zero is given by
n a , ∗ = δ a 2 γ 2 . (95) n_{a,*}=\frac{\delta_a}{2\gamma_2}. \tag*{(95)} n a , ∗ = 2 γ 2 δ a . (95) Linearizing around n a , ∗ n_{a,*} n a , ∗ n = n a , ∗ + x n=n_{a,*}+x n = n a , ∗ + x
A 1 ( a ) ( n ) ≃ − 2 δ a x . (96) A_1^{(a)}(n)\simeq-2\delta_a x. \tag*{(96)} A 1 ( a ) ( n ) ≃ − 2 δ a x . (96) At the drift zero,
A 2 ( a ) ( n a , ∗ ) ≃ 2 ( γ 1 + + μ a ) n a , ∗ + 8 γ 2 n a , ∗ 2 = δ a γ 2 ( 3 γ 1 + − μ a ) . (97) A_2^{(a)}(n_{a,*})\simeq2(\gamma_1^+ +\mu_a)n_{a,*}+8\gamma_2 n_{a,*}^2
=\frac{\delta_a}{\gamma_2}(3\gamma_1^+-\mu_a). \tag*{(97)} A 2 ( a ) ( n a , ∗ ) ≃ 2 ( γ 1 + + μ a ) n a , ∗ + 8 γ 2 n a , ∗ 2 = γ 2 δ a ( 3 γ 1 + − μ a ) . (97) The resulting Gaussian approximation has variance
Var ( n a ) ≃ A 2 ( a ) ( n a , ∗ ) 4 δ a = 3 γ 1 + − μ a 4 γ 2 . (98) \operatorname{Var}(n_a)\simeq\frac{A_2^{(a)}(n_{a,*})}{4\delta_a}
=\frac{3\gamma_1^+-\mu_a}{4\gamma_2}. \tag*{(98)} Var ( n a ) ≃ 4 δ a A 2 ( a ) ( n a , ∗ ) = 4 γ 2 3 γ 1 + − μ a . (98) The Gaussian approximation requires
n a , ∗ 2 ≫ Var ( n a ) , (99) n_{a,*}^2\gg\operatorname{Var}(n_a), \tag*{(99)} n a , ∗ 2 ≫ Var ( n a ) , (99) or, parametrically,
δ a 2 ≫ γ 2 ( 3 γ 1 + − μ a ) . (100) \delta_a^2\gg\gamma_2(3\gamma_1^+-\mu_a). \tag*{(100)} δ a 2 ≫ γ 2 ( 3 γ 1 + − μ a ) . (100) Therefore,
⟨ a † a ⟩ 0 ≃ δ a 2 γ 2 . (101) \langle a^\dagger a\rangle_0\simeq\frac{\delta_a}{2\gamma_2}. \tag*{(101)} ⟨ a † a ⟩ 0 ≃ 2 γ 2 δ a . (101) For oscillator b b b
⟨ b † b ⟩ 0 ≃ δ b 2 Γ 2 , (102) \langle b^\dagger b\rangle_0\simeq\frac{\delta_b}{2\Gamma_2}, \tag*{(102)} ⟨ b † b ⟩ 0 ≃ 2 Γ 2 δ b , (102) and
Var ( n b ) ≃ 3 Γ 1 + − μ b 4 Γ 2 . (103) \operatorname{Var}(n_b)\simeq\frac{3\Gamma_1^+-\mu_b}{4\Gamma_2}. \tag*{(103)} Var ( n b ) ≃ 4 Γ 2 3 Γ 1 + − μ b . (103) Explicit local susceptibility
We next calculate the local zero-frequency susceptibility in the same large-occupation approximation. We write
ρ a = ρ a ( 0 ) + F a ρ a ( 1 ) + O ( F a 2 ) , L a ρ a ( 0 ) = 0 , (104) \rho_a=\rho_a^{(0)}+F_a\rho_a^{(1)}+O(F_a^2),\qquad\mathcal{L}_a\rho_a^{(0)}=0, \tag*{(104)} ρ a = ρ a ( 0 ) + F a ρ a ( 1 ) + O ( F a 2 ) , L a ρ a ( 0 ) = 0 , (104) where F a F_a F a
L a ρ a ( 1 ) = − [ a † − a , ρ a ( 0 ) ] . (105) \mathcal{L}_a\rho_a^{(1)}=-[a^\dagger-a,\rho_a^{(0)}]. \tag*{(105)} L a ρ a ( 1 ) = − [ a † − a , ρ a ( 0 ) ] . (105) Define the nearest off-diagonal matrix elements
x n ( a ) = ⟨ n ∣ ρ a ( 1 ) ∣ n − 1 ⟩ , q n ( a ) = n x n ( a ) . (106) x_n^{(a)}=\langle n|\rho_a^{(1)}|n-1\rangle,\qquad q_n^{(a)}=\sqrt{n}x_n^{(a)}. \tag*{(106)} x n ( a ) = ⟨ n ∣ ρ a ( 1 ) ∣ n − 1 ⟩ , q n ( a ) = n x n ( a ) . (106) The local susceptibility is
⟨ a ⟩ = χ a F a , χ a = ∑ n = 1 ∞ q n ( a ) . (107) \langle a\rangle=\chi_aF_a,\qquad\chi_a=\sum_{n=1}^{\infty}q_n^{(a)}. \tag*{(107)} ⟨ a ⟩ = χ a F a , χ a = n = 1 ∑ ∞ q n ( a ) . (107) Projecting the first-order equation onto the ⟨ n ∣ ⋯ ∣ n − 1 ⟩ \langle n|\cdots|n-1\rangle ⟨ n ∣ ⋯ ∣ n − 1 ⟩
γ 1 + [ 2 n ( n − 1 ) x n − 1 ( a ) − ( 2 n + 1 ) x n ( a ) ] + μ a [ 2 n ( n + 1 ) x n + 1 ( a ) − ( 2 n − 1 ) x n ( a ) ] + γ 2 [ 2 ( n + 1 ) n ( n + 2 ) x n + 2 ( a ) − 2 ( n − 1 ) 2 x n ( a ) ] + n ( p n − 1 ( a ) − p n ( a ) ) = 0. (108) \gamma_1^+\left[2\sqrt{n(n-1)}x_{n-1}^{(a)}-(2n+1)x_n^{(a)}\right]
+\mu_a\left[2\sqrt{n(n+1)}x_{n+1}^{(a)}-(2n-1)x_n^{(a)}\right]
+\gamma_2\left[2(n+1)\sqrt{n(n+2)}x_{n+2}^{(a)}-2(n-1)^2x_n^{(a)}\right]
+\sqrt{n}\left(p_{n-1}^{(a)}-p_n^{(a)}\right)=0. \tag*{(108)} γ 1 + [ 2 n ( n − 1 ) x n − 1 ( a ) − ( 2 n + 1 ) x n ( a ) ] + μ a [ 2 n ( n + 1 ) x n + 1 ( a ) − ( 2 n − 1 ) x n ( a ) ] + γ 2 [ 2 ( n + 1 ) n ( n + 2 ) x n + 2 ( a ) − 2 ( n − 1 ) 2 x n ( a ) ] + n ( p n − 1 ( a ) − p n ( a ) ) = 0. (108) The source term is proportional to p n − 1 ( a ) − p n ( a ) p_{n-1}^{(a)}-p_n^{(a)} p n − 1 ( a ) − p n ( a )
⟨ n ∣ [ a † − a , ρ a ( 0 ) ] ∣ n − 1 ⟩ = n ( p n − 1 ( a ) − p n ( a ) ) . \langle n|[a^\dagger-a,\rho_a^{(0)}]|n-1\rangle
=\sqrt{n}\left(p_{n-1}^{(a)}-p_n^{(a)}\right). ⟨ n ∣ [ a † − a , ρ a ( 0 ) ] ∣ n − 1 ⟩ = n ( p n − 1 ( a ) − p n ( a ) ) . In the large-occupation limit, the recurrence relation reduces to a differential equation for the first angular-harmonic coherence envelope q a ( n ) q_a(n) q a ( n )
K a q a ( n ) − Λ a q a ( n ) = − p a ′ ( n ) + 1 2 p a ′ ′ ( n ) , (109) \mathcal{K}_a q_a(n)-\Lambda_a q_a(n)=-p_a'(n)+\frac{1}{2}p_a''(n), \tag*{(109)} K a q a ( n ) − Λ a q a ( n ) = − p a ′ ( n ) + 2 1 p a ′′ ( n ) , (109) where
K a q = − ∂ n [ A 1 ( a ) q ] + 1 2 ∂ n 2 [ A 2 ( a ) q ] , (110) \mathcal{K}_a q=-\partial_n\left[A_1^{(a)}q\right]+\frac{1}{2}\partial_n^2\left[A_2^{(a)}q\right], \tag*{(110)} K a q = − ∂ n [ A 1 ( a ) q ] + 2 1 ∂ n 2 [ A 2 ( a ) q ] , (110) and
Λ a = 3 γ 1 + − μ a . (111) \Lambda_a=3\gamma_1^+-\mu_a. \tag*{(111)} Λ a = 3 γ 1 + − μ a . (111) The population distribution p a ( n ) p_a(n) p a ( n )
K a p a ( n ) = 0. (112) \mathcal{K}_a p_a(n)=0. \tag*{(112)} K a p a ( n ) = 0. (112) A pointwise ansatz q a ( n ) = χ a p a ( n ) q_a(n)=\chi_a p_a(n) q a ( n ) = χ a p a ( n ) (109) , because the driving source contains higher Hermite components. Instead, we retain the leading population mode and the first nonvanishing odd Hermite component generated by the source. Denoting the corresponding projection by P 1 P_1 P 1
P 1 [ K a q a − Λ a q a ] = P 1 [ − p a ′ + 1 2 p a ′ ′ ] . (113) P_1\left[\mathcal{K}_a q_a-\Lambda_a q_a\right]
=P_1\left[-p_a'+\frac{1}{2}p_a''\right]. \tag*{(113)} P 1 [ K a q a − Λ a q a ] = P 1 [ − p a ′ + 2 1 p a ′′ ] . (113) Using the Gaussian population distribution with
σ a 2 = 3 γ 1 + − μ a 4 γ 2 , \sigma_a^2=\frac{3\gamma_1^+-\mu_a}{4\gamma_2}, σ a 2 = 4 γ 2 3 γ 1 + − μ a , the projected source and response give the closure relation
γ 2 Λ a δ a χ a = 1. (114) \frac{\gamma_2\Lambda_a}{\delta_a}\chi_a=1. \tag*{(114)} δ a γ 2 Λ a χ a = 1. (114) Therefore, the explicit leading-order susceptibility of oscillator a a a
χ a = δ a γ 2 ( 3 γ 1 + − μ a ) . (115) \chi_a=\frac{\delta_a}{\gamma_2(3\gamma_1^+-\mu_a)}. \tag*{(115)} χ a = γ 2 ( 3 γ 1 + − μ a ) δ a . (115) Using
δ a = γ 1 + − γ 1 − − G , μ a = γ 1 − + G , \delta_a=\gamma_1^+ - \gamma_1^- - G,\qquad\mu_a=\gamma_1^-+G, δ a = γ 1 + − γ 1 − − G , μ a = γ 1 − + G , we obtain
χ a = γ 1 + − γ 1 − − G γ 2 ( 3 γ 1 + − γ 1 − − G ) . (116) \chi_a=\frac{\gamma_1^+ - \gamma_1^- - G}{\gamma_2(3\gamma_1^+ - \gamma_1^- - G)}. \tag*{(116)} χ a = γ 2 ( 3 γ 1 + − γ 1 − − G ) γ 1 + − γ 1 − − G . (116) The same calculation for oscillator b b b
χ b = δ b Γ 2 ( 3 Γ 1 + − μ b ) . (117) \chi_b=\frac{\delta_b}{\Gamma_2(3\Gamma_1^+ - \mu_b)}. \tag*{(117)} χ b = Γ 2 ( 3 Γ 1 + − μ b ) δ b . (117) Equivalently,
χ b = Γ 1 + − Γ 1 − − G Γ 2 ( 3 Γ 1 + − Γ 1 − − G ) . (118) \chi_b=\frac{\Gamma_1^+ - \Gamma_1^- - G}{\Gamma_2(3\Gamma_1^+ - \Gamma_1^- - G)}. \tag*{(118)} χ b = Γ 2 ( 3 Γ 1 + − Γ 1 − − G ) Γ 1 + − Γ 1 − − G . (118) The explicit formulas (C43) and (C45) require
δ a > 0 , δ b > 0 , (119) \delta_a>0,\qquad\delta_b>0, \tag*{(119)} δ a > 0 , δ b > 0 , (119) together with the large-occupation conditions
δ a γ 2 ≫ 1 , δ b Γ 2 ≫ 1. (120) \frac{\delta_a}{\gamma_2}\gg1,\qquad\frac{\delta_b}{\Gamma_2}\gg1. \tag*{(120)} γ 2 δ a ≫ 1 , Γ 2 δ b ≫ 1. (120) The above results should not be extrapolated to the critical point δ a = δ b = 0 \delta_a=\delta_b=0 δ a = δ b = 0
Weak coupling between the oscillators
The exact first-moment equations are
d d t ⟨ a ⟩ = Ω + δ a ⟨ a ⟩ + ( g − G ) ⟨ b ⟩ − 2 γ 2 ⟨ a † a 2 ⟩ , (121) \frac{d}{dt}\langle a\rangle=\Omega+\delta_a\langle a\rangle+(g-G)\langle b\rangle-2\gamma_2\langle a^\dagger a^2\rangle, \tag*{(121)} d t d ⟨ a ⟩ = Ω + δ a ⟨ a ⟩ + ( g − G ) ⟨ b ⟩ − 2 γ 2 ⟨ a † a 2 ⟩ , (121) d d t ⟨ b ⟩ = δ b ⟨ b ⟩ − ( g + G ) ⟨ a ⟩ − 2 Γ 2 ⟨ b † b 2 ⟩ . (122) \frac{d}{dt}\langle b\rangle=\delta_b\langle b\rangle-(g+G)\langle a\rangle-2\Gamma_2\langle b^\dagger b^2\rangle. \tag*{(122)} d t d ⟨ b ⟩ = δ b ⟨ b ⟩ − ( g + G ) ⟨ a ⟩ − 2 Γ 2 ⟨ b † b 2 ⟩ . (122) Within the local linear-response approximation, oscillator a a a
F a = Ω + ( g − G ) ⟨ b ⟩ , (123) F_a=\Omega+(g-G)\langle b\rangle, \tag*{(123)} F a = Ω + ( g − G ) ⟨ b ⟩ , (123) whereas oscillator b b b
F b = − ( g + G ) ⟨ a ⟩ . (124) F_b=-(g+G)\langle a\rangle. \tag*{(124)} F b = − ( g + G ) ⟨ a ⟩ . (124) Using the explicit local susceptibilities, we obtain
⟨ a ⟩ = χ a [ Ω + ( g − G ) ⟨ b ⟩ ] , (125) \langle a\rangle=\chi_a\left[\Omega+(g-G)\langle b\rangle\right], \tag*{(125)} ⟨ a ⟩ = χ a [ Ω + ( g − G ) ⟨ b ⟩ ] , (125) ⟨ b ⟩ = − χ b ( g + G ) ⟨ a ⟩ . (126) \langle b\rangle=-\chi_b(g+G)\langle a\rangle. \tag*{(126)} ⟨ b ⟩ = − χ b ( g + G ) ⟨ a ⟩ . (126) Solving the two linear equations yields
⟨ a ⟩ s s = χ a Ω 1 + ( g 2 − G 2 ) χ a χ b (127) \langle a\rangle_{\mathrm{ss}}=\frac{\chi_a\Omega}{1+(g^2-G^2)\chi_a\chi_b} \tag*{(127)} ⟨ a ⟩ ss = 1 + ( g 2 − G 2 ) χ a χ b χ a Ω (127) ⟨ b ⟩ s s = − ( g + G ) χ a χ b Ω 1 + ( g 2 − G 2 ) χ a χ b . (128) \langle b\rangle_{\mathrm{ss}}=-\frac{(g+G)\chi_a\chi_b\Omega}{1+(g^2-G^2)\chi_a\chi_b}. \tag*{(128)} ⟨ b ⟩ ss = − 1 + ( g 2 − G 2 ) χ a χ b ( g + G ) χ a χ b Ω . (128) Substituting Eqs. (C43) and (C45) into the above equations, we obtain
⟨ a ⟩ s s = γ 1 + − γ 1 − − G γ 2 ( 3 γ 1 + − γ 1 − − G ) Ω 1 + ( g 2 − G 2 ) ( γ 1 + − γ 1 − − G ) ( Γ 1 + − Γ 1 − − G ) γ 2 Γ 2 ( 3 γ 1 + − γ 1 − − G ) ( 3 Γ 1 + − Γ 1 − − G ) . (129) \langle a\rangle_{\mathrm{ss}}=
\frac{\dfrac{\gamma_1^+ - \gamma_1^- - G}{\gamma_2(3\gamma_1^+ - \gamma_1^- - G)}\Omega}
{1+\dfrac{(g^2-G^2)(\gamma_1^+ - \gamma_1^- - G)(\Gamma_1^+ - \Gamma_1^- - G)}
{\gamma_2\Gamma_2(3\gamma_1^+ - \gamma_1^- - G)(3\Gamma_1^+ - \Gamma_1^- - G)}}. \tag*{(129)} ⟨ a ⟩ ss = 1 + γ 2 Γ 2 ( 3 γ 1 + − γ 1 − − G ) ( 3 Γ 1 + − Γ 1 − − G ) ( g 2 − G 2 ) ( γ 1 + − γ 1 − − G ) ( Γ 1 + − Γ 1 − − G ) γ 2 ( 3 γ 1 + − γ 1 − − G ) γ 1 + − γ 1 − − G Ω . (129) ⟨ b ⟩ s s = − ( g + G ) ( γ 1 + − γ 1 − − G ) ( Γ 1 + − Γ 1 − − G ) Ω γ 2 Γ 2 ( 3 γ 1 + − γ 1 − − G ) ( 3 Γ 1 + − Γ 1 − − G ) [ 1 + ( g 2 − G 2 ) ( γ 1 + − γ 1 − − G ) ( Γ 1 + − Γ 1 − − G ) γ 2 Γ 2 ( 3 γ 1 + − γ 1 − − G ) ( 3 Γ 1 + − Γ 1 − − G ) ] − 1 . (130) \langle b\rangle_{\mathrm{ss}}=
-\frac{(g+G)(\gamma_1^+ - \gamma_1^- - G)(\Gamma_1^+ - \Gamma_1^- - G)\Omega}
{\gamma_2\Gamma_2(3\gamma_1^+ - \gamma_1^- - G)(3\Gamma_1^+ - \Gamma_1^- - G)}
\left[
1+\frac{(g^2-G^2)(\gamma_1^+ - \gamma_1^- - G)(\Gamma_1^+ - \Gamma_1^- - G)}
{\gamma_2\Gamma_2(3\gamma_1^+ - \gamma_1^- - G)(3\Gamma_1^+ - \Gamma_1^- - G)}
\right]^{-1}. \tag*{(130)} ⟨ b ⟩ ss = − γ 2 Γ 2 ( 3 γ 1 + − γ 1 − − G ) ( 3 Γ 1 + − Γ 1 − − G ) ( g + G ) ( γ 1 + − γ 1 − − G ) ( Γ 1 + − Γ 1 − − G ) Ω [ 1 + γ 2 Γ 2 ( 3 γ 1 + − γ 1 − − G ) ( 3 Γ 1 + − Γ 1 − − G ) ( g 2 − G 2 ) ( γ 1 + − γ 1 − − G ) ( Γ 1 + − Γ 1 − − G ) ] − 1 . (130) When the coupling is completely nonreciprocal (G = g G=g G = g
⟨ a ⟩ s s ≃ γ 1 + − γ 1 − − G γ 2 ( 3 γ 1 + − γ 1 − − G ) Ω (131) \langle a\rangle_{\mathrm{ss}}\simeq
\frac{\gamma_1^+ - \gamma_1^- - G}{\gamma_2(3\gamma_1^+ - \gamma_1^- - G)}\Omega \tag*{(131)} ⟨ a ⟩ ss ≃ γ 2 ( 3 γ 1 + − γ 1 − − G ) γ 1 + − γ 1 − − G Ω (131) and
⟨ b ⟩ s s ≃ − ( g + G ) ( γ 1 + − γ 1 − − G ) ( Γ 1 + − Γ 1 − − G ) γ 2 Γ 2 ( 3 γ 1 + − γ 1 − − G ) ( 3 Γ 1 + − Γ 1 − − G ) Ω . (132) \langle b\rangle_{\mathrm{ss}}\simeq
-\frac{(g+G)(\gamma_1^+ - \gamma_1^- - G)(\Gamma_1^+ - \Gamma_1^- - G)}
{\gamma_2\Gamma_2(3\gamma_1^+ - \gamma_1^- - G)(3\Gamma_1^+ - \Gamma_1^- - G)}\Omega. \tag*{(132)} ⟨ b ⟩ ss ≃ − γ 2 Γ 2 ( 3 γ 1 + − γ 1 − − G ) ( 3 Γ 1 + − Γ 1 − − G ) ( g + G ) ( γ 1 + − γ 1 − − G ) ( Γ 1 + − Γ 1 − − G ) Ω. (132) Occupation-number corrections
The exact number equations are
d d t ⟨ n a ⟩ = Ω ⟨ a + a † ⟩ + 2 δ a ⟨ n a ⟩ + 2 γ 1 + − 4 γ 2 ⟨ n a ( n a − 1 ) ⟩ + ( g − G ) ⟨ a † b + a b † ⟩ , (133) \frac{d}{dt}\langle n_a\rangle
=\Omega\langle a+a^\dagger\rangle+2\delta_a\langle n_a\rangle+2\gamma_1^+
-4\gamma_2\langle n_a(n_a-1)\rangle+(g-G)\langle a^\dagger b+ab^\dagger\rangle, \tag*{(133)} d t d ⟨ n a ⟩ = Ω ⟨ a + a † ⟩ + 2 δ a ⟨ n a ⟩ + 2 γ 1 + − 4 γ 2 ⟨ n a ( n a − 1 )⟩ + ( g − G ) ⟨ a † b + a b † ⟩ , (133) d d t ⟨ n b ⟩ = 2 δ b ⟨ n b ⟩ + 2 Γ 1 + − 4 Γ 2 ⟨ n b ( n b − 1 ) ⟩ − ( g + G ) ⟨ a † b + a b † ⟩ . (134) \frac{d}{dt}\langle n_b\rangle
=2\delta_b\langle n_b\rangle+2\Gamma_1^+
-4\Gamma_2\langle n_b(n_b-1)\rangle-(g+G)\langle a^\dagger b+ab^\dagger\rangle. \tag*{(134)} d t d ⟨ n b ⟩ = 2 δ b ⟨ n b ⟩ + 2 Γ 1 + − 4 Γ 2 ⟨ n b ( n b − 1 )⟩ − ( g + G ) ⟨ a † b + a b † ⟩ . (134) The interoscillator coherence is of first order in the coupling. It enters the number equations multiplied by a second coupling factor. Therefore,
Δ n a , Δ n b = O ( g 2 , G 2 , g G ) . (135) \Delta n_a,\Delta n_b=O(g^2,G^2,gG). \tag*{(135)} Δ n a , Δ n b = O ( g 2 , G 2 , g G ) . (135) The coherent drive changes the occupations only at order Ω 2 \Omega^2 Ω 2
⟨ a † a ⟩ s s = γ 1 + − γ 1 − − G 2 γ 2 + O ( g 2 , G 2 , g G , Ω 2 ) (136) \langle a^\dagger a\rangle_{\mathrm{ss}}=\frac{\gamma_1^+-\gamma_1^--G}{2\gamma_2}+O(g^2,G^2,gG,\Omega^2) \tag*{(136)} ⟨ a † a ⟩ ss = 2 γ 2 γ 1 + − γ 1 − − G + O ( g 2 , G 2 , g G , Ω 2 ) (136) ⟨ b † b ⟩ s s = Γ 1 + − Γ 1 − − G 2 Γ 2 + O ( g 2 , G 2 , g G , Ω 2 ) (137) \langle b^\dagger b\rangle_{\mathrm{ss}}=\frac{\Gamma_1^+-\Gamma_1^--G}{2\Gamma_2}+O(g^2,G^2,gG,\Omega^2) \tag*{(137)} ⟨ b † b ⟩ ss = 2 Γ 2 Γ 1 + − Γ 1 − − G + O ( g 2 , G 2 , g G , Ω 2 ) (137) critical regime
This Appendix studies the response of mode b b b Γ 2 \Gamma_2 Γ 2
We consider two bosonic modes satisfying
[ a , a † ] = [ b , b † ] = 1 , [ a , b ] = [ a , b † ] = 0 , (138) [a,a^\dagger]=[b,b^\dagger]=1,\qquad[a,b]=[a,b^\dagger]=0, \tag*{(138)} [ a , a † ] = [ b , b † ] = 1 , [ a , b ] = [ a , b † ] = 0 , (138) with density operator evolving according to
ρ ˙ = [ X , ρ ] + γ 1 + D [ a † ] ρ + γ 1 − D [ a ] ρ + γ 2 D [ a 2 ] ρ + Γ 1 + D [ b † ] ρ + Γ 1 − D [ b ] ρ + Γ 2 D [ b 2 ] ρ + G D [ a + b ] ρ , (139) \dot{\rho}=[X,\rho]+\gamma_1^+\mathcal{D}[a^\dagger]\rho+\gamma_1^-\mathcal{D}[a]\rho+\gamma_2\mathcal{D}[a^2]\rho+\Gamma_1^+\mathcal{D}[b^\dagger]\rho+\Gamma_1^-\mathcal{D}[b]\rho+\Gamma_2\mathcal{D}[b^2]\rho+G\mathcal{D}[a+b]\rho, \tag*{(139)} ρ ˙ = [ X , ρ ] + γ 1 + D [ a † ] ρ + γ 1 − D [ a ] ρ + γ 2 D [ a 2 ] ρ + Γ 1 + D [ b † ] ρ + Γ 1 − D [ b ] ρ + Γ 2 D [ b 2 ] ρ + G D [ a + b ] ρ , (139) where
X = Ω ( a † − a ) + g ( a † b − a b † ) , (140) X=\Omega(a^\dagger-a)+g(a^\dagger b-ab^\dagger), \tag*{(140)} X = Ω ( a † − a ) + g ( a † b − a b † ) , (140) D [ L ] ρ = 2 L ρ L † − { L † L , ρ } . (141) \mathcal{D}[L]\rho=2L\rho L^\dagger-\{L^\dagger L,\rho\}. \tag*{(141)} D [ L ] ρ = 2 L ρ L † − { L † L , ρ } . (141) We take Ω ∈ R \Omega\in\mathbb{R} Ω ∈ R g > 0 g>0 g > 0 X † = − X X^\dagger=-X X † = − X H = i X H=iX H = i X [ X , ρ ] = − i [ H , ρ ] [X,\rho]=-i[H,\rho] [ X , ρ ] = − i [ H , ρ ]
The linear rates obey the critical condition
γ 1 + − γ 1 − = Γ 1 + − Γ 1 − = G . (142) \gamma_1^+-\gamma_1^-=\Gamma_1^+-\Gamma_1^-=G. \tag*{(142)} γ 1 + − γ 1 − = Γ 1 + − Γ 1 − = G . (142) The reciprocal and nonreciprocal cases are compared at fixed γ 1 + , Γ 1 + , γ 2 , g , Ω \gamma_1^+,\Gamma_1^+,\gamma_2,g,\Omega γ 1 + , Γ 1 + , γ 2 , g , Ω Γ 2 \Gamma_2 Γ 2 G = g G=g G = g
0 < g ≤ min ( γ 1 + , Γ 1 + ) . (143) 0<g\leq\min(\gamma_1^+,\Gamma_1^+). \tag*{(143)} 0 < g ≤ min ( γ 1 + , Γ 1 + ) . (143) Nonreciprocal coupling at G = g G=g G = g
For G = g G=g G = g
ρ ˙ = ( L a + L b + L c a s ) ρ , (144) \dot{\rho}=(\mathcal{L}_a+\mathcal{L}_b+\mathcal{L}_{\mathrm{cas}})\rho, \tag*{(144)} ρ ˙ = ( L a + L b + L cas ) ρ , (144) where
L a ρ = [ Ω ( a † − a ) , ρ ] + γ 1 + ( D [ a ] + D [ a † ] ) ρ + γ 2 D [ a 2 ] ρ , (145) \mathcal{L}_a\rho=[\Omega(a^\dagger-a),\rho]+\gamma_1^+\bigl(\mathcal{D}[a]+\mathcal{D}[a^\dagger]\bigr)\rho+\gamma_2\mathcal{D}[a^2]\rho, \tag*{(145)} L a ρ = [ Ω ( a † − a ) , ρ ] + γ 1 + ( D [ a ] + D [ a † ] ) ρ + γ 2 D [ a 2 ] ρ , (145) L b ρ = Γ 1 + ( D [ b ] + D [ b † ] ) ρ + Γ 2 D [ b 2 ] ρ , (146) \mathcal{L}_b\rho=\Gamma_1^+\bigl(\mathcal{D}[b]+\mathcal{D}[b^\dagger]\bigr)\rho+\Gamma_2\mathcal{D}[b^2]\rho, \tag*{(146)} L b ρ = Γ 1 + ( D [ b ] + D [ b † ] ) ρ + Γ 2 D [ b 2 ] ρ , (146) L c a s ρ = − 2 g ( [ b † , a ρ ] + [ ρ a † , b ] ) . (147) \mathcal{L}_{\mathrm{cas}}\rho=-2g\bigl([b^\dagger,a\rho]+[\rho a^\dagger,b]\bigr). \tag*{(147)} L cas ρ = − 2 g ( [ b † , a ρ ] + [ ρ a † , b ] ) . (147) The cascade term follows from
[ g ( a † b − a b † ) , ρ ] + g ( D [ a + b ] − D [ a ] − D [ b ] ) ρ = − 2 g ( [ b † , a ρ ] + [ ρ a † , b ] ) . (148) [g(a^\dagger b-ab^\dagger),\rho]+g\bigl(\mathcal{D}[a+b]-\mathcal{D}[a]-\mathcal{D}[b]\bigr)\rho=-2g\bigl([b^\dagger,a\rho]+[\rho a^\dagger,b]\bigr). \tag*{(148)} [ g ( a † b − a b † ) , ρ ] + g ( D [ a + b ] − D [ a ] − D [ b ] ) ρ = − 2 g ( [ b † , a ρ ] + [ ρ a † , b ] ) . (148) Tracing over mode b b b
ρ ˙ a = L a ρ a . (149) \dot{\rho}_a=\mathcal{L}_a\rho_a. \tag*{(149)} ρ ˙ a = L a ρ a . (149) Consequently, the stationary state of mode a a a ρ a , s s \rho_{a,\mathrm{ss}} ρ a , ss Γ 2 \Gamma_2 Γ 2
α a ( Ω ) ≡ Tr ( a ρ a , s s ) . (150) \alpha_a(\Omega)\equiv\operatorname{Tr}(a\rho_{a,\mathrm{ss}}). \tag*{(150)} α a ( Ω ) ≡ Tr ( a ρ a , ss ) . (150) The mean force acting on mode b b b
F ( Ω ) = − 2 g α a ( Ω ) . (151) F(\Omega)=-2g\alpha_a(\Omega). \tag*{(151)} F ( Ω ) = − 2 g α a ( Ω ) . (151) For a real drive and a symmetry-preserving stationary state, α a \alpha_a α a F F F
The exact mean equation for mode b b b
d d t ⟨ b ⟩ = F − 2 Γ 2 ⟨ b † b 2 ⟩ . (152) \frac{d}{dt}\langle b\rangle=F-2\Gamma_2\langle b^\dagger b^2\rangle. \tag*{(152)} d t d ⟨ b ⟩ = F − 2 Γ 2 ⟨ b † b 2 ⟩ . (152) In particular, the force F F F Γ 2 \Gamma_2 Γ 2
Large-occupation estimate for the upstream (a a a
An estimate of α a \alpha_a α a
ρ ˙ = γ 1 + ( D [ a ] + D [ a † ] ) ρ + γ 2 D [ a 2 ] ρ . (153) \dot{\rho}=\gamma_1^+\bigl(\mathcal{D}[a]+\mathcal{D}[a^\dagger]\bigr)\rho+\gamma_2\mathcal{D}[a^2]\rho. \tag*{(153)} ρ ˙ = γ 1 + ( D [ a ] + D [ a † ] ) ρ + γ 2 D [ a 2 ] ρ . (153) The number-state transition rates are
n → n + 1 : 2 γ 1 + ( n + 1 ) , (154) n\to n+1:\quad2\gamma_1^+(n+1), \tag*{(154)} n → n + 1 : 2 γ 1 + ( n + 1 ) , (154) n → n − 1 : 2 γ 1 + n , (155) n\to n-1:\quad2\gamma_1^+n, \tag*{(155)} n → n − 1 : 2 γ 1 + n , (155) n → n − 2 : 2 γ 2 n ( n − 1 ) . (156) n\to n-2:\quad2\gamma_2n(n-1). \tag*{(156)} n → n − 2 : 2 γ 2 n ( n − 1 ) . (156) For γ 2 / γ 1 + ≪ 1 \gamma_2/\gamma_1^+\ll1 γ 2 / γ 1 + ≪ 1
n ∼ γ 1 + γ 2 ≫ 1. (157) n\sim\sqrt{\frac{\gamma_1^+}{\gamma_2}}\gg1. \tag*{(157)} n ∼ γ 2 γ 1 + ≫ 1. (157) The first two Kramers–Moyal coefficients are
A ( n ) = 2 γ 1 + − 4 γ 2 n 2 + O ( n ) , (158) A(n)=2\gamma_1^+-4\gamma_2n^2+O(n), \tag*{(158)} A ( n ) = 2 γ 1 + − 4 γ 2 n 2 + O ( n ) , (158) B ( n ) = 4 γ 1 + n + 8 γ 2 n 2 + O ( n ) . (159) B(n)=4\gamma_1^+n+8\gamma_2n^2+O(n). \tag*{(159)} B ( n ) = 4 γ 1 + n + 8 γ 2 n 2 + O ( n ) . (159) Here the Fokker–Planck equation is written as
∂ t p ( n , t ) = − ∂ n [ A ( n ) p ] + 1 2 ∂ n 2 [ B ( n ) p ] . (160) \partial_t p(n,t)=-\partial_n[A(n)p]+\frac{1}{2}\partial_n^2[B(n)p]. \tag*{(160)} ∂ t p ( n , t ) = − ∂ n [ A ( n ) p ] + 2 1 ∂ n 2 [ B ( n ) p ] . (160) At large occupation, the leading terms are
A ( n ) ≃ 2 γ 1 + − 4 γ 2 n 2 , B ( n ) ≃ 4 γ 1 + n . (161) A(n)\simeq2\gamma_1^+-4\gamma_2n^2,\qquad B(n)\simeq4\gamma_1^+n. \tag*{(161)} A ( n ) ≃ 2 γ 1 + − 4 γ 2 n 2 , B ( n ) ≃ 4 γ 1 + n . (161) The zero-current stationary solution is then given by
p 0 ( n ) ≃ 2 γ 2 π γ 1 + exp ( − γ 2 n 2 γ 1 + ) , n ≥ 0. (162) p_0(n) \simeq2\sqrt{\frac{\gamma_2}{\pi\gamma_1^+}}\exp\left(-\frac{\gamma_2 n^2}{\gamma_1^+}\right), \qquad n \ge0. \tag*{(162)} p 0 ( n ) ≃ 2 π γ 1 + γ 2 exp ( − γ 1 + γ 2 n 2 ) , n ≥ 0. (162) Consequently,
⟨ n ⟩ 0 ≃ γ 1 + γ 2 , (163) \langle n\rangle_0 \simeq\sqrt{\frac{\gamma_1^+}{\gamma_2}}, \tag*{(163)} ⟨ n ⟩ 0 ≃ γ 2 γ 1 + , (163) ⟨ n 2 ⟩ 0 ≃ γ 1 + 2 γ 2 . (164) \langle n^2\rangle_0 \simeq\frac{\gamma_1^+}{2\gamma_2}. \tag*{(164)} ⟨ n 2 ⟩ 0 ≃ 2 γ 2 γ 1 + . (164) The exact stationary number balance instead yields
⟨ n ( n − 1 ) ⟩ 0 = γ 1 + 2 γ 2 . (165) \langle n(n-1)\rangle_0 = \frac{\gamma_1^+}{2\gamma_2}. \tag*{(165)} ⟨ n ( n − 1 ) ⟩ 0 = 2 γ 2 γ 1 + . (165) The difference between ⟨ n 2 ⟩ 0 \langle n^2\rangle_0 ⟨ n 2 ⟩ 0 ⟨ n ( n − 1 ) ⟩ 0 \langle n(n-1)\rangle_0 ⟨ n ( n − 1 ) ⟩ 0
For a weak real drive, the leading large-occupation equation of motion for the amplitude reads
d α = ( Ω − 2 γ 2 ∣ α ∣ 2 α ) d t + κ a ( d W x + i d W y ) , (166) d\alpha= \left(\Omega-2\gamma_2|\alpha|^2\alpha\right)dt+\sqrt{\kappa_a}\left(dW_x+i\,dW_y\right), \tag*{(166)} d α = ( Ω − 2 γ 2 ∣ α ∣ 2 α ) d t + κ a ( d W x + i d W y ) , (166) where d W x , d W y dW_x,dW_y d W x , d W y E [ d W i d W j ] = δ i j d t \mathbb{E}[dW_i\,dW_j]=\delta_{ij}\,dt E [ d W i d W j ] = δ ij d t α = x + i y \alpha=x+iy α = x + i y
P Ω ( x , y ) = 1 Z Ω exp [ − γ 2 γ 1 + ( x 2 + y 2 ) 2 + 2 Ω x γ 1 + ] . (167) P_\Omega(x,y)=\frac{1}{Z_\Omega}\exp\left[-\frac{\gamma_2}{\gamma_1^+}(x^2+y^2)^2+\frac{2\Omega x}{\gamma_1^+}\right]. \tag*{(167)} P Ω ( x , y ) = Z Ω 1 exp [ − γ 1 + γ 2 ( x 2 + y 2 ) 2 + γ 1 + 2Ω x ] . (167) At zero drive, rotational symmetry gives ⟨ x 2 ⟩ 0 = 1 2 ⟨ ∣ α ∣ 2 ⟩ 0 \langle x^2\rangle_0=\frac{1}{2}\langle|\alpha|^2\rangle_0 ⟨ x 2 ⟩ 0 = 2 1 ⟨ ∣ α ∣ 2 ⟩ 0 Ω \Omega Ω
α a = ⟨ a ⟩ ≃ 2 Ω γ 1 + ⟨ x 2 ⟩ 0 ≃ Ω π γ 1 + γ 2 . (168) \alpha_a=\langle a\rangle\simeq\frac{2\Omega}{\gamma_1^+}\langle x^2\rangle_0
\simeq\frac{\Omega}{\sqrt{\pi\gamma_1^+\gamma_2}}. \tag*{(168)} α a = ⟨ a ⟩ ≃ γ 1 + 2Ω ⟨ x 2 ⟩ 0 ≃ π γ 1 + γ 2 Ω . (168) Define the upstream zero-drive susceptibility by
χ a ≡ ∂ α a ∂ Ω ∣ Ω = 0 . (169) \chi_a\equiv\left.\frac{\partial\alpha_a}{\partial\Omega}\right|_{\Omega=0}. \tag*{(169)} χ a ≡ ∂ Ω ∂ α a Ω = 0 . (169) Then
N a ≡ ⟨ a † a ⟩ ≃ γ 1 + π γ 2 , χ a ≃ 1 π γ 1 + γ 2 . (170) N_a\equiv\langle a^\dagger a\rangle\simeq\sqrt{\frac{\gamma_1^+}{\pi\gamma_2}},\qquad
\chi_a\simeq\frac{1}{\sqrt{\pi\gamma_1^+\gamma_2}}. \tag*{(170)} N a ≡ ⟨ a † a ⟩ ≃ π γ 2 γ 1 + , χ a ≃ π γ 1 + γ 2 1 . (170) The weak-drive condition for this estimate is
∣ Ω ∣ ≪ ( γ 1 + ) 3 / 4 γ 2 1 / 4 . (171) |\Omega|\ll(\gamma_1^+)^{3/4}\gamma_2^{1/4}. \tag*{(171)} ∣Ω∣ ≪ ( γ 1 + ) 3/4 γ 2 1/4 . (171) Exact zero-drive heating rate at Γ 2 = 0 \Gamma_2=0 Γ 2 = 0
At Ω = 0 \Omega=0 Ω = 0 Γ 2 = 0 \Gamma_2=0 Γ 2 = 0 ρ a , 0 \rho_{a,0} ρ a , 0
M ( t ) = Tr b [ b ρ ( t ) ] . (172) M(t)=\operatorname{Tr}_b[b\rho(t)]. \tag*{(172)} M ( t ) = Tr b [ b ρ ( t )] . (172) Once the upstream marginal has reached ρ a , 0 \rho_{a,0} ρ a , 0
M ˙ = L a , 0 M − 2 g a ρ a , 0 . (173) \dot{M}=\mathcal{L}_{a,0}M-2g\,a\rho_{a,0}. \tag*{(173)} M ˙ = L a , 0 M − 2 g a ρ a , 0 . (173) Since Tr ( a ρ a , 0 ) = 0 \operatorname{Tr}(a\rho_{a,0})=0 Tr ( a ρ a , 0 ) = 0
M ⊥ = − 2 g ∫ 0 ∞ e L a , 0 t ( a ρ a , 0 ) d t . (174) M_\perp=-2g\int_0^\infty e^{\mathcal{L}_{a,0}t}(a\rho_{a,0})\,dt. \tag*{(174)} M ⊥ = − 2 g ∫ 0 ∞ e L a , 0 t ( a ρ a , 0 ) d t . (174) The intermode coherence is
⟨ a † b ⟩ = Tr ( a † M ⊥ ) = − 2 g ∫ 0 ∞ ⟨ a † ( t ) a ( 0 ) ⟩ 0 d t . (175) \begin{aligned}
\langle a^\dagger b\rangle=\operatorname{Tr}(a^\dagger M_\perp) \\
&=-2g\int_0^\infty\langle a^\dagger(t)a(0)\rangle_0\,dt. \tag*{(175)}
\end{aligned} ⟨ a † b ⟩ = Tr ( a † M ⊥ ) = − 2 g ∫ 0 ∞ ⟨ a † ( t ) a ( 0 ) ⟩ 0 d t . (175) where the last equality follows from the quantum regression theorem. Any stationary component proportional to ρ a , 0 \rho_{a,0} ρ a , 0 ⟨ a † ⟩ 0 = 0 \langle a^\dagger\rangle_0=0 ⟨ a † ⟩ 0 = 0
At G = g G=g G = g Γ 2 = 0 \Gamma_2=0 Γ 2 = 0 b b b
d N b d t = − 2 g ( ⟨ a † b ⟩ + ⟨ a b † ⟩ ) + 2 Γ 1 + , (176) \frac{dN_b}{dt}=-2g\left(\langle a^\dagger b\rangle+\langle ab^\dagger\rangle\right)+2\Gamma_1^+, \tag*{(176)} d t d N b = − 2 g ( ⟨ a † b ⟩ + ⟨ a b † ⟩ ) + 2 Γ 1 + , (176) where N b = ⟨ b † b ⟩ N_b=\langle b^\dagger b\rangle N b = ⟨ b † b ⟩
lim t → ∞ d N b d t = 2 Γ 1 + + 8 g 2 Re ∫ 0 ∞ ⟨ a † ( t ) a ( 0 ) ⟩ 0 d t \lim_{t\to\infty}\frac{dN_b}{dt}
=2\Gamma_1^++8g^2\operatorname{Re}\int_0^\infty\langle a^\dagger(t)a(0)\rangle_0\,dt t → ∞ lim d t d N b = 2 Γ 1 + + 8 g 2 Re ∫ 0 ∞ ⟨ a † ( t ) a ( 0 ) ⟩ 0 d t ≡ 2 K . (177) \equiv2K. \tag*{(177)} ≡ 2 K . (177) The exact zero-drive heating coefficient is therefore
K = κ b + 4 g 2 Re ∫ 0 ∞ ⟨ a † ( t ) a ( 0 ) ⟩ 0 d t (178) K=\kappa_b+4g^2\operatorname{Re}\int_0^\infty\langle a^\dagger(t)a(0)\rangle_0\,dt \tag*{(178)} K = κ b + 4 g 2 Re ∫ 0 ∞ ⟨ a † ( t ) a ( 0 ) ⟩ 0 d t (178) or, using κ b = Γ 1 + \kappa_b=\Gamma_1^+ κ b = Γ 1 +
K = Γ 1 + + 4 g 2 Re ∫ 0 ∞ ⟨ a † ( t ) a ( 0 ) ⟩ 0 d t . (179) K=\Gamma_1^++4g^2\operatorname{Re}\int_0^\infty\langle a^\dagger(t)a(0)\rangle_0\,dt. \tag*{(179)} K = Γ 1 + + 4 g 2 Re ∫ 0 ∞ ⟨ a † ( t ) a ( 0 ) ⟩ 0 d t . (179) Eq. (D42) is an exact expression within the quantum master equation. For the nonlinear upstream problem containing γ 2 D [ a 2 ] \gamma_2\mathcal{D}[a^2] γ 2 D [ a 2 ] γ 1 + \gamma_1^+ γ 1 + γ 2 \gamma_2 γ 2 g g g Γ 1 + \Gamma_1^+ Γ 1 +
For numerical evaluation, solve the Poisson equation
L a , 0 Y = − a ρ a , 0 , Tr Y = 0 , (180) \mathcal{L}_{a,0}Y=-a\rho_{a,0},\qquad\operatorname{Tr}Y=0, \tag*{(180)} L a , 0 Y = − a ρ a , 0 , Tr Y = 0 , (180) on the trace-zero subspace. Then
K = Γ 1 + + 4 g 2 Re Tr ( a † Y ) . (181) K=\Gamma_1^++4g^2\operatorname{Re}\operatorname{Tr}(a^\dagger Y). \tag*{(181)} K = Γ 1 + + 4 g 2 Re Tr ( a † Y ) . (181) This gives the exact zero-drive long-time heating rate. It does not, by itself, imply an exact white-noise description of the downstream dynamics.
Effective downstream stochastic description
To obtain a local stochastic description at small positive Γ 2 \Gamma_2 Γ 2 τ a \tau_a τ a
d β = ( F − 2 Γ 2 ∣ β ∣ 2 β ) d t + K ( d W x + i d W y ) , (182) d\beta= (F - 2\Gamma_2|\beta|^2\beta)\,dt + \sqrt{K}\,(dW_x + i\,dW_y), \tag*{(182)} d β = ( F − 2 Γ 2 ∣ β ∣ 2 β ) d t + K ( d W x + i d W y ) , (182) where β = x + i y \beta= x+iy β = x + i y
d d t E [ ∣ β ∣ 2 ] = 2 K . (183) \frac{d}{dt}\mathbb{E}[|\beta|^2] = 2K. \tag*{(183)} d t d E [ ∣ β ∣ 2 ] = 2 K . (183) Eq. (D45) neglects operator-ordering corrections and noise generated by two-photon loss. The latter produces quadrature diffusion of order Γ 2 ∣ β ∣ 2 \Gamma_2|\beta|^2 Γ 2 ∣ β ∣ 2 K K K
At finite drive, the coarse-grained noise may be anisotropic. More generally, with
r = ( x , y ) T , \mathbf{r} = (x,y)^{\mathsf{T}}, r = ( x , y ) T , one may write
d r = [ F − 2 Γ 2 ∣ r ∣ 2 r ] d t + B Ω d W , (184) d\mathbf{r} = [\mathbf{F} - 2\Gamma_2|\mathbf{r}|^2\mathbf{r}]\,dt + B_\Omega\,d\mathbf{W}, \tag*{(184)} d r = [ F − 2 Γ 2 ∣ r ∣ 2 r ] d t + B Ω d W , (184) where
F = ( Re F , Im F ) T , Q Ω = B Ω B Ω T . (185) \mathbf{F} = (\operatorname{Re}F,\operatorname{Im}F)^{\mathsf{T}},\qquad Q_\Omega= B_\Omega B_\Omega^{\mathsf{T}}. \tag*{(185)} F = ( Re F , Im F ) T , Q Ω = B Ω B Ω T . (185) At zero drive, phase symmetry gives Q 0 = K I 2 Q_0 = KI_2 Q 0 = K I 2 K K K
For real F F F
V ( x , y ) = Γ 2 2 ( x 2 + y 2 ) 2 − F x . (186) V(x,y) = \frac{\Gamma_2}{2}(x^2+y^2)^2 - Fx. \tag*{(186)} V ( x , y ) = 2 Γ 2 ( x 2 + y 2 ) 2 − F x . (186) The diffusion coefficient of each real quadrature in the Fokker–Planck equation is
D q u a d = K 2 . (187) D_{\mathrm{quad}} = \frac{K}{2}. \tag*{(187)} D quad = 2 K . (187) The stationary density is therefore
P F ( x , y ) = 1 Z F exp [ − Γ 2 K ( x 2 + y 2 ) 2 + 2 F x K ] . (188) P_F(x,y) = \frac{1}{Z_F}\exp\left[-\frac{\Gamma_2}{K}(x^2+y^2)^2+\frac{2Fx}{K}\right]. \tag*{(188)} P F ( x , y ) = Z F 1 exp [ − K Γ 2 ( x 2 + y 2 ) 2 + K 2 F x ] . (188) Define the dimensionless variables
ℓ = ( K Γ 2 ) 1 / 4 , h = 2 ∣ F ∣ K 3 / 4 Γ 2 1 / 4 . (189) \ell= \left(\frac{K}{\Gamma_2}\right)^{1/4},\qquad h = \frac{2|F|}{K^{3/4}\Gamma_2^{1/4}}. \tag*{(189)} ℓ = ( Γ 2 K ) 1/4 , h = K 3/4 Γ 2 1/4 2∣ F ∣ . (189) Integration in polar coordinates gives
∣ ⟨ b ⟩ ∣ ≃ ℓ ∫ 0 ∞ s 2 e − s 4 I 1 ( h s ) d s ∫ 0 ∞ s e − s 4 I 0 ( h s ) d s , (190) |\langle b\rangle| \simeq\ell\frac{\int_0^\infty s^2e^{-s^4}I_1(hs)\,ds}{\int_0^\infty se^{-s^4}I_0(hs)\,ds}, \tag*{(190)} ∣ ⟨ b ⟩ ∣ ≃ ℓ ∫ 0 ∞ s e − s 4 I 0 ( h s ) d s ∫ 0 ∞ s 2 e − s 4 I 1 ( h s ) d s , (190) and
N b ≃ ℓ 2 ∫ 0 ∞ s 3 e − s 4 I 0 ( h s ) d s ∫ 0 ∞ s e − s 4 I 0 ( h s ) d s , (191) N_b \simeq\ell^2\frac{\int_0^\infty s^3e^{-s^4}I_0(hs)\,ds}{\int_0^\infty se^{-s^4}I_0(hs)\,ds}, \tag*{(191)} N b ≃ ℓ 2 ∫ 0 ∞ s e − s 4 I 0 ( h s ) d s ∫ 0 ∞ s 3 e − s 4 I 0 ( h s ) d s , (191) where I ν I_\nu I ν F F F ⟨ b ⟩ \langle b\rangle ⟨ b ⟩ F F F
At F = 0 F=0 F = 0 n = ∣ β ∣ 2 n=|\beta|^2 n = ∣ β ∣ 2
p b ( n ) ≃ 2 Γ 2 π K exp ( − Γ 2 n 2 K ) , n ≥ 0. (192) p_b(n) \simeq2\sqrt{\frac{\Gamma_2}{\pi K}}\exp\left(-\frac{\Gamma_2n^2}{K}\right),\qquad n\geq0. \tag*{(192)} p b ( n ) ≃ 2 π K Γ 2 exp ( − K Γ 2 n 2 ) , n ≥ 0. (192) Consequently,
N b ≃ K π Γ 2 . (193) N_b \simeq\sqrt{\frac{K}{\pi\Gamma_2}}. \tag*{(193)} N b ≃ π Γ 2 K . (193) Expanding Eq. (D51) to first order in F F F
⟨ b ⟩ ≃ 2 F K ⟨ x 2 ⟩ F = 0 = F K ⟨ ∣ β ∣ 2 ⟩ F = 0 = F π K Γ 2 = − 2 g α a π K Γ 2 . (194) \langle b\rangle\simeq\frac{2F}{K}\langle x^2\rangle_{F=0}=\frac{F}{K}\langle|\beta|^2\rangle_{F=0}
=\frac{F}{\sqrt{\pi K\Gamma_2}}=-\frac{2g\alpha_a}{\sqrt{\pi K\Gamma_2}}. \tag*{(194)} ⟨ b ⟩ ≃ K 2 F ⟨ x 2 ⟩ F = 0 = K F ⟨ ∣ β ∣ 2 ⟩ F = 0 = π K Γ 2 F = − π K Γ 2 2 g α a . (194) The downstream zero-drive susceptibility is defined by
χ b ≡ ∂ ⟨ b ⟩ ∂ Ω ∣ Ω = 0 . (195) \chi_b \equiv\left.\frac{\partial\langle b\rangle}{\partial\Omega}\right|_{\Omega=0}. \tag*{(195)} χ b ≡ ∂ Ω ∂ ⟨ b ⟩ Ω = 0 . (195) Using F = − 2 g α a F=-2g\alpha_a F = − 2 g α a
χ b ≃ − 2 g χ a π K Γ 2 . (196) \chi_b \simeq-\frac{2g\chi_a}{\sqrt{\pi K\Gamma_2}}. \tag*{(196)} χ b ≃ − π K Γ 2 2 g χ a . (196) The linear-response result requires
∣ F ∣ ≪ K 3 / 4 Γ 2 1 / 4 . (197) |F| \ll K^{3/4}\Gamma_2^{1/4}. \tag*{(197)} ∣ F ∣ ≪ K 3/4 Γ 2 1/4 . (197) Fixed nonzero drive
Now keep Ω ≠ 0 \Omega\ne0 Ω = 0 Γ 2 → 0 + \Gamma_2\to0^+ Γ 2 → 0 + α a \alpha_a α a Γ 2 \Gamma_2 Γ 2 F = − 2 g α a F=-2g\alpha_a F = − 2 g α a α a ≠ 0 \alpha_a\ne0 α a = 0
In this regime, the stationary distribution concentrates near the deterministic fixed point
F = 2 Γ 2 ∣ β 0 ∣ 2 β 0 . (198) F = 2\Gamma_2|\beta_0|^2\beta_0. \tag*{(198)} F = 2 Γ 2 ∣ β 0 ∣ 2 β 0 . (198) Thus,
β 0 = F ∣ F ∣ ( ∣ F ∣ 2 Γ 2 ) 1 / 3 . (199) \beta_0 = \frac{F}{|F|}\left(\frac{|F|}{2\Gamma_2}\right)^{1/3}. \tag*{(199)} β 0 = ∣ F ∣ F ( 2 Γ 2 ∣ F ∣ ) 1/3 . (199) For Ω > 0 \Omega>0 Ω > 0 α a > 0 \alpha_a>0 α a > 0 F = − 2 g α a < 0 F=-2g\alpha_a<0 F = − 2 g α a < 0
⟨ b ⟩ ≃ − ( g α a Γ 2 ) 1 / 3 , N b ≃ ( g α a Γ 2 ) 2 / 3 . (200) \langle b\rangle\simeq-\left(\frac{g\alpha_a}{\Gamma_2}\right)^{1/3},\qquad N_b\simeq\left(\frac{g\alpha_a}{\Gamma_2}\right)^{2/3}. \tag*{(200)} ⟨ b ⟩ ≃ − ( Γ 2 g α a ) 1/3 , N b ≃ ( Γ 2 g α a ) 2/3 . (200) More generally,
∣ ⟨ b ⟩ ∣ ∝ Γ 2 − 1 / 3 , N b ∝ Γ 2 − 2 / 3 . (201) |\langle b\rangle| \propto\Gamma_2^{-1/3},\qquad N_b\propto\Gamma_2^{-2/3}. \tag*{(201)} ∣ ⟨ b ⟩ ∣ ∝ Γ 2 − 1/3 , N b ∝ Γ 2 − 2/3 . (201) Homodyne measurement of the drive amplitude
We now consider the estimation of the real drive amplitude Ω \Omega Ω b b b
Homodyne quadrature and error propagation
For a local-oscillator phase θ \theta θ
X θ = b e − i θ + b † e i θ 2 . (202) X_{\theta}=\frac{b e^{-i\theta}+b^{\dagger}e^{i\theta}}{\sqrt{2}}. \tag*{(202)} X θ = 2 b e − i θ + b † e i θ . (202) For a stationary state depending on Ω \Omega Ω
Δ Ω θ ( 1 ) = Var ( X θ ) ∣ ∂ Ω ⟨ X θ ⟩ ∣ . (203) \Delta\Omega_{\theta}^{(1)}=\frac{\sqrt{\operatorname{Var}(X_{\theta})}}{\left|\partial_{\Omega}\langle X_{\theta}\rangle\right|}. \tag*{(203)} Δ Ω θ ( 1 ) = ∣ ∂ Ω ⟨ X θ ⟩ ∣ Var ( X θ ) . (203) The optimal homodyne angle is the one for which the measured quadrature is parallel to the response vector
∂ Ω ⟨ b ⟩ = ∣ ∂ Ω ⟨ b ⟩ ∣ e i θ o p t . \partial_{\Omega}\langle b\rangle=\left|\partial_{\Omega}\langle b\rangle\right|e^{i\theta_{\mathrm{opt}}}. ∂ Ω ⟨ b ⟩ = ∣ ∂ Ω ⟨ b ⟩ ∣ e i θ opt . Thus,
X o p t = b e − i θ o p t + b † e i θ o p t 2 . (204) X_{\mathrm{opt}}=\frac{b e^{-i\theta_{\mathrm{opt}}}+b^{\dagger}e^{i\theta_{\mathrm{opt}}}}{\sqrt{2}}. \tag*{(204)} X opt = 2 b e − i θ opt + b † e i θ opt . (204) For a real drive and a symmetry-preserving stationary state, both ⟨ b ⟩ \langle b\rangle ⟨ b ⟩ ∂ Ω ⟨ b ⟩ \partial_{\Omega}\langle b\rangle ∂ Ω ⟨ b ⟩
X = X 0 = b + b † 2 . (205) X=X_{0}=\frac{b+b^{\dagger}}{\sqrt{2}}. \tag*{(205)} X = X 0 = 2 b + b † . (205) In the effective classical description β = x + i y \beta=x+iy β = x + i y X ≃ 2 x X\simeq\sqrt{2}x X ≃ 2 x
Zero-drive linear-response sensitivity
At zero drive, the downstream stationary distribution is isotropic. Hence,
⟨ X ⟩ Ω = 0 = 0 , Var ( X ) Ω = 0 ≃ 2 ⟨ x 2 ⟩ 0 = N b . (206) \langle X\rangle_{\Omega=0}=0,\qquad\operatorname{Var}(X)_{\Omega=0}\simeq2\langle x^{2}\rangle_{0}=N_b. \tag*{(206)} ⟨ X ⟩ Ω = 0 = 0 , Var ( X ) Ω = 0 ≃ 2 ⟨ x 2 ⟩ 0 = N b . (206) Using N b ≃ K π Γ 2 N_b\simeq\sqrt{\frac{K}{\pi\Gamma_2}} N b ≃ π Γ 2 K ∂ ⟨ X ⟩ ∂ Ω ∣ Ω = 0 = 2 χ b \left.\frac{\partial\langle X\rangle}{\partial\Omega}\right|_{\Omega=0}=\sqrt{2}\chi_b ∂ Ω ∂ ⟨ X ⟩ Ω = 0 = 2 χ b
Δ Ω h o m ( 1 ) ≃ N b 2 ∣ χ b ∣ . (207) \Delta\Omega_{\mathrm{hom}}^{(1)}\simeq\frac{\sqrt{N_b}}{\sqrt{2}|\chi_b|}. \tag*{(207)} Δ Ω hom ( 1 ) ≃ 2 ∣ χ b ∣ N b . (207) Homodyne measurement of the drive amplitude
Using Eq. (196) ,
Δ Ω h o m ( 1 ) ≃ π 1 / 4 K 3 / 4 Γ 2 1 / 4 2 2 g ∣ χ a ∣ . (208) \Delta\Omega_{\mathrm{hom}}^{(1)}\simeq\frac{\pi^{1/4}K^{3/4}\Gamma_2^{1/4}}{2\sqrt{2}g|\chi_a|}. \tag*{(208)} Δ Ω hom ( 1 ) ≃ 2 2 g ∣ χ a ∣ π 1/4 K 3/4 Γ 2 1/4 . (208) Therefore,
Δ Ω h o m ( 1 ) ∝ Γ 2 1 / 4 . (209) \Delta\Omega_{\mathrm{hom}}^{(1)}\propto\Gamma_2^{1/4}. \tag*{(209)} Δ Ω hom ( 1 ) ∝ Γ 2 1/4 . (209) This result is valid only inside the zero-drive linear-response window
∣ F ∣ ≪ K 3 / 4 Γ 2 1 / 4 . (210) |F|\ll K^{3/4}\Gamma_2^{1/4}. \tag*{(210)} ∣ F ∣ ≪ K 3/4 Γ 2 1/4 . (210) For fixed nonzero Ω \Omega Ω Γ 2 → 0 + \Gamma_2\to0^{+} Γ 2 → 0 +
Fixed nonzero drive
We next consider fixed Ω ≠ 0 \Omega\ne0 Ω = 0 F ≠ 0 F\ne0 F = 0 Γ 2 \Gamma_2 Γ 2
r 0 = ∣ β 0 ∣ = ( ∣ F ∣ 2 Γ 2 ) 1 / 3 . (211) r_0=|\beta_0|=\left(\frac{|F|}{2\Gamma_2}\right)^{1/3}. \tag*{(211)} r 0 = ∣ β 0 ∣ = ( 2 Γ 2 ∣ F ∣ ) 1/3 . (211) The mean amplitude is aligned with the force. The optimal homodyne quadrature is therefore the quadrature parallel to F F F
The longitudinal fluctuation relaxation rate is
λ ∥ = 6 Γ 2 r 0 2 . (212) \lambda_{\parallel}=6\Gamma_2r_0^2. \tag*{(212)} λ ∥ = 6 Γ 2 r 0 2 . (212) Since the diffusion coefficient of the longitudinal real quadrature is K / 2 K/2 K /2
Var ( δ x ) ≃ K 12 Γ 2 r 0 2 . (213) \operatorname{Var}(\delta x)\simeq\frac{K}{12\Gamma_2r_0^2}. \tag*{(213)} Var ( δ x ) ≃ 12 Γ 2 r 0 2 K . (213) Therefore,
Var ( X ∥ ) ≃ 2 Var ( δ x ) = K 6 Γ 2 r 0 2 . (214) \operatorname{Var}(X_{\parallel})\simeq2\operatorname{Var}(\delta x)=\frac{K}{6\Gamma_2r_0^2}. \tag*{(214)} Var ( X ∥ ) ≃ 2 Var ( δ x ) = 6 Γ 2 r 0 2 K . (214) The mean quadrature is
⟨ X ∥ ⟩ ≃ 2 r 0 . (215) \langle X_{\parallel}\rangle\simeq\sqrt{2}r_0. \tag*{(215)} ⟨ X ∥ ⟩ ≃ 2 r 0 . (215) Since F F F Ω \Omega Ω
∂ r 0 ∂ Ω = r 0 3 F ∂ F ∂ Ω . (216) \frac{\partial r_0}{\partial\Omega}=\frac{r_0}{3F}\frac{\partial F}{\partial\Omega}. \tag*{(216)} ∂ Ω ∂ r 0 = 3 F r 0 ∂ Ω ∂ F . (216) with the appropriate real sign convention. In the weak-drive upstream approximation,
∂ F ∂ Ω = − 2 g χ a . (217) \frac{\partial F}{\partial\Omega}=-2g\chi_a. \tag*{(217)} ∂ Ω ∂ F = − 2 g χ a . (217) The single-sample homodyne uncertainty is therefore
Δ Ω ∥ ( 1 ) ≃ K / ( 6 Γ 2 r 0 2 ) 2 ∣ ∂ Ω r 0 ∣ = K 12 Γ 2 r 0 2 ∣ ∂ Ω r 0 ∣ . (218) \Delta\Omega_{\parallel}^{(1)} \simeq\frac{\sqrt{K/(6\Gamma_2 r_0^2)}}{\sqrt{2}\lvert\partial_{\Omega}r_0\rvert}
= \frac{\sqrt{K}}{\sqrt{12\Gamma_2 r_0^2}\lvert\partial_{\Omega}r_0\rvert}. \tag*{(218)} Δ Ω ∥ ( 1 ) ≃ 2 ∣ ∂ Ω r 0 ∣ K / ( 6 Γ 2 r 0 2 ) = 12 Γ 2 r 0 2 ∣ ∂ Ω r 0 ∣ K . (218) Using Eq. (216) ,
Δ Ω ∥ ( 1 ) ≃ ∣ F ∣ K 2 3 Γ 2 r 0 2 ∣ ∂ Ω F ∣ . (219) \Delta\Omega_{\parallel}^{(1)} \simeq\frac{\lvert F\rvert\sqrt{K}}{2\sqrt{3\Gamma_2}\,r_0^2\lvert\partial_{\Omega}F\rvert}. \tag*{(219)} Δ Ω ∥ ( 1 ) ≃ 2 3 Γ 2 r 0 2 ∣ ∂ Ω F ∣ ∣ F ∣ K . (219) Since r 0 2 ∝ Γ 2 − 2 / 3 r_0^2 \propto\Gamma_2^{-2/3} r 0 2 ∝ Γ 2 − 2/3
Δ Ω ∥ ( 1 ) ∝ Γ 2 1 / 6 . (220) \Delta\Omega_{\parallel}^{(1)} \propto\Gamma_2^{1/6}. \tag*{(220)} Δ Ω ∥ ( 1 ) ∝ Γ 2 1/6 . (220) The fixed-drive relaxation time scales as
τ ∥ ∼ ( Γ 2 r 0 2 ) − 1 ∼ Γ 2 − 1 / 3 ∣ F ∣ − 2 / 3 . (221) \tau_{\parallel} \sim\left(\Gamma_2r_0^2\right)^{-1} \sim\Gamma_2^{-1/3}\lvert F\rvert^{-2/3}. \tag*{(221)} τ ∥ ∼ ( Γ 2 r 0 2 ) − 1 ∼ Γ 2 − 1/3 ∣ F ∣ − 2/3 . (221) Consequently, in the long-measurement regime T ≫ τ ∥ ( Γ 2 ) T \gg\tau_{\parallel}(\Gamma_2) T ≫ τ ∥ ( Γ 2 )
Δ Ω ∥ ( T ) ≃ 2 τ ∥ T Δ Ω ∥ ( 1 ) ∝ Γ 2 0 . (222) \Delta\Omega_{\parallel}(T) \simeq\sqrt{\frac{2\tau_{\parallel}}{T}}\,\Delta\Omega_{\parallel}^{(1)} \propto\Gamma_2^0. \tag*{(222)} Δ Ω ∥ ( T ) ≃ T 2 τ ∥ Δ Ω ∥ ( 1 ) ∝ Γ 2 0 . (222) Again, the Γ 2 0 \Gamma_2^0 Γ 2 0 Γ 2 \Gamma_2 Γ 2 T T T
Homodyne estimation of Ω \Omega Ω G = 0 G=0 G = 0 Γ 2 → 0 \Gamma_2 \to0 Γ 2 → 0
We consider the estimation of the real driving amplitude Ω \Omega Ω b b b G = 0 G=0 G = 0 Γ 2 = 0 \Gamma_2=0 Γ 2 = 0
Model and assumptions
The Lindblad master equation at G = 0 G=0 G = 0
ρ ˙ = − i [ H , ρ ] + γ 1 D [ a ] ρ + γ 1 D [ a † ] ρ + γ 2 D [ a 2 ] ρ + κ b D [ b ] ρ + κ b D [ b † ] ρ + Γ 2 D [ b 2 ] ρ . (223) \dot{\rho}=-i[H,\rho]+\gamma_1\mathcal{D}[a]\rho+\gamma_1\mathcal{D}[a^\dagger]\rho+\gamma_2\mathcal{D}[a^2]\rho+\kappa_b\mathcal{D}[b]\rho+\kappa_b\mathcal{D}[b^\dagger]\rho+\Gamma_2\mathcal{D}[b^2]\rho. \tag*{(223)} ρ ˙ = − i [ H , ρ ] + γ 1 D [ a ] ρ + γ 1 D [ a † ] ρ + γ 2 D [ a 2 ] ρ + κ b D [ b ] ρ + κ b D [ b † ] ρ + Γ 2 D [ b 2 ] ρ . (223) with
H = i Ω ( a † − a ) + i g ( a † b − a b † ) , (224) H=i\Omega(a^\dagger-a)+ig(a^\dagger b-ab^\dagger), \tag*{(224)} H = i Ω ( a † − a ) + i g ( a † b − a b † ) , (224) where g ≠ 0 g\ne0 g = 0
D [ L ] ρ = 2 L ρ L † − L † L ρ − ρ L † L . (225) \mathcal{D}[L]\rho=2L\rho L^\dagger-L^\dagger L\rho-\rho L^\dagger L. \tag*{(225)} D [ L ] ρ = 2 L ρ L † − L † L ρ − ρ L † L . (225) Here γ 2 > 0 \gamma_2>0 γ 2 > 0 a a a Γ 2 = 0 \Gamma_2=0 Γ 2 = 0
Exact displacement relation at Γ 2 = 0 \Gamma_2=0 Γ 2 = 0
Introduce the displaced mode
c = b + Ω g , b = c − Ω g . (226) c=b+\frac{\Omega}{g},\qquad b=c-\frac{\Omega}{g}. \tag*{(226)} c = b + g Ω , b = c − g Ω . (226) Substituting Eq. (226) into Eq. (224) gives
H = i Ω ( a † − a ) + i g [ a † ( c − Ω g ) − a ( c † − Ω g ) ] = i g ( a † c − a c † ) . (227) H=i\Omega(a^\dagger-a)+ig\left[a^\dagger\left(c-\frac{\Omega}{g}\right)-a\left(c^\dagger-\frac{\Omega}{g}\right)\right]
=ig(a^\dagger c-ac^\dagger). \tag*{(227)} H = i Ω ( a † − a ) + i g [ a † ( c − g Ω ) − a ( c † − g Ω ) ] = i g ( a † c − a c † ) . (227) Thus, the drive disappears from the Hamiltonian in the displaced variables.
The balanced linear gain and loss of mode b b b s = Ω / g s=\Omega/g s = Ω/ g
D [ c − s ] ρ = D [ c ] ρ + s [ c † − c , ρ ] , (228) \mathcal{D}[c-s]\rho=\mathcal{D}[c]\rho+s[c^\dagger-c,\rho], \tag*{(228)} D [ c − s ] ρ = D [ c ] ρ + s [ c † − c , ρ ] , (228) D [ c † − s ] ρ = D [ c † ] ρ + s [ c − c † , ρ ] . (229) \mathcal{D}[c^\dagger-s]\rho=\mathcal{D}[c^\dagger]\rho+s[c-c^\dagger,\rho]. \tag*{(229)} D [ c † − s ] ρ = D [ c † ] ρ + s [ c − c † , ρ ] . (229) Consequently,
D [ b ] ρ + D [ b † ] ρ = D [ c ] ρ + D [ c † ] ρ . (230) \mathcal{D}[b]\rho+\mathcal{D}[b^\dagger]\rho=\mathcal{D}[c]\rho+\mathcal{D}[c^\dagger]\rho. \tag*{(230)} D [ b ] ρ + D [ b † ] ρ = D [ c ] ρ + D [ c † ] ρ . (230) At Γ 2 = 0 \Gamma_2=0 Γ 2 = 0 a a a c c c ρ 0 \rho_0 ρ 0
ρ s s ( Ω ) = D b ( − Ω g ) ρ 0 D b † ( − Ω g ) , (231) \rho_{\mathrm{ss}}(\Omega)=D_b\left(-\frac{\Omega}{g}\right)\rho_0D_b^\dagger\left(-\frac{\Omega}{g}\right), \tag*{(231)} ρ ss ( Ω ) = D b ( − g Ω ) ρ 0 D b † ( − g Ω ) , (231) where
D b ( α ) = exp ( α b † − α ∗ b ) . (232) D_b(\alpha)=\exp(\alpha b^\dagger-\alpha^*b). \tag*{(232)} D b ( α ) = exp ( α b † − α ∗ b ) . (232) Eq. (231) is exact and does not require a Gaussian approximation. It relies on the balanced linear gain and loss of mode b b b Γ 2 = 0 \Gamma_2=0 Γ 2 = 0
At zero drive, the master equation is invariant under
a ↦ a e i ϕ , b ↦ b e i ϕ . (233) a\mapsto ae^{i\phi},\qquad b\mapsto be^{i\phi}. \tag*{(233)} a ↦ a e i ϕ , b ↦ b e i ϕ . (233) For a phase-invariant steady state,
⟨ b ⟩ 0 = 0 , ⟨ b 2 ⟩ 0 = 0. (234) \langle b\rangle_0=0,\qquad\langle b^2\rangle_0=0. \tag*{(234)} ⟨ b ⟩ 0 = 0 , ⟨ b 2 ⟩ 0 = 0. (234) Define the zero-drive occupation
n b , 0 = ⟨ b † b ⟩ 0 . (235) n_{b,0}=\langle b^\dagger b\rangle_0. \tag*{(235)} n b , 0 = ⟨ b † b ⟩ 0 . (235) Using Eq. (231) , we obtain
⟨ b ⟩ s s = − Ω g , ∂ ⟨ b ⟩ s s ∂ Ω = − 1 g . (236) \langle b\rangle_{\mathrm{ss}}=-\frac{\Omega}{g},\qquad\frac{\partial\langle b\rangle_{\mathrm{ss}}}{\partial\Omega}=-\frac{1}{g}. \tag*{(236)} ⟨ b ⟩ ss = − g Ω , ∂ Ω ∂ ⟨ b ⟩ ss = − g 1 . (236) The second moments are
⟨ b † b ⟩ s s = n b , 0 + Ω 2 g 2 , ⟨ b 2 ⟩ s s = Ω 2 g 2 . (237) \langle b^\dagger b\rangle_{\mathrm{ss}}=n_{b,0}+\frac{\Omega^2}{g^2},\qquad\langle b^2\rangle_{\mathrm{ss}}=\frac{\Omega^2}{g^2}. \tag*{(237)} ⟨ b † b ⟩ ss = n b , 0 + g 2 Ω 2 , ⟨ b 2 ⟩ ss = g 2 Ω 2 . (237) The displacement contribution to the occupation must therefore be distinguished from the fluctuation noise. Consider the homodyne quadrature
X θ = b e − i θ + b † e i θ 2 . (238) X_\theta= \frac{b e^{-i\theta} + b^\dagger e^{i\theta}}{\sqrt{2}}. \tag*{(238)} X θ = 2 b e − i θ + b † e i θ . (238) Its mean and response are
⟨ X θ ⟩ s s = − 2 Ω g cos θ , ∂ ⟨ X θ ⟩ s s ∂ Ω = − 2 g cos θ . (239) \langle X_\theta\rangle_{\mathrm{ss}} = -\frac{\sqrt{2}\Omega}{g}\cos\theta,\qquad\frac{\partial\langle X_\theta\rangle_{\mathrm{ss}}}{\partial\Omega} = -\frac{\sqrt{2}}{g}\cos\theta. \tag*{(239)} ⟨ X θ ⟩ ss = − g 2 Ω cos θ , ∂ Ω ∂ ⟨ X θ ⟩ ss = − g 2 cos θ . (239) A displacement does not alter the centered fluctuations. Since the zero-drive steady state is phase invariant, it follows that
Var ( X θ ) = n b , 0 + 1 2 . (240) \operatorname{Var}(X_\theta) = n_{b,0} + \frac{1}{2}. \tag*{(240)} Var ( X θ ) = n b , 0 + 2 1 . (240) The variance is independent of both Ω \Omega Ω θ \theta θ
θ o p t = 0 ( m o d π ) . (241) \theta_{\mathrm{opt}} = 0 \pmod{\pi}. \tag*{(241)} θ opt = 0 ( mod π ) . (241) Accordingly, we measure
X = b + b † 2 . (242) X = \frac{b+b^\dagger}{\sqrt{2}}. \tag*{(242)} X = 2 b + b † . (242) According to the error propagation formula, the measurement uncertainty of Ω \Omega Ω
δ Ω = Var ( X ) ∣ ∂ ⟨ b ⟩ s s ∂ Ω ∣ = ∣ g ∣ 2 n b , 0 + 1 2 . (243) \delta\Omega= \frac{\sqrt{\operatorname{Var}(X)}}{\left|\frac{\partial\langle b\rangle_{\mathrm{ss}}}{\partial\Omega}\right|} = \frac{|g|}{\sqrt{2}}\sqrt{n_{b,0}+\frac{1}{2}}. \tag*{(243)} δ Ω = ∂ Ω ∂ ⟨ b ⟩ ss Var ( X ) = 2 ∣ g ∣ n b , 0 + 2 1 . (243) This result does not require Gaussian-distributed measurement outcomes. To express n b , 0 n_{b,0} n b , 0 Γ 2 = 0 \Gamma_2=0 Γ 2 = 0
n a = ⟨ a † a ⟩ 0 , n b = ⟨ b † b ⟩ 0 = n b , 0 , C = ⟨ a † b ⟩ 0 . (244) n_a=\langle a^\dagger a\rangle_0,\qquad n_b=\langle b^\dagger b\rangle_0=n_{b,0},\qquad C=\langle a^\dagger b\rangle_0. \tag*{(244)} n a = ⟨ a † a ⟩ 0 , n b = ⟨ b † b ⟩ 0 = n b , 0 , C = ⟨ a † b ⟩ 0 . (244) The master equation gives
n ˙ a = 2 γ 1 + − 4 γ 2 ⟨ a † 2 a 2 ⟩ + 2 g Re C , (245) \dot n_a=2\gamma_1^+-4\gamma_2\langle a^{\dagger2}a^2\rangle+2g\operatorname{Re}C, \tag*{(245)} n ˙ a = 2 γ 1 + − 4 γ 2 ⟨ a † 2 a 2 ⟩ + 2 g Re C , (245) n ˙ b = 2 Γ 1 + − 2 g Re C , (246) \dot n_b=2\Gamma_1^+-2g\operatorname{Re}C, \tag*{(246)} n ˙ b = 2 Γ 1 + − 2 g Re C , (246) C ˙ = g ( n b − n a ) − 2 γ 2 ⟨ a † 2 a b ⟩ . (247) \dot C=g(n_b-n_a)-2\gamma_2\langle a^{\dagger2}ab\rangle. \tag*{(247)} C ˙ = g ( n b − n a ) − 2 γ 2 ⟨ a † 2 ab ⟩ . (247) At steady state, Eqs. (E23) and (E24) imply
Re C = Γ 1 + g , ⟨ a † 2 a 2 ⟩ 0 = γ 1 + + Γ 1 + 2 γ 2 . (248) \operatorname{Re}C=\frac{\Gamma_1^+}{g},\qquad\langle a^{\dagger2}a^2\rangle_0=\frac{\gamma_1^++\Gamma_1^+}{2\gamma_2}. \tag*{(248)} Re C = g Γ 1 + , ⟨ a † 2 a 2 ⟩ 0 = 2 γ 2 γ 1 + + Γ 1 + . (248) However, Eq. (E25) contains a fourth-order correlation. The second-order moment equations are therefore not closed, and Eq. (E26) alone does not determine n b , 0 n_{b,0} n b , 0
Gaussian moment-closure approximation
To obtain an explicit parameter-dependent expression, we adopt a zero-mean, phase-invariant Gaussian moment closure:
⟨ a † 2 a 2 ⟩ ≃ 2 n a 2 , ⟨ a † 2 a b ⟩ ≃ 2 n a C . (249) \langle a^{\dagger2}a^2\rangle\simeq2n_a^2,\qquad\langle a^{\dagger2}ab\rangle\simeq2n_aC. \tag*{(249)} ⟨ a † 2 a 2 ⟩ ≃ 2 n a 2 , ⟨ a † 2 ab ⟩ ≃ 2 n a C . (249) The first relation and Eq. (E26) give
n a G a u s s = γ 1 + + Γ 1 + 4 γ 2 . (250) n_a^{\mathrm{Gauss}}=\sqrt{\frac{\gamma_1^++\Gamma_1^+}{4\gamma_2}}. \tag*{(250)} n a Gauss = 4 γ 2 γ 1 + + Γ 1 + . (250) The steady-state equation for C C C
g ( n b − n a ) − 4 γ 2 n a C = 0. (251) g(n_b-n_a)-4\gamma_2n_aC=0. \tag*{(251)} g ( n b − n a ) − 4 γ 2 n a C = 0. (251) Within this closure, C C C
C G a u s s = Γ 1 + g . (252) C^{\mathrm{Gauss}}=\frac{\Gamma_1^+}{g}. \tag*{(252)} C Gauss = g Γ 1 + . (252) Consequently,
n b , 0 G a u s s = ( 1 + 4 Γ 1 + γ 2 g 2 ) γ 1 + + Γ 1 + 4 γ 2 . (253) n_{b,0}^{\mathrm{Gauss}}=\left(1+\frac{4\Gamma_1^+\gamma_2}{g^2}\right)\sqrt{\frac{\gamma_1^++\Gamma_1^+}{4\gamma_2}}. \tag*{(253)} n b , 0 Gauss = ( 1 + g 2 4 Γ 1 + γ 2 ) 4 γ 2 γ 1 + + Γ 1 + . (253) Substituting Eq. (E31) into Eq. (E21) yields
δ Ω G a u s s = 1 2 [ ( g 2 + 4 Γ 1 + γ 2 ) γ 1 + + Γ 1 + 4 γ 2 + g 2 2 ] 1 / 2 . (254) \delta\Omega_{\mathrm{Gauss}}=\frac{1}{\sqrt{2}}\left[\left(g^2+4\Gamma_1^+\gamma_2\right)\sqrt{\frac{\gamma_1^++\Gamma_1^+}{4\gamma_2}}+\frac{g^2}{2}\right]^{1/2}. \tag*{(254)} δ Ω Gauss = 2 1 ( g 2 + 4 Γ 1 + γ 2 ) 4 γ 2 γ 1 + + Γ 1 + + 2 g 2 1/2 . (254) For equal linear rates, γ 1 + = Γ 1 + = κ \gamma_1^+=\Gamma_1^+=\kappa γ 1 + = Γ 1 + = κ
δ Ω G a u s s = 1 2 [ ( g 2 + 4 κ γ 2 ) κ 2 γ 2 + g 2 2 η ] 1 / 2 . (255) \delta\Omega_{\mathrm{Gauss}}=\frac{1}{\sqrt{2}}\left[\left(g^2+4\kappa\gamma_2\right)\sqrt{\frac{\kappa}{2\gamma_2}}+\frac{g^2}{2\eta}\right]^{1/2}. \tag*{(255)} δ Ω Gauss = 2 1 [ ( g 2 + 4 κ γ 2 ) 2 γ 2 κ + 2 η g 2 ] 1/2 . (255) The limit g → 0 g\to0 g → 0
