A quantum ticking self-oscillator using delayed feedback

First studied in 1994 [20], time-delay in quantum feedback has not enjoyed much progress primarily due to the significant challenge faced to obtain a valid master equation. When the delay time is not infinitesimal, the quantum dynamics becomes non-Markovian [21], which can be quite problematic. There are, to date only a few methods discovered to describe a quantum master equation which incorporates time delayed feedback. These methods involve the use of cascaded quantum channels and can use matrix product states (MPS) to help cope with the very large Hilbert spaces that appear [21, 22]. By considering the time-delayed quantum system in a cascaded way, where the current system is driven (affected) in a temporally directed fashion by its past dynamics, a master equation in piece-wise form can be derived [23]. In this section, we utilise and generalise the cascaded theory to get the master equation for quantum self-oscillator with time delay.

We can first write the Langevin equation of the cavity mode $\hat{a}$ as:

$$\eqalign{\dot{\hat{a}}(t)&=-\kappa\hat{a}(t)-\sqrt{\kappa_{1}}\nu_{1}(t)-\sqrt{\kappa_{2}}\nu_{2}(t),\cr\tilde{\nu}_{1}(t)&=\sqrt{\kappa_{1}}\hat{a}(t)+\nu_{1}(t),\\ \tilde{\nu}_{2}(t)&=\sqrt{\kappa_{2}}\hat{a}(t)+\nu_{2}(t),}$$

(1)

where $\nu_{1}(t)$ ($\nu_{2}(t)$) are the input fields of the left (right) partial mirrors with $\tilde{\nu}_{1}(t)$ ($\tilde{\nu}_{2}(t)$) being the output fields of left (right) partial mirrors. The decay rates of these mirrors are $\kappa_{1}$ ($\kappa_{2}$) respectively, with $\kappa=\frac{\kappa_{1}+\kappa_{2}}{2}$ being the average decay of the cavity.

We insert a quantum amplifier into the feedback loop, so we get $G>1$. It is well known that such amplification cannot be achieved without incorporating noise and continuous-wave driving [24]. We can express the returning feedback signal $\nu_{2}(t)$ following a delay, and amplification on the signal $\tilde{\nu}_{1}(t-\tau)$, as:

$$\nu_{2}(t)=e^{i\phi}\left(\sqrt{G}\,\tilde{\nu}_{1}(t-\tau)+\sqrt{G-1}\,\nu_{amp}^{\dagger}(t-\tau)\right),$$

(2)

where $\phi$ is the phase change in the feedback loop while $\nu_{amp}$ is the amplifier noise field [24]. Substituting Equation (2) into (1), we obtain

$$\eqalign{\dot{\hat{a}}(t)&=-\kappa\hat{a}(t)-e^{i\phi}\sqrt{G\kappa_{1}\kappa_{2}}\,\hat{a}(t-\tau)-\sqrt{\kappa_{1}}\,\nu_{1}(t)-\sqrt{G\kappa_{2}}\,\nu_{1}(t-\tau)\\ &-e^{i\phi}\sqrt{(G-1)\kappa_{2}}\,\nu_{amp}^{\dagger}(t-\tau).}$$

(3)

The noise field of the amplifier appears in the form of $\nu_{amp}^{\dagger}$ in (3), which will inject thermal noise into the system due to the relations $d\nu_{amp}^{\dagger}(t)d\nu_{amp}(t)=\bar{N}_{amp}dt$ and $d\nu_{amp}(t)d\nu_{amp}^{\dagger}(t)=(\bar{N}_{amp}+1)dt$ (see A for details). If we assume the input field of the left mirror $\nu_{1}(t)$ and the quantum amplifier noise field $\nu_{amp}(t)$ are all vacuum fields, taking the expectation values of the above quantum DDE (3) leads to:

$$\langle\dot{\hat{a}}(t)\rangle=-\kappa\langle\hat{a}(t)\rangle-e^{i\phi}\sqrt{G\kappa_{1}\kappa_{2}}\,\langle\hat{a}(t-\tau)\rangle.$$

(4)

This quantum DDE has exactly the same form to the following general classical DDE [18]

$$\dot{x}(t)=\alpha x(t)+\beta x(t-\tau),$$

(5)

where $\alpha\in\mathbb{R},\beta\in\mathbb{R}$ are the usual system parameters.

