Parametrically driving a quantum oscillator into exceptionality

The mathematical objects employed in physical theories do not always behave well. Einstein’s theory of space and time allows for spacetime singularities and Van Hove singularities arise in condensed matter physics, while intensity, phase and polarization singularities pervade wave physics. Within dissipative systems governed by matrices, singularities occur at the exceptional points in parameter space whereby some eigenvalues and eigenvectors coalesce simultaneously. However, the nature of exceptional points arising in quantum systems described within an open quantum systems approach has been much less studied. Here we consider a quantum oscillator driven parametrically and subject to loss. This squeezed system exhibits an exceptional point in the dynamical equations describing its first and second moments, which acts as a borderland between two phases with distinctive physical consequences. In particular, we discuss how the populations, correlations, squeezed quadratures and optical spectra crucially depend on being above or below the exceptional point. We also remark upon the presence of a dissipative phase transition at a critical point, which is associated with the closing of the Liouvillian gap. Our results invite the experimental probing of quantum resonators under two-photon driving, and perhaps a reappraisal of exceptional and critical points within dissipative quantum systems more generally.


Introduction
Conventionally, Hermitian physics has reigned supreme.Non-Hermitian extensions to standard theories were typically considered as mere perturbations, perhaps allowing for some small amount of dissipation to be accounted for.However, when a non-Hermitian matrix permits degeneracies in both its spectrum and its eigenfunctions at some special point -a so-called exceptional point (EP) [1] -everything changes.Indeed, the physics arising in the vicinity of an EP may be significantly different to anything occuring in the Hermitian version of the theory [2,3].
Inspired by the recent studies of nonlinear, Kerr-like resonators with two-photon driving [23][24][25][26][27][28][29][30][31], here we study a parametric driven-dissipative quantum oscillator.We place an emphasis on the squeezing and the phase-space representation of the quantum state of the oscillator, as well as both the first moments and the second moments of the system (including their steady state and transient behaviours).We link the EP arising in each case to a relevant observable quantity, including the mean populations [32,33], optical spectrum, degrees of coherence and squeezed quadratures.We also highlight the importance of a critical point in the system, which is associated with a dynamical instability and a dissipative phase transition.The effects of anharmonicities on the oscillator are briefly discussed, mostly in relation to its effect on the closing of the Liouvillian gap and the presence of the dynamical instability, since a full treatment of such a nonlinearity can be found in Refs.[23][24][25][26][27][28][29][30][31].
The studied driven-dissipative parametric system may be realized in modern quantum optical laboratories with judicious use of cavities and traps [34][35][36][37], or by designing certain types of nonlinear quantum circuits [38][39][40].The discovery of EPs in such a parametric system would join a growing list of quantum systems, mostly involving superconducting qubits but also including trapped-ion systems, which exhibit exceptional point physics [41][42][43][44][45][46][47][48].Already, EPs have been experimentally shown to be useful for the control of quantum states and even for the enhancement of quantum heat engines [49][50][51].
The parametric driven-dissipative quantum harmonic oscillator model, as sketched in Fig. 1 (a), may be described by the following Hamiltonian Ĥ (here and throughout we take = 1) [52] where, without loss of generality, the oscillator-driving detuning ∆ ≥ 0, the two-excitation driving amplitude Ω ≥ 0, and the driving phase −π ≤ θ ≤ π (see the Supplementary Information for more details).The creation and annihilation ladder operators  1) is possible with the aid of the bosonic Bogoliubov operator β, defined in terms of the squeezing parameter φ as which satisfies the commutator [β, β † ] = 1, and which is valid for sufficiently small driving amplitudes such that Ω < ∆ [53,54].This Bogoliubov transformation leads to the diagonalized form of the Hamiltonian Ĥ, complete with Bogoliubov mode eigenfrequency ω, as follows This brief analysis reveals a squeezed energy ladder with the eigenfrequencies E n = nω, which correspond to the squeezed number states |n, φ = S φ |n , where the squeezing operator S φ = exp 1 2 φe −iθ bb − 1 2 φe iθ b † b † .The squeezed oscillator levels E n are separated by the spacing ∆ at small driving amplitudes (namely, the harmonic oscillator limit with Ω ∆), while the inter-level spacing is vanishing in the limit of large driving amplitudes (Ω → ∆).Notably, there is so-called "spectral collapse" for Ω > ∆, where the eigenfrequencies E n become complex [in keeping with the construction of the squeezing parameter φ in Eq. ( 2)], however we do not enter this regime within the closed system version of our theory as solely described by Eq. (1).
However, the squared counterparts 0, φ| X2 |0, φ = e −2φ /2 and 0, φ| P 2 |0, φ = e 2φ /2, highlighting the quadrature squeezing as governed by the squeezing parameter φ [cf.Eq. ( 2)].These closed system results follow from the inverse relation of Eq. ( 2), where b = cosh (φ) β − e iθ sinh (φ) β † .We include dissipation in the model via an open quantum systems approach.Employing the Born, Markov and secular approximations, the associated quantum master equation for the system's reduced density matrix ρ may be given by [55,56] where γ ≥ 0 is the dissipation rate and the Hamiltonian Ĥ is given by Eq. ( 1).This quantum model includes the unitary evolution of the density matrix ρ via the the first term on the right-hand-side of Eq. ( 4) -the Liouville-von Neumann equation -while non-unitary evolution is accounted for by the second term in Eq. ( 4), which upgrades Eq. ( 4) to the form of a Gorini-Kossakowski-Sudarshan-Lindblad equation.We take γ as a phenomenological parameter, but in the case of modelling a thermal bath γ is the zero-temperature decay rate, and we show in the Supplementary Information that the effects of nonzero temperatures are rather negligible (for example, the location of the EP is unchanged).Using the property Tr (Oρ) = O , which leads to the mean value O of any operator O, the averages of the first moments of the system (like b for example) may be readily found.For the presented parametric model, as defined by Eq. ( 4) supplemented with Eq. ( 1), this procedure leads to the following Schrödinger-like dynamical equation where ψ, the two-dimensional column vector which collects the first moments, and the dynamical matrix H are defined via The two complex eigenfrequencies ω ± arising from the 2 × 2 matrix effective Hamiltonian H read which generalizes the closed Hamiltonian Ĥ description of Eq. ( 3), with its wholly real energies ω and its restriction Ω < ∆.
In particular, we have introduced the effective decay rate Γ, defined via the relation [cf. the definition of the frequency ω from Eq. ( 3)] which captures the physics of system in its dissipative phase with purely imaginary eigenvalues ω ± .The two normalized eigenvectors α ± , corresponding to the eigenfrequencies ω ± , are given by Clearly, both Eq. ( 7) and Eq. ( 9) present a point of coalescence -a so-called EP -at the particular driving amplitude Ω = Ω EP , defined by At this special point the eigenfrequency ω = 0 [cf.Eq. ( 3)], and the effective Hamiltonian H in Eq. ( 5) becomes defective [2,3].Since the EP of Eq. ( 10) has arisen from the Hamiltonian-like analysis of the dynamical equation in Eq. ( 5), it can be considered to be a so-called 'Hamiltonian EP'.Using the eigenfrequencies of Eq. ( 7), we plot in Fig. 1 (b, c) the real and imaginary parts of ω ± as a function of the drive amplitude Ω (for an example case with the detuning ∆ = 3γ/2).The EP is marked by the vertical, dashed grey line in both panels and highlights the bifurcations in the complex eigenvalues ω ± , as well as the guarding of the boundary between two drastically different physical phases.In what follows, we shall consider the consequences of this EP for the physical responses of the squeezed oscillator system.

