Steady states, squeezing, and entanglement in intracavity triplet down conversion

TL;DR

This paper analyzes intracavity triplet down conversion using quantum phase-space methods, revealing steady states, squeezing, and entanglement, with quantum effects prominent near the threshold and validated by simulations.

Contribution

It introduces a quantum phase-space analysis of intracavity triplet down conversion, highlighting quantum effects and steady-state properties near the threshold.

Findings

Quantum effects are significant just above the pumping threshold.

Regimes of measurable squeezing and entanglement are identified.

Good agreement with Monte Carlo simulations for low populations.

Abstract

Triplet down conversion, the process of converting one high-energy photon into three low-energy photons, may soon be experimentally feasible due to advances in optical resonator technology. We use quantum phase-space techniques to analyse the process of degenerate intracavity triplet down conversion by solving stochastic differential equations within the truncated positive-P representation. The time evolution of both intracavity mode populations are simulated, and the resulting steady-states are examined as a function of the pump intensity. Quantum effects are most pronounced in the region immediately above the semi-classical pumping threshold, where our numerical results differ significantly from semi-classical predictions. Regimes of measurable squeezing and bipartite entanglement are identified from steady-state spectra of the cavity output fields. We validate the truncated positive-P description against Monte Carlo wave function simulations, finding good agreement for low mode populations.

Equations

$${\hat{H}}_{\text{int}}=i\hslash \frac{\kappa }{3}\left[{\hat{a}}^{†3}\hat{b}-{\hat{a}}^3{\hat{b}}^{†}\right],$$

$${\hat{H}}_{\mathrm{pump}}=i\hslash \left[{ϵ}_b{\hat{b}}^{†}-{ϵ}_b^{\ast }\hat{b}\right],$$

$$\mathcal{L}\rho ={\gamma }_a[2\hat{a}\rho {\hat{a}}^{†}-{\hat{a}}^{†}\hat{a}\rho -\rho {\hat{a}}^{†}\hat{a}]+{\gamma }_b[2\hat{b}\rho {\hat{b}}^{†}-{\hat{b}}^{†}\hat{b}\rho -\rho {\hat{b}}^{†}\hat{b}],$$

$$\dot{\rho }=-\frac{i}{\hslash }[{\hat{H}}_{\mathrm{sys}},\rho ]+\frac{1}{2}\mathcal{L}\rho ,$$

$$\displaystyle\begin{split}\frac{\partial P}{\partial t}=\Bigg{\{}&-\Bigg{[}\frac{\partial}{\partial\alpha}(\kappa\alpha^{\!+2}\beta)+\frac{\partial}{\partial\alpha^{\!+}}(\kappa\alpha^{2}\beta^{+})+\frac{\partial}{\partial\beta}\left(-\frac{\kappa}{3}\alpha^{3}\right)\\ &+\frac{\partial}{\partial\beta^{+}}\left(-\frac{\kappa}{3}\alpha^{\!+3}\right)-\frac{\partial}{\partial\alpha}(\gamma_{a}\alpha)-\frac{\partial}{\partial\alpha^{\!+}}(\gamma_{a}\alpha^{\!+})\\ &+\frac{\partial}{\partial\beta}(\epsilon_{b})+\frac{\partial}{\partial\beta^{+}}(\epsilon_{b}^{*})\Bigg{]}\\ &+\frac{1}{2}\Bigg{[}\frac{\partial^{2}}{\partial\alpha^{2}}(2\kappa\alpha^{\!+}\beta)+\frac{\partial^{2}}{\partial\alpha^{\!+2}}(2\kappa\alpha\beta^{+})\Bigg{]}\\ &-\frac{1}{6}\Bigg{[}\frac{\partial^{3}}{\partial\alpha^{3}}(2\kappa\beta^{+})+\frac{\partial^{3}}{\partial\alpha^{\!+3}}(2\kappa\beta)\Bigg{]}\Bigg{\}}\\ \times&\ P(\alpha,\alpha^{\!+}\!,\beta,\beta^{+}),\end{split}$$

$$\begin{array}{rl}d\alpha & =\left(-{\gamma }_a\alpha +\kappa {\alpha }^{+2}\beta \right)dt+\sqrt{2\kappa {\alpha }^{+}\beta }d{W}_1(t),\\ d{\alpha }^{+}& =\left(-{\gamma }_a{\alpha }^{+}+\kappa {\alpha }^2{\beta }^{+}\right)dt+\sqrt{2\kappa \alpha {\beta }^{+}}d{W}_2(t),\\ d\beta & =\left({ϵ}_b-{\gamma }_b\beta -\frac{\kappa }{3}{\alpha }^3\right)dt,\\ d{\beta }^{+}& =\left({ϵ}_b^{\ast }-{\gamma }_b{\beta }^{+}-\frac{\kappa }{3}{\alpha }^{+3}\right)dt,\end{array}$$