This classical DDE has been extensively explored [25] and we will investigate its properties below. In particular, we will obtain conditions on the system parameters $\alpha,\beta$, and $\tau$ to obtain a self-sustained oscillation of $x(t)$, and hence $\langle\hat{a}(t)\rangle$ in Equation (4).

Systems with time-delay feedback are usually described using DDEs in the form of Equation (5). Such DDEs have been widely studied in various fields such as population biology [26], physiology [27], epidemiology [28], economics [25], neural networks [29], and electronic circuits [30] and thus provide an ideal basis to explore SSOs.

If we assume the solution of the Equation (5) has the exponential form $x(t)=e^{\lambda t}$, then substituting this solution into Equation (5) gives us the characteristic equation (CE), $\lambda-\alpha-\beta e^{-\lambda\tau}=0$ which we can rewrite as

$$z-\alpha^{\prime}-\beta^{\prime}e^{-z}=0,$$

(6)

where we have multiplied the CE by $\tau$, and denote $z=\lambda\tau,\alpha^{\prime}=\alpha\tau,\beta^{\prime}=\beta\tau$. The solutions of Equation (6) for the characteristic values of $z$ can be written as $z=\Re{z}+i\Im{z}$, which means that the solution to Equation (5) is $x(t)=\textrm{Exp}\left[\frac{\Re{z}}{\tau}t\right]\textrm{Exp}\left[i\frac{\Im{z}}{\tau}t\right]$. The behavior $\textrm{Exp}\left[\frac{\Re{z}}{\tau}t\right]$ describes whether $x(t)$ is increasing or decaying based on the sign of $\Re{z}$, while latter $\textrm{Exp}\left[i\frac{\Im{z}}{\tau}t\right]$, describes oscillations of $x(t)$. With this notation, Equation (6) becomes:

$$\eqalign{\Re{z}-\alpha^{\prime}-\beta^{\prime}e^{-\Re{z}}\cos{\Im{z}}&=0,\\ \Im{z}+\beta^{\prime}e^{-\Re{z}}\sin{\Im{z}}&=0.}$$

(7)

To ensure $x(t)$ oscillates, we need to guarantee that Equation (6) has purely imaginary solutions, that is, $\Re{z}=0$. The values of $(\alpha^{\prime},\beta^{\prime})$ satisfying this fall into separate classes depending on the argument of $z$, with the parameters $(\alpha^{\prime},\beta^{\prime})$ in sets of one-dimensional curves, which we denote as $C_{j}$ with