Results
Similar to the calculation leading to the coupled dynamical equations of Eq. ( 5), the second moments of the system (that is quantities like b † b , the average population of the oscillator) are defined through the non-homogeneous equation where the three-dimensional column vector Ψ, which gathers the second moments, the driving term P and the 3 × 3 matrix M read [cf.Eq. ( 6)] The three complex eigenvalues λ + , λ − and λ 3 of the dynamical matrix M are given by [cf.Eq. ( 7)] where ω and Γ are defined in Eq. ( 3) and Eq. ( 8) respectively.The three corresponding (and unnormalized) eigenvectors β 3 , β + and β − read (we assume Ω < Ω EP for the moment, the results for Ω > Ω EP are of course similar) The simultaneous coalescing of all three eigenvalues λ i and eigenvectors β i at the certain point Ω = Ω EP reveals that the EP in the second moments occurs at exactly the same point in parameter space as the EP in the first moments [cf.Eq. ( 13) and Eq. ( 14) with Eq. ( 10)].This suggests that the 'Hamiltonian EP' [arising from Eq. ( 5)] and 'Liouvillian EP' [arising from Eq. ( 11)] are identical for this simple system.Furthermore, the second moments EP may be classed as being of third-order (due to its emergence with the coalescing of three objects), while the first moments are associated with a second-order EP (only two objects coalesce in this case) [57,58].The higher-order EP associated with Eq. ( 11) is graphed in Fig. 2 (a, b) using the results of Eq. ( 13), where the real parts of the eigenvalues λ i are displayed in the upper panel, and the corresponding imaginary parts in the lower panel.This analysis shows the importance of exceptional point physics throughout the different levels of description of open quantum systems (for example with increasing large moments).