$$\begin{array}{rl}\frac{\partial \alpha }{\partial t}& =\kappa {\alpha }^{\ast 2}\beta -{\gamma }_a\alpha ,\\ \frac{\partial \beta }{\partial t}& ={ϵ}_b-{\gamma }_b\beta -\frac{\kappa }{3}{\alpha }^3,\end{array}$$

$$|{\alpha }_s{|}^4-\frac{3|{ϵ}_b|}{\kappa }|{\alpha }_s|+\frac{3{\gamma }_a{\gamma }_b}{{\kappa }^2}=0,$$

$${\beta }_s=\frac{{ϵ}_b-\kappa {\alpha }_s^3/3}{{\gamma }_b},$$

$$\kappa =\frac{3\hslash {\epsilon }_0{\chi }^{(3)}\sqrt{{\omega }_a^3{\omega }_b}}{4\sqrt{{\epsilon }_a^3{\epsilon }_b}}\frac{\delta [m_b-3m_a]}{L}\sigma ,$$

$$\sigma =\frac{\int u_a^3(x,y)u_b^{\ast }(x,y)dxdy}{{\left(\int |u_a(x,y){|}^{2}dxdy\right)}^{\frac{3}{2}}{\left(\int |u_b(x,y){|}^{2}dxdy\right)}^{\frac{1}{2}}}.$$

$$\mathbf{\alpha }(t)={\overline{\mathbf{\alpha }}}_s+\delta \mathbf{\alpha }(t).$$

$$d\delta \mathbf{\alpha }=-A\delta \mathbf{\alpha }dt+Bd\mathbf{W},$$

$$A=\left[\begin{array}{cccc}{\gamma }_a& -2\kappa {\alpha }_s^{\ast }{\beta }_s& -\kappa {\alpha }_s^{\ast 2}& 0\\ -2\kappa {\alpha }_s{\beta }_s^{\ast }& {\gamma }_a& 0& -\kappa {\alpha }_s^2\\ \kappa {\alpha }_s^2& 0& {\gamma }_b& 0\\ 0& \kappa {\alpha }_s^{\ast 2}& 0& {\gamma }_b\end{array}\right],$$

$$B{B}^T=D=\left[\begin{array}{cccc}2\kappa {\alpha }_s^{\ast }{\beta }_s& 0& 0& 0\\ 0& 2\kappa {\alpha }_s{\beta }_s^{\ast }& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right].$$

$$G(\tau )=⟨\delta \mathbf{\alpha }(t+\tau )\delta {\mathbf{\alpha }}^{†}(t)⟩,$$

$$S(\omega )=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }{e}^{-i\omega \tau }G(\tau )d\tau ,$$

$$S(\omega )=(A+i\omega I{)}^{-1}D(A^{†}-i\omega I{)}^{-1},$$

$$S^q(\omega )=QS(\omega )Q^T,\text{where}Q=\left[\begin{array}{cccc}1& 1& 0& 0\\ -i& i& 0& 0\\ 0& 0& 1& 1\\ 0& 0& -i& i\end{array}\right].$$

$$\begin{array}{rl}\mathrm{C}[{\hat{X}}_i,{\hat{X}}_j](\omega )& ={\delta }_{ij}+\sqrt{{\gamma }_i{\gamma }_j}(S_{2i-1,2j-1}^q+S_{2j-1,2i-1}^q),\\ \mathrm{C}[{\hat{Y}}_i,{\hat{Y}}_j](\omega )& ={\delta }_{ij}+\sqrt{{\gamma }_i{\gamma }_j}(S_{2i,2j}^q+S_{2j,2i}^q),\end{array}$$

$${\mathrm{DS}}_{\pm }\equiv \mathrm{V}({\hat{X}}_i\pm {\hat{X}}_j)+\mathrm{V}({\hat{Y}}_i\mp {\hat{Y}}_j)\ge 4.$$

$$\overline{\alpha +{\alpha }^{+}}$$

$$\overline{1+2\alpha {\alpha }^{+}+{\alpha }^2+{\alpha }^{+2}}$$

Full text

I Introduction