$$\left\{\eqalign{(\alpha^{\prime},\beta^{\prime})\in C_{0}&,\Im{z}\in(0,\pi)\\ (\alpha^{\prime},\beta^{\prime})\in C_{1}&,\Im{z}\in(\pi,2\pi)\\ &\vdots\\ (\alpha^{\prime},\beta^{\prime})\in C_{j}&,\Im{z}\in(j\pi,(j+1)\pi)}\right.$$

as shown in Figure 2. We consider $(\alpha^{\prime},\beta^{\prime})\in C_{0}$ from here on and these solutions display the greatest level of robustness. In the right of Figure 2 we plot the SSO of Equation (5) for $(\alpha^{\prime},\beta^{\prime})=(-3.573,-4.375)$, where we consider two initial states, $x(t=0)=1$ and $x(t=0)=0.5$ with $x(t<0)=0$. We observe that by choosing parameters on the line $C_{0}$ we can obtain perfect self-oscillation for all time where the oscillation amplitude depends on the initial value.

From studying the classical DDEs (5), we find a one-dimensional set of parameters $(\alpha^{\prime},\beta^{\prime})$ that displays self-sustained oscillation in $x(t)$. When $(\alpha^{\prime},\beta^{\prime})\in C_{0}$, we can always find a critical delay time $\tau_{cr}=\arccos(-\frac{\alpha}{\beta})/\sqrt{\beta^{2}-\alpha^{2}}$ to achieve self-oscillation if $|\beta|>|\alpha|$. Returning now to the quantum expectation in Equation (4) for $\langle\hat{a}(t)\rangle$, we observe that in the case $\phi=0$, for any set of quantum parameters $(\kappa_{1},\kappa_{2},G)$, we can find a critical delay $\tau_{cr}$ such that the quantum system has pure self-oscillation in $\langle\hat{a}(t)\rangle$. The field will oscillate forever, with the only the CW drive required by the quantum amplifier $G$, as long as we can guarantee the delay time is set as $\tau_{cr}$.

Although the time-delayed quantum Langevin equations give an elegant description of the system, it is unclear how to solve them numerically. An alternative method is to derive a time-delayed master equation which can be solved numerically. One approach to derive such a master equation with time-delayed feedback is based on cascaded channel theory [31, 32]. This has been previously been applied to model quantum systems with time-delayed feedback when amplification is absent [23]. The basic idea is to copy the quantum system into multiple identical new systems, each of them representing the evolution of the original system only within a time interval $\tau$ - the delay time [21]. This we schematically depict in Figure 1-(b), where $S_{1},S_{2},\cdots,S_{m}$ are $m$ copies of the system. For $t\in[0,\tau]$, only the first copy $S_{1}$ evolves, which represents the evolution of the original system before the delayed feedback signal appears. When $t\in[\tau,2\tau]$, the second copy $S_{2}$ is involved to represent the evolution of current system, and the delayed feedback signal will be represented as the output of the previous copy $S_{1}$. Similarly, we can use a new copy $S_{i}$ to represent the most current system, which will be driven by its previous copy. By writing the master equation of this cascaded chain, we can derive the master equation of the time-delayed quantum system. However, the cascaded theory in [19] needs to be extended before we can apply it to our case. The reason is that we consider an amplifier in the feedback loop and this quantum amplifier noise that effects subsequent temporal dynamics will be amplified step by step.

To obtain such a master equation, we first start from the coupling Langevin equation. Then by defining all the input fields for each system in the cascaded chain, we can derive the Ito white noise quantum stochastic differential equation as:

$$\eqalign{d\hat{o}_{m}=-\frac{i}{\hbar}[\hat{o}_{m},H_{\textrm{tot}}]dt-\sum_{j=1}^{m}\frac{\gamma_{j}}{2}(\bar{N}_{j}+1)\left(2\hat{a}_{j}^{\dagger}\hat{o}_{m}\hat{a}_{j}-\hat{o}_{m}\hat{a}_{j}^{\dagger}\hat{a}_{j}-\hat{a}_{j}^{\dagger}\hat{a}_{j}\hat{o}_{m}\right)dt\\ -\sum_{j=1}^{m}\frac{\gamma_{j}}{2}\bar{N}_{j}\left(2\hat{a}_{j}\hat{o}_{m}\hat{a}_{j}^{\dagger}-\hat{o}_{m}\hat{a}_{j}\hat{c}_{j}^{\dagger}-\hat{a}_{j}\hat{a}_{j}^{\dagger}\hat{o}_{m}\right)dt\\ +\sum_{j=1}^{m}\left\{-[\hat{o}_{m},\hat{a}_{j}^{\dagger}]\sqrt{G\gamma_{j-1}\gamma_{n}}\hat{a}_{j-1}+\sqrt{G\gamma_{j-1}\gamma_{j}}\hat{a}_{j-1}^{\dagger}[\hat{o}_{m},\hat{a}_{j}]\right\}dt\\ +\sum_{j=1}^{m}\{-[\hat{o}_{m},\hat{a}_{j}^{\dagger}]\sqrt{\gamma_{j}}\left(\epsilon_{in}(j,t)dt+dB(j,t)\right)\cr+\sqrt{\gamma_{j}}\left(\epsilon^{*}_{in}(j,t)dt+dB^{\dagger}(j,t)\right)[\hat{o}_{m},\hat{a}_{j}]\}.}$$

(8)

Here we assume there are $m$ systems in the cascaded chain, and $H_{\textrm{tot}}$ is the sum over all system Hamiltonians. The input fields into the $j$-th system in the chain is defined as $db_{j}^{in}(t)=\epsilon_{in}(j,t)dt+dB(j,t)$, where $\epsilon_{in}(j,t)$ is the coherent part and $dB(j,t)$ is the white noise part. We use $\gamma_{j}$ as the decay rate of the $j$-th system, and all the operators have the form of:

$$\mathcal{O}_{j}=\underbrace{I\otimes\cdots\otimes I}_{j\textit{ times}}\otimes\mathcal{O}\otimes\underbrace{I\otimes\cdots\otimes I}_{(m-j-1)\textit{ times}}.$$

(9)

The operator $\mathcal{O}$ is in the Hilbert space of each individual system, and the operator $\mathcal{O}_{j}$ resides in the whole Hilbert space which is a tensor product of all the temporally separate systems. The details to obtain Equation (8) can be found in the B. Using this Ito equation, a master equation can be derived for the density operator $\rho_{m}(t)$ for the chained system, that is, $\rho_{m}(t)\in\otimes_{j=1}^{m}\mathcal{H}_{S}(j)$, where $\mathcal{H}_{S}(j)$ is the Hilbert space of individual systems in the chain. The master equation can be obtained by setting $\langle d\hat{o}(t)\rho_{m}\rangle=\langle d\hat{o}\rho_{m}(t)\rangle$, to obtain the evolution of the entire historical chain up to the $m$-th time interval $t\in[(m-1)\tau,m\tau]$:

$$\eqalign{\frac{d\rho_{m}}{dt}&=\frac{i}{\hbar}[\rho,H_{\textrm{tot}}]+\sum_{j=1}^{m}\left\{\gamma_{j}(\bar{N}_{j}+1)\mathcal{D}[\hat{a}_{j}]\rho+\gamma_{j}\bar{N}_{j}\mathcal{D}[\hat{a}_{j}^{\dagger}]\rho\right\}\\ &+\sum_{j=2}^{m}\sqrt{G\gamma_{j-1}\gamma_{j}}\left\{[\hat{a}_{j},\hat{a}_{j-1}\rho]+[\rho\hat{a}_{j-1}^{\dagger},\hat{a}_{j}]\right\}\\ &-\sum_{j=1}^{m}\left\{\sqrt{\gamma_{j}}\left[\hat{a}_{j}^{\dagger},\rho\right]\epsilon_{in}(j,t)-\sqrt{\gamma_{j}}\epsilon^{*}_{in}(j,t)\left[\hat{a}_{j},\rho\right]\right\}\\ }$$

(10)

with the initial state at $t=0$ being:

$$\rho_{m}(0)=\underbrace{\rho_{S}(0)\otimes\cdots\otimes\rho_{S}(0)}_{m\textit{ times}}.$$

(11)

This master equation has exactly the same form to the general cascaded theory in [19], but with the modified mean ancestor bath occupation of the fields $\bar{N}_{j}$, the detail of which can be found in B. If we remove the amplifier, by setting $G=1$, and let the input and amplifier noise fields be vacuum fields, Equation (10) reduces exactly to that derived in [33].

Armed with the general cascaded description incorporating amplification, we can now return to our original time-delayed system shown in Fig 1-(a). Similar to the idea in Figure 1-(b), we consider the time-delayed optical system in a cascaded way. This can be shown in Figure 1-(c), where we use one copy of the ancestor cavity to represent its evolution in different time intervals of length $\tau$. The output of one copy will be amplified and used to drive the next copy in the chain. To numerically simulate this we will use $\hat{a}_{0}$ to represent the current system, and $\hat{a}_{m}$ to represent older systems, where larger $m$ means older iterates. The master equation of the time-delayed optical system shown in Figure 1-(a) when $t\in[m\tau,(m+1)\tau]$ is thus:

$$\eqalign{\frac{d\rho}{dt}&=\mathcal{L}_{m}\rho,t\in[m\tau,(m+1)\tau],\textrm{with}\\ \mathcal{L}_{m}\rho&=-\frac{i}{\hbar}\sum_{j=0}^{m}[H_{j},\rho]+\sum_{j=0}^{m}\bar{N}_{m-j}\kappa_{1}\left(\mathcal{D}[\hat{a}_{j}]\rho+\mathcal{D}[\hat{a}_{j}^{\dagger}]\rho\right)\\ &+(\kappa_{1}+\kappa_{2})\sum_{j=0}^{m}\mathcal{D}[\hat{a}_{j}]\rho-\sqrt{G\kappa_{1}\kappa_{2}}\sum_{j=1}^{m}\left\{[\hat{a}_{j-1}^{\dagger},\hat{a}_{j}\rho]+[\rho\hat{a}_{j}^{\dagger},\hat{a}_{j-1}]\right\},}$$

(12)

where

$$\bar{N}_{m-j}=G\bar{N}+(G-1)(\bar{N}_{amp}+1),\textrm{for }j=0,1,\cdots,m,$$

(13)

where $\bar{N}$ is the bath occupation of the input field to the initial system, and $\bar{N}_{amp}$ is the bath occupation of the noise input to the amplifier. If we assume the input field at the initial time interval and the noise input to the amplifier are both vacuum field, this will simplify $\bar{N}_{m-j}$ as:

$$\bar{N}_{m-j}=G-1,\textrm{for }j=0,1,\cdots,m.$$

(14)

Therefore, we can identify $\hat{a}_{j}$ in Equation (12) with the operator $\hat{a}(t-j\tau)$ in the time-delay Langevin equation (4), and the expectation value of the most current time step we aim to calculate will correspond to the zeroth system $\langle\hat{a}_{0}\rangle$.

In this section, we will employ two approaches to simulate the evolution of the expectation values for the current system annihilation operator, $\langle\hat{a}_{0}\rangle$ and occupation $\langle\hat{a}_{0}^{\dagger}\hat{a}_{0}\rangle$. We observe that the $\langle\hat{a}_{0}\rangle$ exhibits indefinite oscillation without decay, albeit at the expense of an increasing energy.

The first approach uses the piece-wise master equation (12) to obtain $\langle\hat{a}_{0}\rangle={\rm Tr}(\hat{a}_{0}\rho)$. If we truncate the Hilbert space for each system in the chain to be $N_{trunc}$, the total dimension of the Hilbert space of the complete cascaded chain, when one evolves over $m$ delay intervals, grows as $N_{trunc}^{m}$. The computational cost to follow this grows to be prohibitive even when $m\sim 5-10$. This is the principle computational bottleneck for this approach. Such large Hilbert spaces motivates some researchers to explore using MPS to address this [22]. The increasing Hilbert space dimension required becomes quickly very expensive as regards computer storage, and typically we can only simulate for a few time steps. Furthermore, increasing the truncation $N_{trunc}$ will further increase the storage resources. Therefore, we propose a second approach to simulate $\langle\hat{a}_{0}\rangle$ for extended durations, which we call the mean-field approach. For a given time interval $t\in[m\tau,(m+1)\tau]$, we first derive a set of coupled ordinary differential equations (ODEs) in terms of $d\langle\hat{a}_{j}\rangle/dt,j=0,\cdots m$ in each time step from $[0,\tau]$ to $[m\tau,(m+1)\tau)]$. By solving this set of ODEs, we obtain the evolution of the most current system annihilation operator $\langle\hat{a}_{0}(t)\rangle$. We also wish to observe the evolution of photon number $\langle\hat{a}_{0}^{\dagger}\hat{a}_{0}\rangle$, with time. For this we will need to obtain sets of ODEs involving product moments between two systems. Although the number of required ODEs, $N_{Eqns}$, grows in time like $\sim O(m^{2})$, when integrating over $m$-delay periods, the required storage is still much smaller than that needed to solve the full quantum master equation. We perform several simulations using both the full quantum master equation and the above mentioned mean-field approach.