Average populations.
The mean population of the oscillator n(t) = b † b is contained within the first element of Ψ [cf.Eq. ( 12)] and is found by formally solving the equation of motion defined in Eq. (11).At long time scales (t → ∞), the steady state population n(∞) = lim t→∞ b † b readily follows from Eq. ( 11) by noting that in this limit ∂ t Ψ = 0 and so its steady state solution where ω is defined in Eq. ( 3) and Γ in Eq. (8).As a population, the denominator of Eq. ( 15) should be non-negative.This suggests that the presented system is stable for driving amplitudes Ω < Ω c , where the critical driving frequency Ω c is defined as Otherwise there is no steady state -the competition from the driving and the dissipation is conclusively won by the driving and the population increases in time without bound.In practice, this continual climbing of the infinite and bosonic energy ladder of the oscillator can be tamed by either truncating the oscillator or by considering anharmoncities (both cases are discussed later on, although we are not so concerned with driving amplitudes satisfying Ω ≥ Ω c since the EP has already been passed by this stage).Notably, the result of Eq. ( 16) was foreshadowed by the complex eigenvalues ω ± provided in Eq. ( 7).In particular, in the regime of Ω > ∆ one of the eigensolutions (ω − ) corresponds to exponentially growing behaviour in time when γ/2 < Γ, consistent with the relation of Eq. ( 16).
The steady state phase diagram of the system implied by Eq. ( 16), as a function of the detuning ∆ and the driving amplitude Ω, is shown in Fig. 1 (d).The green region occurs below the EP (Ω < ∆), while the red region arises above the EP (Ω > ∆).No steady state is able to form in the blue region, as suggested by Eq. ( 16), and hence the population here is unbounded with increasing time.Hence for vanishing detuning ∆ γ the steady state most likely exists in the 'above EP' red region, within the bound 0 < Ω < γ/2.Contrarily, with large detuning ∆ γ the steady state is most likely to be in the 'below EP' green region ∆ > Ω, with a smaller chance of being in the slither of 'above EP' red region when ∆ < Ω < ∆ 1 + γ 2 /8∆ 2 .The diagram of Fig. 1 (d) thus acts as a kind of map for the parametric oscillator, highlighting both the exceptional point at Ω = Ω EP and the critical point at Ω = Ω c .
The full solution of the equation of motion given in Eq. ( 11) leads to the oscillator population n(t), including both transient and steady state parts, which in general is given by for driving amplitudes Ω < Ω EP , where the frequency ω is defined in Eq. ( 3), and where n(0) is the population of the oscillator at the initial time t = 0. Similarly, for driving amplitudes Ω > Ω EP the analogous result, using the damping rate Γ as given in Eq. ( 8), is (18) The solution of Eq. ( 17) displays characteristic sinusoidal and cosinusoidal Rabi-like oscillations, along with an exponential decay with the time constant γ, until the driving amplitude Ω overcomes the detuning ∆ and the nonoscillatory solution of Eq. ( 18) supersedes it.However, when the driving amplitude is exactly Ω = Ω EP [cf.Eq. ( 10)] the solution of Eq. ( 17) is drastically reconstructed into the much simpler form Most notably, while the exponential decay is unaffected, the expression of Eq. ( 19) features both linear and quadratic terms in the dimensionless time γt (instead of this quantity only appearing in trigonometric or hyperbolic functions) due to the nature of the EP.We plot the average population n(t) in Fig. 2 (c) for several values of driving amplitude Ω, where the detuning is fixed at ∆ = 3γ/2 (so that Ω EP = 3γ/2 and Ω c = 5/2γ 1.58γ from Eq. ( 10) and Eq. ( 16) respectively).For small driving amplitudes Ω -well below the EP (dark green line) -the population of the oscillator never becomes significant and monotonically decreases to a low plateau.With increasing drive amplitudes Ω a non-monotonic behaviour develops in the manner of Rabi Populations of the parametric driven-dissipative oscillator.Panel (a): the real parts of the eigenvalues λi, arising from the second moments matrix M, as a function of the drive amplitude Ω (in units of the loss rate γ) [cf.Eq. ( 13)].Panel (b): the corresponding imaginary parts of λi.Dashed grey lines: the exceptional point (where Ω = ∆).Panel (c): the mean population of the oscillator n(t), as a function of time t (in units of γ −1 ) with the initial mean population n(0) = 1 due to the oscillator being in its first excited state [cf.Eq. ( 17)].We consider several values of Ω, including below the exceptional point (green lines) and above the exceptional point (red line).The result exactly at the exceptional point is given by the orange line [cf.Eq. ( 19)].Thin grey lines: intermediate values of Ω.In this figure, we consider the case of the detuning ∆ = 3γ/2, so that ΩEP = 3γ/2 [cf.Eq. ( 10)] and Ωc = 5/2γ 1.58γ [cf.Eq. ( 16)].
oscillations (lime green line), which is more pronounced for larger drives, for example approaching the EP (orange line).Above the EP (red line) the population surges to increasingly large steady state plateaus without oscillation.Finally, above the critical amplitude Ω c no steady state is formed such that a form of dynamical instability has arisen due to the driving overpowering the dissipation.