Intracavity pair down conversion, whereby pump photons are parametrically converted into pairs of photons inside a $\chi^{(2)}$ nonlinear medium, is a well studied process due to its many applications, particularly for allowing widely-tunable light sources—even as wide as two octaves Olsen (2017); Kolker et al. (2018); Olsen (2018); Li and Olsen (2018)—and squeezed light generation Fürst et al. (2011); Peano et al. (2015); Andersen et al. (2016). The degenerate case, $2\omega\rightarrow\omega+\omega$, where the generated photons are of the same polarization, is notable for being the leader in squeezed light generation, with a record of more than $15~{}\mathrm{dB}$ squeezing Vahlbruch et al. (2016). Above the optical parametric oscillation threshold, when spontaneous pair production is frequent enough to seed a large-scale transfer of energy between the modes, the same system still exhibits squeezing Fürst et al. (2011), as well as correlation in the output beams Takei et al. (2005); Jing et al. (2003). Degenerate triplet down conversion, $3\omega\rightarrow\omega+\omega+\omega$, is the third-order ($\chi^{(3)}$) equivalent of pair down conversion, and the inverse of third-harmonic generation Cavanna et al. (2016). This process has not yet been observed experimentally, neither above nor below threshold. Emerging resonator technologies combining low losses with large optical overlaps and nonlinearities are bringing experimental realization closer to fruition, motivating deeper theoretical investigation.

In this work we use a truncated form of Drummond and Gardiner’s positive-P representation Drummond and Gardiner (1980) to degenerate triplet down conversion in quantum phase space. The resulting equations are solved numerically to determine the time evolution of the intracavity fields of a damped, driven optical cavity, shown schematically in Fig. 1. Of particular interest are the resulting steady-state fields, which are examined over a range of pump intensities, including a comparison to analytical semi-classical predictions Bajer (1991); Felbinger et al. (1998). We observe good agreement except in the region immediately above the semi-classical pumping threshold where critical fluctuations are signficant. We examine this region for a range of initial conditions, and weak injected signals in the low-energy mode. We use a linearized fluctuation analysis to numerically identify regimes of squeezing and bipartite entanglement in the output fields. We perform a comparison of the time evolution given by our phase space methods with that of a more exact Monte Carlo wave function simulation in a number state basis. The results demonstrate the validity of the truncated positive-P method for this system.

The structure of this paper is as follows: in Sec. II we develop the equations of motion for truncated positive-$P$ and semiclassical approaches, and present the known semi-classical steady state solutions. In Sec. III we estimate experimental nonlinearities to justify our choice of simulations parameters, and solve for steady states of the system comparing positive-$P$ and semiclassical approaches. In Sec. IV we present numerical results for measures of squeezing and entanglement in the system, and in Sec. V we discuss our results and conclude.

II System and Equations of Motion

The system we consider is that of degenerate intracavity triplet down conversion, shown schematically in Fig. 1. The interaction is described by the Hamiltonian

[TABLE]

where $\hat{b}^{\dagger}$, $\hat{b}$ [$\hat{a}^{\dagger}$, $\hat{a}$] are the bosonic creation and annihilation operators for the high [low] energy photons, and $\kappa\propto\chi^{(3)}$ is the effective non-linearity of the cavity medium, which can be taken to be real and positive with no loss of generality. The high-energy mode is driven by a classical laser field described by the additional Hamiltonian

[TABLE]

where $\epsilon_{b}$ is the pump amplitude of the high-energy field. Loss of both high and low energy photons from the cavity is represented by the Liouvillian superoperator acting on the system density matrix

[TABLE]

where $\gamma_{a}$ [$\gamma_{b}$] is the cavity loss rate of the low [high] energy mode. The overall time evolution of the density matrix is given by the Linblad master equation Lindblad (1976); Daley (2014)

[TABLE]

where $\hat{H}_{\mathrm{sys}}=\hat{H}_{\mathrm{int}}+\hat{H}_{\mathrm{pump}}$. Using the positive-P operator correspondences Gardiner and Zoller (2014) we map (4) onto a partial differential equation for the positive-P distribution

[TABLE]

where the positive-P variables $\alpha$ and $\alpha^{\!+}$ are independent and complex-valued. They correspond to the operators $\hat{a}$ and $\hat{a}^{\dagger}$ in the sense that expectation values of normally ordered operator moments can be calculated via $\langle\hat{a}^{\dagger m}\hat{a}^{n}\rangle=\int d^{2}\alpha\ d^{2}\alpha^{\!+}\ \alpha^{\!+m}\alpha^{n}P(\alpha,\alpha^{\!+}\!,\beta,\beta^{+})$, and analogously for the $b$ mode. Due to the cubic terms in (1), third-order derivatives exist in (5), meaning it cannot be interpreted as a true Fokker–Plank equation (FPE). However, a system size expansion of (5) shows that the third-order terms scale with $n^{-1/2}$, where $n$ is the number of photons in the cavity. We can hence reasonably expect these terms to become negligible for large mode populations. This approximation is further examined in Appendix B. Truncating the third-order terms from (5) results in a true FPE. Since the diffusion matrix is positive definite, the FPE can be mapped to a set of coupled stochastic differential equations (SDEs) using the rules of Itô calculus:

[TABLE]

where the Wiener increments $dW_{i}$ are Gaussian distributed with $\overline{dW_{i}(t)}=0$ and $\overline{dW_{i}(t)dW_{\!j}(t)}=dt$. Stochastically averaging over trajectories of (6) yields moments approximately consistent with (5), i.e. $\overline{\alpha^{\!+m}\alpha^{n}}\rightarrow\langle\hat{a}^{\dagger m}\hat{a}^{n}\rangle$. These Itô SDEs are related to a set of equivalent Langevin equations, which can be obtained by dividing both sides of (6) by $dt$. The benefit of the Itô SDEs is that the Wiener increments are mathematically rigourously defined, allowing for well defined numerical integration. Each of the terms present has an easily interpretable physical meaning. The $\epsilon$ terms correspond to a constant pump, which continuously increases the population of the $b$ mode. The $\gamma$ terms generate damping by removing photons from each mode at a rate proportional to the mode population. The last non-stochastic term in each equation describes the change in the fields due to the $\chi^{(3)}$ interaction, proportional to the effective interaction strength $\kappa$. The interaction leads to a decrease in the high energy mode population, and an increase in the low energy mode. We lack a factor of $\frac{1}{3}$ in the $\alpha$ and $\alpha^{\!+}$ equations because down conversion of one $b$ photon leads to three $a$ photons. The appearance of e.g. the product $\kappa\alpha\beta^{+}$ in the stochastic term indicates that the stochisticity originates from the quantum nonlinear interaction.

It should be noted that, using the techniques developed in Plimak et al. (2003), it is possible to develop stochastic difference equations corresponding to (5), which, although still approximate, do not require truncation of third-order terms. However, these equations are unsatisfactory for practical use, as convergence cannot be proven and they suffer from severe numerical instabilities Olsen et al. (2002).

A further, semi-classical, approximation is obtained by neglecting the noise terms in (6). The convenient choice of initial condition $\alpha^{\!+}(0)=\alpha_{0}^{*}$ [$\beta^{+}(0)=\beta_{0}^{*}$], where $\alpha_{0}=\alpha(0)$ [$\beta_{0}=\beta(0)$], yields the coupled ordinary differential equations (ODEs)

[TABLE]

with $\alpha^{\!+}$ recovered via $\alpha^{*}\to\alpha^{\!+}$, and analagously for $\beta^{+}$. The steady-states of (7) have been well studied analytically as a function of $\epsilon_{b}$ in the physically relevant $\gamma_{a}\ll\gamma_{b}$ regime Bajer (1991); Felbinger et al. (1998). The trivial solution $\alpha=0$, $\beta=\epsilon_{b}/\gamma_{b}$ is a steady-state at all pump intensities. Steady state solutions with non-zero down converted fields are only realized above the pumping threshold $\epsilon_{b,\mathrm{th}}=4(\gamma_{a}\gamma_{b})^{3/4}/(3\!\sqrt{\kappa})$. The magnitude of these non-trivial solutions is determined by solving

[TABLE]

and each solution is triply degenerate in phase, which is independently determined by solving $e^{i3\theta}=\epsilon_{b}/|\epsilon_{b}|$. For a given set of parameters above the pumping threshold there are two realisable solutions to (8), which take the form of a pair of diverging branches when plotted against $\epsilon_{b}$. The lower branch is unstable, while the upper branch is stable for $\gamma_{a}<\gamma_{b}$ Felbinger et al. (1998). For each steady-state solution for $\alpha$, the corresponding solution for $\beta$ is given by

[TABLE]

determined by the balance between pumping, damping, and down/up conversion.

III Time Evolution and Steady States

III.1 Experimental considerations

We have to make a choice for the regime for our simulations, specifically the ration of damping to nonlinear interaction; this ratio sets the size of signficant quantum effects caused by the triplet down conversion. We motivate our choice of interaction physically by providing an estimate for $\kappa$ in a monolithic resonant cavity. Assume the fields are linearly polarized within a resonant cavity of length $L$ entirely made of a nonlinear material with effective nonlinearity $\chi^{(3)}$ and permittivity $\varepsilon_{a}$ [$\varepsilon_{b}$] for the low- [high-]energy mode. We can then perform a macroscopic quantisation of the displacement fields Drummond and Hillery (2014), and substitute the quantized fields into the interaction Hamiltonian. We find for $\kappa$:

[TABLE]