For simplicity, we can choose the decay rates $\kappa_{1}=\kappa_{2}=\kappa$, the feedback gain as $G=1.5$, which gives us $\alpha/\kappa=-1,\beta/\kappa=-1.225$, and $\kappa\tau\approx 3.592$. This gives the self-oscillation frequency of $\langle\hat{a}_{0}\rangle$ as $\omega_{\tau}/\kappa\approx 0.707$. Due to the gain $G>1$ in the feedback loop, we need to choose a suitable truncation of the Fock space and we explore the convergence of our simulations as a function of this truncation in Figure 3-(a).

From Figure 3-(a) we observe that $\langle\hat{a}_{0}\rangle$ shows oscillation in time due to the feedback gain starting from the second time interval. We also observe good convergence if we choose a Hilbert space truncation $N_{trunc}\geq 12$. Choose $N_{trunc}=12$ we evolve the photon number $\langle n\rangle=\langle\hat{a}_{0}^{\dagger}\hat{a}_{0}\rangle$, as well as the variance of the amplitude quadrature $X$ and phase quadrature $P$, calculated as $\Delta X=\sqrt{\langle X^{2}\rangle-\langle X\rangle^{2}}$ and $\Delta P=\sqrt{\langle P^{2}\rangle-\langle P\rangle^{2}}$. The results are shown in Figure 3-(b), which indicates that the photon number and the quadrature variance increase with time. We have no reason to expect that this behaviour will saturate in time, i.e. we expect the energy of the system to grow indefinitely.