Correlation functions.
The impact of the EP is also felt in the n-th degree of coherence g (n) (τ ), of which we consider the first-order and secondorder degrees in particular [55,56].The normalized first-order correlation function g (1) (τ ), with some delay time τ ≥ 0, may be defined as g (1) (τ ) = lim t→∞ b † (t)b(t + τ ) / b † (t)b(t) .For the presented model, using the results for the first and second moments and the quantum regression formula, it is explicitly given by either a damped-trigonometric expression or a damped-hyperbolic expression (see the Supplementary Information for details) where the frequency ω and the damping rate Γ are defined in Eq. ( 3) and Eq. ( 8) respectively, and where the characteristic time constant γ/2 appears in the decay factor.Clearly, in the regime of smaller driving amplitudes where Ω < Ω EP , the first-order coherence g (1) (τ ) is damped out at long timescales to zero, lim τ →∞ g (1) (τ ) = 0.However, in the opposing regime where Ω > Ω EP , the asymptotics are instead described by which describes a convergent solution only when γ/2 > Γ, or equivalently when the driving amplitude Ω < Ω c [cf. Eq. ( 16)], which corresponds to when the system supports a steady state.Otherwise, a divergence in g (1) (τ ) with increasing time may be prevented by truncating the infinite energy ladder of the harmonic oscillator or by introducing anharmonicities.Notably, the formal definition of the normalized first-order correlation function g (1) (τ ), as given above, ensures that g (1) (0) = 1 for all values of the driving amplitude Ω.
At the intermediate case between those described in Eq. (20), that is when the EP is approached [cf.Eq. ( 10)], the firstorder correlation function g (1) (τ ) collapses [in a similar manner to the population dynamics reconstruction of Eq. ( 19)] into the (τ ) FIG. 3. Correlations of the parametric driven-dissipative oscillator.Panel (a): the degree of first-order coherence g (1) (τ ) as a function of the delay time τ , in units of the inverse damping rate γ −1 [cf.Eq. ( 20)].We consider several values of driving amplitude Ω, including below the exceptional point (green lines), exactly at the exceptional point (orange line) and above the exceptional point (red line), as described in the legend found in panel (b).Panel (b): the degree of second-order coherence g (2) (τ ) as a function of τ [cf.Eq. ( 23)].Dashed grey line: g (2) (τ ) = 1 for coherent light as a guide for the eye.In this figure, we consider the case of the detuning ∆ = 3γ/2, so that ΩEP = 3γ/2 [cf.Eq. ( 10)] and Ωc = 5/2γ 1.58γ [cf.Eq. ( 16)], a case which is given by the blue lines in both panels.
damped-quadratic form g EP We plot g (1) (τ ) as a function of the delay time τ in Fig. 3 (a), which in general shows partial coherence with 0 < |g (1) (τ )| < 1, as opposed to complete coherence with |g (1) (τ )| = 1 or full incoherence with |g (1) (τ )| = 0. We plot the results for driving amplitudes below the EP (green lines), at the EP (orange line), above the EP (red line) and at the critical driving amplitude (blue line).Notably, below the EP (two thickest lines) the correlation function displays characteristic damped oscillations in ωτ , and the exponential damping ensures lim τ →∞ g (1) (τ ) = 0, and thus total incoherence, for all uncritical cases.Above the EP (red line) the behavior becomes, quite predictably, nonoscillatory.However, exactly at the critical driving amplitude Ω = Ω c , Eq. ( 20) reduces to g c (τ ) = 1 and complete coherence (thin blue line) is displayed for all times τ .Temporal light intensity correlations may be measured using the second order correlation function g (2) (τ ), which can likewise be defined in its normalized form as g (2)  2 , for some delay time τ ≥ 0. Bunched light emissions arise from the oscillator when g (2) (0) > g (2) (τ ), while photon antibunching occurs for g (2) (0) < g (2) (τ ).The analytic expression for g (2) (τ ) may be calculated using the quantum regression formula and the first and second moments as (see the Supplementary Information for details) The expression valid for Ω < Ω EP exhibits exponential decay with the time constant γ, such that lim τ →∞ g (2) (τ ) = 1, suggesting a tendency towards Poissonian statistics for large delay times τ .However, for the case of Ω > Ω EP the long delay time asymptotics may be described by which describes a convergent solution only when γ > 2Γ, or equivalently when Ω < Ω c [cf. Eq. ( 16)], which corresponds to the situation when the system is able to support a steady state.At zero time delay (τ = 0), we find from Eq. ( 23) the much simpler expression In this figure, we consider the case of the detuning ∆ = 3γ/2, so that ΩEP = 3γ/2 [cf.Eq. ( 10)] and Ωc = 5/2γ 1.58γ [cf.Eq. ( 16)].
which has the upper bound g (2) (0) (Ω c /Ω) 2 for weak driving amplitudes Ω Ω c , while the lower bound g (2) (0) = 3 for strong drivings tending towards the critical amplitude Ω → Ω c .These results effectively assume zero temperature of the thermal bath encasing the oscillator [cf. the discussion after the quantum master equation of Eq. ( 4)], however the order of limits is important in this case, and the expression of Eq. ( 25) should be replaced in the finite temperature case with a more complicated formula (as given in the Supplementary Information), such that one instead obtains the weak driving limit result lim Ω→0 g (2) (0) = 2 for all temperatures, corresponding to photon bunching.
Exactly at the EP, the following specific form of the correlation function g (2) (τ ) arises [cf.Eq. ( 23) for the marginal cases directly below and above it] where again the EP has reconstructed the response of the system into a damped algebraic one.In Fig. 3 (b), we plot g (2) (τ ) as a function of the delay time τ , using the same colour coding as in panel (a).For drivings below the EP (green lines) we see gentle damped oscillations in the correlation function -always satisfying the bunching inequality g (2) (0) > g (2) (τ ).Exactly at the EP (orange line) we see the damped-quadratic scaling in γτ following Eq.( 26) before a fast washing out of the correlations with large delay times τ → ∞, and a similar pattern occurs for cases above the EP (red line).Finally, at the critical driving amplitude Ω = Ω c (blue line), where Eq. ( 23) reduces to the critical result g c (τ ) = 3, the delay time-independent result is observed.All of these results have been obtained at zero temperature (see the Supplementary Information for a discussion of nonzero temperature).