where $\varepsilon_{0}$ is the vacuum permittivity, $m_{a}$ [$m_{b}$] is the number of field oscillations of the low- [high-]energy mode in one round trip, and $\sigma$ is the transverse modal overlap given by

[TABLE]

Here $u_{a}$ [$u_{b}$] is the transverse mode profile of the displacement field of the low- [high-] energy mode, which is assumed to be constant along the length of the cavity. The Kronecker-delta term $\delta[m_{b}-3m_{a}]/L$ originates under this assumption from the overlap of the mode profiles in the direction of propagation, and asserts that the modes are phasematched.

From the most efficient recent demonstration of doubly-resonant third-harmonic generation we are aware of Surya et al. (2018), we can extract the following parameters: $\omega_{a}=2\pi\times$194\text{,}\mathrm{THz}; $L=2\pi\times$20\text{\,}\mathrm{\SIUnitSymbolMicro m}; $\sigma=$0.43\text{,}{\mathrm{\SIUnitSymbolMicro m}}^{-2} Surya and Tang ([2019](#bib.bib22)); $\chi^{(3)}\approx$1.5\text{\times}{10}^{-20}\text{\,}{\mathrm{m}}^{2}\text{\,}{\mathrm{V}}^{-2}, taking the average of values for silicon nitride Ning et al. (2013) and aluminium nitride Mateen et al. (2018); and $\varepsilon_{a,b}\approx 4\varepsilon_{0}$ Kischkat et al. (2012). From these parameters, we can estimate that $\kappa$ is of the order of magnitude of $100\text{\,}{\mathrm{s}}^{-1}$, which is seven orders of magnitude smaller than their measured loss rate of $\gamma_{a}=\omega_{a}/(2Q_{a})=$1.5\text{\times}{10}^{9}\text{,}{\mathrm{s}}^{-1}$$, where $Q_{a}$ is the $Q$-factor of the low-energy mode. However, numerically we also have a limitation posed by finite resources. Thus, in all of the presented results we set $\kappa=0.001\gamma_{a}$, so that the interaction remains much weaker than dissipation, but its effects can be more easily studied in simulations. A further advantage is that the artifacts of truncation that we consider further in Appendix B are exaggerated by this choice (and our choice of lower photon number), offering a pessimistic assessment of the vality of truncation.

III.2 Simulations

We use the adaptive LambaEulerHeun Rackauckas and Nie (2016) algorithm implemeted in the Julia Bezanson et al. (2017) package DifferentialEquations.jl Rackauckas and Nie (2017) to simulate the truncated positive-P SDE (6). Expectation values of the mode populations $\langle\hat{n}_{a}\rangle=\langle\hat{a}^{\dagger}\hat{a}\rangle$ and $\langle\hat{n}_{b}\rangle=\langle\hat{b}^{\dagger}\hat{b}\rangle$ were calculated by averaging over $10^{6}$ trajectories, which was found to be sufficient for convergence of the second-order moments used to characterize the system throughout this work. For comparison with a semi-classical treatment, mode populations were obtained from the mean-field equations (7) using an ODE solver.

We begin by examining the time evolution of the mode populations for typical experimental parameters. We set the initial state to be the vacuum, since for optical photons the thermal occupation of the modes will be negligible. Due to the nature of the positive-P function as a pseudo-probability distribution, there are an infinite number of corresponding choices of distributions for $\bm{\alpha}\equiv(\alpha,\alpha^{\!+},\beta,\beta^{+})$. The simplest positive-P distribution for the vacuum is a coherent state with zero amplitude, and the distribution is a delta function at $\bm{\alpha}=\mathbf{0}$. Inspection of (6) shows that down conversion will not proceed in this case because the noise term in $\alpha$ [$\alpha^{\!+}$] is proportional to the product $\alpha^{\!+}\beta$ [$\alpha\beta^{+}$], which will always be zero. We might guess that a different representation of the vacuum that involves distribution including $\bm{\alpha}\neq\mathbf{0}$ Olsen and Bradley (2009) could allow for a mechanism that seeds spontaneous triplet production. However, find numerically that this is not the case, consistent with the fact that the physical result should be independent of the particular positive-P representation. The true source of the inhibition of spontaneous processes is that the truncation of third order derivatives in (5) is not valid for small populations, due to the invalidity of the scaling argument given above. Thus, we are unable to describe true spontaneous processes within a truncated positive-P approach. Inclusion of the third-order terms would allow for spontaneous triplet production Olsen et al. (2002), but at prohibitive technical cost. In contrast, the positive-P SDE for spontaneous pair production requires no truncation, and has noise terms proportional to only $\beta$ and $\beta^{+}$, which seed the process Olsen (2004).

Here we introduce a weak injected signal $\epsilon_{a}\ll\epsilon_{b}$ to the equation for $\alpha$ (and $\epsilon_{a}^{*}$ for $\alpha^{\!+}$). This approach is also likely to be used in experiments as the low efficiency of the process all but prohibits triplet production, and it remains unclear whether continuous down conversion is even possible without a seed. Using this aproach, the expected mode populations as given by (6) and (7) are compared in Fig. 2(a). Inclusion of noise significantly reduces the magnitude of the oscillatory behaviour in the transient regime, but both models eventually reach the same steady-state. Removing the injected signal once the steady-state is reached, as shown in Fig. 2(b) and (c), leads to small transients and a slightly altered steady-state. This steady-state is consistent with semi-classical values of $|\alpha_{s}|^{2}$ and $|\beta_{s}|^{2}$, determined analytically from (8) and (9).

The time evolution and eventual steady-state presented in Fig. 2 corresponds to a particular choice of the pump intensity $\epsilon_{b}$. The corresponding steady-states of $|\alpha|$, as given by (6) and by (7), are compared in Fig. 3 as a function of the pump intensity $\epsilon_{b}$. The analytic steady-states of (7), discussed above for $\gamma_{a}\ll\gamma_{b}$, are replicated in Fig. 3 (red lines) for $\gamma_{a}\not\ll\gamma_{b}$, using the numeric steady-state solver from DifferentialEquations.jl Rackauckas and Nie (2017) acting on (7). Direct simulation of (7) produces results consistent with the steady-state solver, but the transition from $\alpha=0$ to the upper branch will depend on the initial state and/or any injected signal. A single realisation is presented in Fig. 3 (black crosses). The steady-states of $|\alpha|$ given by simulations of the truncated positive-P SDE (6) are presented in Fig. 3 for four different initial states (cyan triangles, blue squares, orange diamonds, green circles). As above, these simulations use an injected signal, $\epsilon_{a}$, which is turned off part way through the simulations. The results match the semi-classical steady-states for most pump intensities, except in the region immediately above the semi-classical pumping threshold, referred to henceforth as the transition region. In this region the steady-state populations take on values intermediate to, and form a connection between, the stable semi-classical solutions. Individual trajectories of (6) in the transition region end up at either $\alpha=0$, or stochastically fluctuating about the upper branch; the expectation values presented indicate the ratio of trajectories that go up or down, and are dependent on the initial state (as shown in Fig. 3) and/or the presence of any injected signal (see Fig. 4). It is important to emphasize that only operator moments can be recovered from these simulations, so the individual trajectories are not physically meaningful.

We now briefly consider the steady-states in the case that we do not turn off the injected signal, which has so far only been used to initiate down conversion. Fig. 4 presents the steady-states around the transition region for systems with a range of injected signal intensities. Values given by simulations of (6) are presented alongside the (stable and unstable) semi-classical results of (7), obtained using a numeric steady-state solver. For large injected signals there is only a single stable semi-classical solution at any given pump intensity. As the injected signal is decreased below some critical intensity, this solution morphs into what will eventually be the upper branch, the (unstable) lower branch, and part of the $\alpha=0$ solution. The other part of the $\alpha=0$ solution is in fact a separate solution which is only realized below the critical injected signal. Simulations of (6) are largely in agreement with these semi-classical results. Below the critical injected signal we observe the same behaviour in the transition region as in the $\epsilon_{a}=0$ results from Fig. 3, i.e. steady-states intermediate to the two stable semi-classical solutions. We emphasize that the steady-states will depend on the initial state of the system. In particular, steady-states intermediate to the two semi-classically stable solutions will not always exist; parameters were specifically chosen for Fig. 4 to show that they are possible.

IV Squeezing and entanglement

The quantities of interest in quantum information applications are usually the steady-state fluctuations of the output fields. It is standard quantum optical practice to consider the steady-state spectra, which are calculated via Fourier transform of the two-time covariance matrix. This approach is outlined in Gardiner (2009), having been originally developed for stochastic analysis of chemical reactions by Chaturvedi et al. Chaturvedi and Gardiner (1977). We begin by decomposing the phase space variables into their steady-state expectation values and a fluctuation term,

[TABLE]

Substituting (12) into (6) and linearising, the resulting equation of motion is the multivariate Ornstein–Uhlenbeck process

[TABLE]

with

[TABLE]

and

[TABLE]

The process of linearization assumes the fluctuations are small relative to their respective expectation values, which holds for all pump intensities except for those in the transition region. This region is examined in more detail in Appendix A. The Ornstein–Uhlenbeck processes (13) is stationary, so we can define the two-time covariance matrix as a function of a single variable,

[TABLE]

and in turn define the Spectrum matrix,

[TABLE]

where $\omega$ is an angular frequency. Using standard results Gardiner (2009), we can calculate the Spectrum matrix directly from (14) and (15) via

[TABLE]

where $I$ is the identity matrix. It is important to note that this process is not valid if the eigenvalues of $A$ have any negative real parts, as this will result in fluctuations increasing without bound. However, this is not the case for any of the results presented.

In order to express the spectrum in terms of the canonical quadrature operators $\hat{X}_{a}=\hat{a}^{\dagger}+\hat{a},\hat{Y}_{a}=i(\hat{a}^{\dagger}-\hat{a})$, and those defined analogously for the $b$ mode, we calculate Walls and Milburn (1994)

[TABLE]

The output spectra is calculated using the standard input–output relations Gardiner and Collett (1985), hence we can calculate output spectral variances and covariances via

[TABLE]

where we use $a\rightarrow 1$ and $b\rightarrow 2$ for the $i,j$ indexing, required to be cyclic in all indices.

The non-classical nature of the output fields can be demonstrated in various ways. One of the common measures is whether the uncertainty in one quadrature is less than the vacuum (or coherent state) level, in which case the output is said to be squeezed Caves (1981); Walls (1983). In the present system we find squeezing of the $\hat{Y}$ quadratures, shown in Fig. 5 for the simulation presented in Fig. 2(a). Also of interest is entanglement and inseperability of the two output fields, often demonstrated by violation of the Duan–Simon inequality Duan et al. (2000); Simon (2000), which can be written as

[TABLE]

We find that the modes are entangled, with $D_{-}$ clearly violated for $\omega\rightarrow 0$, as shown in Fig. 6. Violation of the $\mathrm{DS}_{-}$ inequality is also predicted at higher frequencies, but at a magnitude that is unlikely to be experimentally measurable. It should be noted that we are dealing with a non-Gaussian system, so in regions where (21) is satisfied we cannot necessarily claim that the fields are not entangled.

The steady-states obtained from the simulations examined in Fig. 3 were used to investigate squeezing and bipartite entanglement over a range of pump intensities. The linearized fluctuation analysis applied here cannot be applied to the transition region due to the magnitude of the fluctuations. The squeezing results are presented in Fig. 7. Below threshold the low-energy mode is in a vacuum state, and the classically pumped high energy mode is in a coherent state, so the uncertainty is at the vacuum level. Above threshold we observe squeezing at all pump intensities up to $\epsilon_{b}/\gamma_{a}=300$, the largest value examined. The spectra have the qualitative shape of those presented in Fig. 5. As the pump intensity increases, the maximal squeezing moves towards higher frequencies and decreases in magnitude. The key difference between the high- and low-energy modes is that near the transition region the high energy ($b$) mode exhibits maximal squeezing at zero-frequency. The $b$ mode also exhibits stronger squeezing at all pump intensities, and the magnitude of its squeezing decreases slower than that of the low energy mode as the pump intensity increases.

We find entanglement of the modes at all pump intensities above the transition region up to at least $\epsilon_{b}/\gamma_{a}=300$. The $DS_{-}$ inequality is always violated at large and small frequencies, but is satisfied between these two regions, just as in Fig. 6. $DS_{+}$ is never violated. The maximal violation of the $\mathrm{DS}_{-}$ inequality decreases in magnitude as the pumping is increased, indicating less entanglement. The peak in $DS_{-}$ [$DS_{+}$] moves towards larger frequencies and increases [decreases] in magnitude as the pumping is increased.

V Discussion and Conclusions

V.1 Discussion

The requirement to seed the process at $\omega_{a}$ is a significant experimental challenge. The seed must not only be as close as possible to satisfying $\omega_{a}=\omega_{b}/3$, but it must also be phase stable and have the right phase relative to the pump laser; a phase shift of $\pi/3$ causes the gain in the non-trivial solution to turn negative and the signal will decay, similar to the behavior of a degenerate second-order optical parametric amplifier with a pump phase shift of $\pi/2$. This challenge can be mitigated by only seeding the interaction for a time long enough to move the system close to steady state, but short enough that the relative phase varies by less that $\pi/3$. It is worth noting that relative phase stability of the seed and pump does not pose a problem for non-degenerate third-order down conversion (not considered here), as phase wandering is compensated by complementary wandering in the unseeded modes. Although tripartite entanglement can be created via tuned pair downconversion Bradley et al. (2005), the non-degenerate system could also be useful as a more direct source of tripartite quantum correlations.

An alternative approach is to derive the pump from the third harmonic of the seed laser. This automatically phase locks the system provided the optical path length is stable, and additionally has the benefit that the seed laser can act as the local oscillator for an optical homodyne measurement. If the pump were generated by third-harmonic generation, however, it could be amplitude squeezed, and it is unclear how this would impact the phase squeezing created by triplet down conversion.

Our choice of relatively large nonlinearity $\kappa=0.001\gamma_{a}$ enabled a study of significant nonlinearity in the system, including the intermediate steady states observed in the transition region. In the context of making a comparison between truncated positive-$P$ and the number state basis (see Appendix B) our choice of $\kappa$ acts to over-estimate the discrepancy between the two methods, amplifying the role of nonlinearity relative to dissipation, compared to more readily accessible experimental parameters (see Sec. III.1). Despite this choice, we observe good agreement between the truncated positive-$P$ and Monte Carlo wavefunction methods.

V.2 Conclusion

We have performed a numerical investigation of intracavity degenerate triplet down conversion using a truncated positive-P representation. A systematic analysis of the stable steady-states of the expected mode populations found close agreement with semi-classical predictions, except in the region immediately above the pumping threshold. Interestingly, here we find that quantum stochastic simulations predict the possibility of steady states intermediate to the two stable semi-classical solutions.

A steady-state linearized fluctuation analysis was performed and used to calculate the intracavity fluctuation spectrum. From this, regimes of squeezing and bipartite entanglement were identified. Both of these quantities were found to be maximal closest to the pumping threshold. However, the linearisation involved in determining the spectrum is not valid immediately above the pumping threshold. An interesting future direction would be to carry out a systematic study of critical fluctuations, and their measurable extracavity signatures.

Acknowledgements.

AB and HS are supported by the Dodd-Walls Centre for Photonic and Quantum Technologies.

Appendix A Validity of Linearized Fluctuation Analysis

The Ornstein–Uhlenbeck process (13) used to determine the spectra presented in this work was derived by linearising the steady-state equations of motion under the assumption that fluctuations of the phase-space variables are negligible compared to their expectation values. This assumption is examined by considering the quadrature operators for both modes, rather than the phase-space variables themselves, because they are explicitly real-valued quantities. However, this introduces the question as to which quantities to compare to which. For the systems examined in this paper, the steady-state (quadrature) phase-space distributions for each mode are centred on the $X$ axis (i.e. $\langle\hat{Y}\rangle=0$), so we must consider the standard deviation of both the $X$ and $Y$ quadratures relative to $\langle\hat{X}\rangle$. These values are calculated by averaging over stochastic trajectories. For example,

[TABLE]

and then $\Delta X_{a}=\sqrt{\langle\hat{X}_{a}^{2}\rangle-\langle\hat{X}_{a}\rangle^{2}}$.

The results are presented over a range of pump intensities in Fig. 8(a) for the low energy mode, and (b) for the high energy mode. The quadrature fluctuations exhibits the same behaviour in both modes. While $\Delta Y$ is negligible everywhere, $\Delta X$ becomes comparable to $\langle\hat{X}\rangle$ within the transition region. In context of the phase-space variables, this corresponds to some trajectories going towards each of the two distinct stable semi-classical solutions, resulting in a large range of values. For this reason, we note that our method for determining the spectra is not valid within the transition region. It should be noted that just like the expectation value, the standard deviation is dependent on the initial conditions. While the results presented here correspond to a single set of parameters, they are representative of the general behavior of the system.

Appendix B Comparison of Truncated Positive-P and Monte Carlo Wave Function Techniques

The truncation involved in the mapping from (5) to (6) inevitably introduces a level of discrepancy, but will generally be more valid at higher intensities due to the scaling of the third-order terms with respect to the number of photons. This claim is examined via direct simulation of (4) using a Monte Carlo wave function (MCWF) technique Dalibard et al. (1992); Dum et al. (1992); Daley (2014) in a number state basis.

These simulations were implemented using the QuantumOptics.jl package Krämer et al. (2018). Optical phase space techniques are used in these kinds of systems specifically because they scale well with system size, while a number state basis does not. For this reason it is not feasible to run MCWF simulations with the parameters used in this paper. Instead, we compare MCWF and truncated positive-P simulations for a smaller system of $\sim\!100$ photons. The two techniques are compared in Fig. 9. Both above and below threshold the two techniques are in close agreement, although there is a slight difference in the period of the oscillations exhibited in (b).

For this reason, the effect of the truncation involved in deriving the SDEs used throughout this paper can be taken to be negligible. Both above and below threshold the MCWF simulations used Hilbert spaces with the dimensions $N_{a}=270$ and $N_{b}=90$ for the respective modes. It has been confirmed that increasing the size of these Hilbert spaces does not change the result.