In Figure 3, $\langle a\rangle$ exhibits sustained oscillation without any decay, which strongly suggests the absence of phase diffusion in phase space. For the delayed oscillator described by the linear DDE in Equation (5), the trajectories in phase space $(\langle x(t)\rangle,\langle p(t)\rangle)$, actually describe straight lines which is not typically what one considers to be a closed cycle. By modifying the system slightly we can devise periodic dynamics which do exhibit closed curves in phase space more clearly. To do this we include a non-zero frequency term in the Hamiltonian of the non-delayed system, which results in the following modified DDE:

$$\dot{x}(t)=i\omega x(t)+\alpha x(t)+\beta x(t-\tau)\;\;.$$

(15)

We can actually find new coordinates to relate (15), to the original DDE (5). To do this we multiply both sides of (15), by $e^{-i\omega t}$, and let $\tilde{x}(t)=e^{-i\omega t}x(t)$, and then we can rewrite the DDE for $\tilde{x}$ as:

$$\dot{\tilde{x}}(t)=\alpha\tilde{x}(t)+\beta e^{-i\omega\tau}\tilde{x}(t-\tau),$$

(16)

with complex coefficient $\beta^{\prime}=\beta e^{-i\omega\tau}$. By choosing $\omega\tau=2n\pi$ for an integer $n$, we can have $\beta^{\prime}=\beta$. In this scenario a closed curve cycle, (rather than a straight line), can be observed between the quadratures $x(t)$ and $p(t)$ when starting from a complex initial state.