Optical spectrum.
The optical spectrum S(ω), describing the mean number of photons emitted from the oscillator with the frequency ω, is defined via the integral of the population correlator b [55,56,59].We have included a normalization prefactor in this definition so that the spectral integral is equal to unity, that is ∞ −∞ S(ω)dω = 1.We compute the analytic expression for the spectrum S(ω) of the parametric oscillator as (see the Supplementary Information for the derivation) where the key quantities ω and Γ were introduced in Eq. ( 7) and Eq. ( 8) respectively.This spectral result S(ω) is plotted in Fig. 4 as the thick lines.In panel (a), the driving amplitude Ω is below the EP (Ω < Ω EP ), such that the optical spectrum presents a Eq. ( 10)].Panel (d): Ω = 5/2γ 1.58γ, corresponding to the critical point Ωc [cf.Eq. ( 16)].Panel (e): the quadrature variances σ 2 X (thick green line) and σ 2 P (thick orange line) in the steady state, as a function of Ω.The product of the standard deviations σX σP (medium red line) is also shown, as is a guide for the eye at the Robertson-Schrödinger minimum uncertainty of 1/2 (horizontal, dashed grey line).Vertical grey lines: driving amplitudes corresponding to ΩEP and Ωc respectively.Throughout this figure, we consider the case of the detuning ∆ = 3γ/2.characteristic doublet lineshape (thick green line).This response may have been anticipated due to the presence of two distinct real parts of the complex eigenfrequencies ω ± in this regime [cf.Fig. 1 (b)].However, for the case of Fig. 4 (b), where the driving amplitude is above the EP (Ω > Ω EP ), the spectrum displays a singlet structure (thick red line) since there is now only one distinct real part of ω ± , which is essentially the only transition frequency appearing in the system [cf.Fig. 1 (b)].In both panels of Fig. 4, alongside the full spectrum S(ω) we plot its decomposition into two parts S + (ω) (thin cyan lines) and S − (ω) (thin orange lines), where each contribution S ± (ω) corresponds to an allowed transition in the system (see the Supplementary Information for further information).In panel (a) we therefore see one peak centered around ω and the other peak at −ω, while in panel (b) both peaks are at resonance.The presence of a singlet lineshape within the lower expression within Eq. ( 27) can be more easily seen exactly at the EP, where the optical spectrum S(ω) reduces to the more compact result This expression is of the form of a probability density function for a certain type of Student's t-distribution with the degree of freedom ν = 3 (notably, a Lorentzian distribution corresponds to ν = 2, and a Gaussian distribution arises in the limit of ν → ∞) [60].Hence the optical spectrum S(ω) presents one of the most tangible indicators of passing through an EP, due the stark difference in spectral features presented, namely a doublet to singlet transition.

Quantum states.
The quantum states of the system can be analysed within a phase-space formalism [61,62].In particular, the expectation value of the density operator ρ, with respect to the coherent state |α , is the so-called Husimi function Q (α) = α|ρ|α /π.This Q-function acts like a kind of phase-space quasiprobability distribution, since it has the normalization Q (α) d 2 α = 1 and it is non-negative for all quantum states.In Fig. 5 we show the evolution of the Q-function with increasing driving amplitude Ω across the row of panels (in the steady state t → ∞, where the detuning ∆ = 3γ/2 and for a truncated oscillator with N = 40 levels).In panel (a) the driving amplitude Ω = γ/2, which is sufficiently small such that the Q-function is approximately circular since the quantum state is essentially brought into the vacuum state with n = 0 due to the dissipation (approximately, here Q (α) e −|α| 2 /π and its maximum max{Q} = 1/π 0.318).With increasing drive amplitude up to Ω = γ in panel (b), the Q-function becomes increasingly oval-shaped in the first manifestations of the parametric driving inducing squeezing behavior.In panel (c) the situation at the EP is shown (Ω = Ω EP = 3γ/2), and the significance of the squeezing Hamiltonian of Eq. ( 1) is clearly seen through the highly squeezed Q-function (the tilt is caused by the competition between the driving and the dissipation).Finally in panel (d), where the critical point is finally reached (Ω = Ω c = 5/2γ 1.58γ), the Q-function becomes stadium-shaped (and rather diluted) due to the finite truncation (see the Supplementary Information for more details).While the results of Fig. 5 (a-c) are essentially those for an oscillator with an infinite number of levels, the truncation to N = 40 levels is necessary in panel (d) since the driving amplitude there coincides with the critical point in the untruncated limit [where no steady state exists, as suggested by Eq. ( 15)].
The quasiprobability distribution results of Fig. 5 (a-d) can be better understood in conjunction with an analysis of the variances σ 2 X and σ 2 P of the generalized quadrature operators X and P [cf. the discussion above Eq.( 4)].The two quadrature variances σ 2 X = X2 − X 2 and σ 2 P = P 2 − P 2 should then satisfy, in their standard deviation forms σ X and σ P , the Robertson-Schrödinger uncertainty relation σ X σ P ≥ 1/2.In the steady state (t → ∞), we find the following simple expressions for the variances (see the Supplementary Information for the derivation) where the critical driving amplitude Ω c is given by Eq. ( 16), and where Eq. ( 29) is defined for driving amplitudes Ω < Ω c so that a steady state limit exists.We plot the variances σ 2 X (thick green line) and σ 2 P (thick orange line), as a function of Ω in Fig. 5 (e).Most notably, for drivings Ω < Ω EP the position variance σ 2 X < 1/2, showcasing the same steady-state squeezing as displayed pictorially in Fig. 5 (a-c).In particular, Eq. ( 29) implies that the position variance minimum, min{σ 2 X } = [1 + γ/(2Ω c )]/4, occurs when Ω = Ω c (Ω c − γ/2)/∆.For the case of Fig. 5 (e), where the detuning parameter ∆ = 3γ/2, this means that the extremum min{σ 2 X } = (1 + 1/ √ 10)/4 0.329 is realized when the driving amplitude Ω = γ(10 − √ 10)/6 1.14γ, which is the case somewhere between those described in panels (b) and (c).Notably, exactly at Ω = Ω EP the position variance σ 2 X = 1/2, suggesting the role of the EP as the threshold above which steady state squeezing is absent (this coincidence is seemingly accidental and does not hold at sufficiently high temperatures of the thermal bath, as is shown in the Supplementary Information).Finally, the product of the standard deviations σ X σ P (medium red line) has the lower bound 1/2 in the vanishing driving limit Ω → 0 and tends to infinity in the maximum driving limit of Ω → Ω c , marking the onset of the dynamical instability associated with the loss of the steady state.