We present simulation results for these closed curved cycles in Figure 4, for the mean quadratures either by solving their corresponding DDEs Equation (15), or by a full quantum simulation of the delayed dynamics.

In Figure 4-(a), we solve Equation (15) for two different different closed cycles corresponding to two values of $\beta$,for the same initial point. While in Figure 4-(b) we plot the solution of the DDE and the mean values of the quadratures $(\langle x(t)\rangle,\langle p\rangle)$, for the full delayed master equation. For the DDE we evolve forward in time from $t=0$, but we also have to specify the values of the dependent variables for $t<0$, and by exponentially damping them to zero $x,p\sim\exp(\gamma t),\;\;t<0$, we find a very good fit for the full quantum evolution over the first two delay time steps. The trajectories in phase space converge to their corresponding closed cycles. Notably, for the linear DDE, the closed cycles are influenced by decay rates, feedback gain, and initial states.

To explore the closed cycles in linear quantum systems, we include an extra Hamiltonian $-\omega\hbar\hat{a}^{\dagger}\hat{a}$ into the quantum system, which changes the original quantum system described in Equation (4) to

$$\langle\dot{\hat{a}}(t)\rangle=i\omega\langle\hat{a}(t)\rangle-\kappa\langle\hat{a}(t)\rangle-e^{i\phi}\sqrt{G\kappa_{1}\kappa_{2}}\,\langle\hat{a}(t-\tau)\rangle.$$

(17)

which with the scaled time $\kappa t$ has the form

$$\langle\dot{\hat{a}}(t)\rangle=(i\omega/\kappa-1)\langle\hat{a}(t)\rangle-e^{i\phi}\sqrt{G\kappa_{1}\kappa_{2}}/\kappa\,\langle\hat{a}(t-\tau)\rangle.$$

We then simulate our quantum systems using the same set of parameters $\kappa_{1}=\kappa_{2}=\kappa,G=4$ and $\omega\tau=2\pi$ as in the classical DDE. We expect to observe the same closed cycles in the quantum case. However, due to the extreme high dimensionality of full quantum cascaded master equation, we are unable to observe the final closed cycle on conventional computers. The excellent agreement between the quantum simulation and the classical DDE depicted in Figure 4-(b) is sufficient evidence to assert that we would observe closed cycles in quantum systems if we could simulate for a greater number of time steps.

It would be advantageous if one can obtain a similar nonlinear quantum self-oscillator by introducing a nonlinear damping in the quantum optical system. One of the options to achieve this is to add a two-photon absorption (TPA) [35] in our ring cavity. TPA is a nonlinear optical phenomenon that occurs when two photons are simultaneously absorbed. This phenomenon has been detected experimentally in many optical materials such as europium-doped crystal [36], cesium vapor [37], and Cds [38]. The TPA used in our setup is to absorb two photons from the cavity $\hat{a}$ at the same time, which will change the master equation of the optical system as:

$$\eqalign{\frac{d\rho}{dt}&=\mathcal{L}_{m}\rho,t\in[m\tau,(m+1)\tau],\textrm{where}\\ \mathcal{L}_{m}\rho&=-\frac{i}{\hbar}\sum_{j=0}^{m}[H_{j},\rho]+\sum_{j=0}^{m}\bar{N}_{m-j}\kappa_{1}\left(\mathcal{D}[\hat{a}_{j}]\rho+\mathcal{D}[\hat{a}_{j}^{\dagger}]\rho\right)+\gamma_{non}\sum_{j=0}^{m}\mathcal{D}[\hat{a}_{j}^{2}]\\ &+(\kappa_{1}+\kappa_{2})\sum_{j=0}^{m}\mathcal{D}[\hat{a}_{j}]\rho-\sqrt{G\kappa_{1}\kappa_{2}}\sum_{j=1}^{m}\left\{[\hat{a}_{j-1}^{\dagger},\hat{a}_{j}\rho]+[\rho\hat{a}_{j}^{\dagger},\hat{a}_{j-1}]\right\}\\ }$$

(19)

The DDE of $\langle\hat{a}\rangle$ of this nonlinear system will be:

$$\langle\dot{\hat{a}}(t)\rangle=-\kappa\langle\hat{a}(t)\rangle-e^{i\phi}\sqrt{G\kappa_{1}\kappa_{2}}\langle\hat{a}(t-\tau)\rangle-\gamma_{non}\langle\hat{a}^{\dagger}\hat{a}\hat{a}\rangle,$$

(20)

which is similar to the nonlinear DDE (18). However, there is a crucial difference with $\langle\hat{a}^{\dagger}\hat{a}\hat{a}\rangle\nsim x^{3}(t)$. This difference will inexorably introduce phase diffusion into the oscillation of $\langle\hat{a}\rangle$ in nonlinear quantum systems.

With the master equation (19), we can choose $G=1.2,\kappa_{1}=\kappa_{2}=\kappa,\gamma_{non}/\kappa=1$ and any delay time larger than $\tau_{cr}$. Here we choose $\kappa\tau=\kappa\left(\tau_{cr}+1\right)\approx 6.284$, to simulate the evolution of $\langle\hat{a}\rangle$. Unlike the results in classical Equation (18), the evolution of $\langle\hat{a}\rangle$ in the nonlinear quantum system shows dissipative oscillation, primarily due to quantum phase diffusion arising from the $\gamma_{non}\langle\hat{a}^{\dagger}\hat{a}\hat{a}\rangle$ term.

Due to the computer memory limitations, we can only simulate $4$ time steps in this case. We are still interested in how to simulate over longer times, even though we cannot see a perfect self-oscillation of $\langle\hat{a}\rangle$. Therefore, we also apply a similar mean-field approximation approach to explore the extended time evolution. Like the linear quantum case, we can focus on the ODEs of $\langle\hat{a}\rangle$ instead of the whole master equation, but the nonlinear case will be much more complex due to the nonlinear damping, since ODEs for $\langle\hat{a}^{\dagger}\hat{a}\hat{a}\rangle$ will be involved with higher order operator products. This results in an infinite number of ODEs involving these higher order operator products. We can truncate these equations by expanding the expectation of higher order operator products in terms of lower order products. This is the original concept of cumulants [39, 40], and we expect that truncating to higher order quantum correlations will result in more precision in the simulation. The details to apply the cumulants can be found in the C.

When we increase to larger time steps, the number of ODEs will increase. However, the resources required is far fewer than those needed to solve the full quantum master equation (19). Since the numbers of ODEs required increase rapidly with both time and truncation level, we need to balance the accuracy of this mean-field approximation and the number of ODEs required. We examine the accuracy of this method in Figure 6-(c). We observe that convergence of truncation of order $3$ seems excellent. We set the order as $3$ to simulate the effect of feedback gain $G$ on $\langle\hat{a}\rangle$ in Figure 6-(a). The results in Figure 6-(b), indicate that $\langle\hat{a}^{\dagger}\hat{a}\rangle$, $\Delta X$ and $\Delta P$ no longer rises in an unbounded fashion but both appear to saturate. We also examine the convergence of photon number of different expansion order to the full quantum master equation in Figure 6-(d). Finally, from extensive numerical analysis of this nonlinear SSOs system we could not find conditions which permitted indefinite SSO of $\langle\hat{a}(t)\rangle$ without phase diffusion.