Liouvillian gap.
The Liouvillian eigenmatrix L, corresponding to the model of Eq. ( 4) truncated into its finite-dimensional matrix form ∂ t ρ = Lρ, possesses a certain number of eigenvalues (see the Supplementary Information for more details).From this, one may obtain the Liouvillian spectral gap (formally, the Liouvillian eigenvalue with the smallest real part after discounting any zero eigenvalues) [63][64][65].A closing of the Liouvillian gap at some critical value of a system parameter signifies a dissipative phase transition, examples of which have recently been observed in several photonic platforms [66][67][68].We plot the Liouvillian gap as a function of the drive amplitude Ω for the studied squeezed oscillator in Fig. 6 (a).We show results for an increasingly large number of levels N of the oscillator (thin cyan-blue lines), such that in the untruncated limit -that is, the infinite limit N → ∞ as is written in Eq. ( 1) -the Liouvillian gap closes.This closing of the gap occurs at the critical amplitude Ω c (vertical grey line), as was previously defined in Eq. ( 16), and which signifies a dissipative phase transition.Furthermore, the purity P = Tr(ρ) 2 of the density matrix ρ in the steady state suggests that while below Ω c the state is relatively pure (P 1), sharply above Ω c the state is maximally mixed (P 1/N ), as is shown in the Supplementary Information.In the panel directly below, Fig. 6 (d), we plot the average population of the oscillator in the steady state, lim t→∞ b † b , again as a function of Ω and for different level numbers N (thin yellow-red lines) due to the truncation.These restricted results approach the behaviour of the infinite-level oscillator result [thick grey line, Eq. ( 15)] with higher N , revealing that the critical amplitude Ω c is indeed linked to the disappearance of a steady state population (or at least a saturation for finite N cases) as well as a dissipative phase transition.
A note on interactions.Anharmonic deviations to the harmonic oscillator described in Eq. ( 1) can be modelled by adding a Kerr-like energy to the Hamiltonian, Ĥ → Ĥ + ĤI , with the help of the nonlinear term where the on-site interaction strength U ≥ 0 is considered to be repulsive.Carrying out an analysis of the Liouvillian gap of this anharmonic model, again for a spread of different truncations to N energy levels (pink-purple lines), leads to the results of Fig. 6 (b).Notably, unlike the harmonic case of panel (a), the Liouvillian gap in panel (b) seems to decrease with increasing N without a convergence to some critical value of the driving amplitude Ω in the large N thermodynamic limit.This statement is supported by Fig. 6 (c, f).In panel (c), where the minimum of the Liouvillian gap is shown as a function of the restriction parameter N , for the anharmonic case (pink circles) and harmonic case (cyan circles), we see how the Liouvillian gap is closing as a power function and exponentially respectively (coloured lines) via the presented fittings.Meanwhile in panel (f), which displays the value of the drive amplitude Ω at which min {Liouvillian gap} occurs as a function of N , the fittings suggest that the harmonic case (cyan line) approaches Ω c via a power function, while the anharmonic case (pink line) seems to diverge exponentially.This apparent absence of a dissipative phase transition at some critical system parameter for the anharmonic oscillator comes with the guarantee of a steady state population in this case.Indeed, a semiclassical analysis [69] of this nonlinear model [cf.Eq. ( 1) with Eq. ( 30)] leads to the following expression for the steady state population n I (∞) of the anharmonic oscillator (see the Supplementary Information for the derivation) In this figure, we consider the case of the detuning ∆ = γ, so that Ωc = √ 5γ/2 1.12γ [cf.Eq. ( 16)] for the harmonic case.
where the critical frequency Ω c is defined in Eq. ( 16) (cf.Eq. ( 15) for the corresponding harmonic oscillator population result).Clearly, it is necessary (within this semiclassical analysis) for the driving amplitude to be sufficiently strong Ω > Ω c , in order to ensure a physical steady state with n (∞) > 0. The population of Eq. ( 31) is plotted in Fig. 6 (e) as the thick grey line, along with results for various truncations N of the anharmonic oscillator (thin green lines).This panel suggests that Eq. ( 31) is indeed a reasonable approximation in the thermodynamic limit N → ∞ of a large number of excitations.We conclude that while interactions in the form of Eq. ( 30), at least for perturbative values of U , do not strongly modify the majority of our results, this anharmonicity is important for the extinction of the critical closing of the Liouvillian gap and hence for the absence of a dissipative phase transition.

Discussion
In conclusion, we have studied one of the simplest driven-dissipative quantum models in which exceptional points can arisethat of a quantum harmonic oscillator with parametric driving.We have revealed that the exceptional point is of second-order in the first moments of the system, which impacts upon the lineshape of the optical spectrum and the character of the first-order degree of coherence.In the second moments of the system the exceptional point instead appears at the third-order, which influences both the dynamics of the populations and the behaviour of the second-order degree of coherence.Furthermore, the exceptional point is shown to coincide with last remnants of quantum squeezing in the steady state, which perhaps highlights the importance of this kind of physics for quantum states, as described within an open quantum systems approach.We have also discussed the occurrence of a critical point for the parametric oscillator which is associated with a dynamical instability, the phenomenon of a dissipative phase transition (as adjudicated by the closing of the Liouvillian gap), and the impact of small anharmonicities and truncations of the oscillator.Overall, we hope that our results can stimulate further experimental work, with platforms including photonic cavities [70,71] and quantum circuits [72,73], which can observe the predicted non-Hermitian quantum physics through the insightful lens of exceptional points and critical points.

Methods
We use methods from theoretical quantum optics as outlined in the main text and as discussed in detail in the Supplementary Information.
b † and b act on a number state |n as follows: b † |n = √ n + 1|n and b|n = √ n|n−1 , and they obey the bosonic commutation relation [b, b † ] = 1.The diagonalization of the Hamiltonian Ĥ of Eq. (

FIG. 5 .
FIG. 5.The phase-space quasiprobability distribution and variances of the parametric driven-dissipative oscillator.Panels (a-d): the Husimi function Q (α) = α|ρ|α /π, or the expectation value of the density operator ρ with respect to the coherent state |α , for several values of the driving amplitude Ω.The Q-function is calculated in the steady state (t → ∞) for a truncated oscillator with N = 40 levels.Panel (a): Ω = γ/2, where γ is the loss rate.Panel (b): Ω = γ.Panel (c): Ω = 3γ/2, corresponding to the exceptional point ΩEP [cf.Eq. (10)].Panel (d): Ω = 5/2γ 1.58γ, corresponding to the critical point Ωc [cf.Eq. (16)].Panel (e): the quadrature variances σ 2 X (thick green line) and σ 2 P (thick orange line) in the steady state, as a function of Ω.The product of the standard deviations σX σP (medium red line) is also shown, as is a guide for the eye at the Robertson-Schrödinger minimum uncertainty of 1/2 (horizontal, dashed grey line).Vertical grey lines: driving amplitudes corresponding to ΩEP and Ωc respectively.Throughout this figure, we consider the case of the detuning ∆ = 3γ/2.

Ωc= √ 5γ/ 2 FIG. 6 .
FIG.6.The effect of interactions on the parametric driven-dissipative oscillator.Panel (a): the Liouvillian gap as a function of the drive amplitude Ω (in units of the loss rate γ) for the harmonic case (U = 0).Results are shown for various cases when the oscillator is truncated to have a finite number of levels N .Vertical line: the gap closes as Ω → Ωc with N → ∞ [cf.Eq. (16)].Panel (b): as for panel (a), but for a typical anharmonic case (U = γ/100).Panel (c): the minimum value of the Liouvillian gap as a function of the system size N .The algebraic scaling for the harmonic case (cyan line) is γA/N B , with {A, B} = {1.505,0.933}.The anharmonic case (pink line) has the exponential scaling γA/B N , with {A, B} = {1941, 1.230}.Panel (d): the mean steady state population of the oscillator limt→∞ b † b as a function of Ω, computed for the harmonic case (U = 0).Thick grey line: the untruncated oscillator expression of Eq.(15).Panel (e): as for panel (d), but for the considered anharmonic case (U = γ/100).Thick grey line: the semiclassical expression of Eq.(31).Panel (f): the value of Ω/γ at which min {Liouvillian gap} occurs as a function of N .The harmonic case (cyan line) approaches Ωc via the algebraic scaling Ω = Ωc + γA/N B , with {A, B} = {5.282,1.333}.The anharmonic case (pink line) diverges exponentially like Ω = γAB N , with {A, B} = {1.189,1.0024}.In this figure, we consider the case of the detuning ∆ = γ, so that Ωc = √ 5γ/2 1.12γ [cf.Eq. (16)] for the harmonic case.