*Article* **Quantum Behavior of a** *PT* **-Symmetric Two-Mode System with Cross-Kerr Nonlinearity**

**Jan Peˇrina Jr. 1,\* and Antonín Lukš <sup>2</sup>**


Received: 12 July 2019 ; Accepted: 31 July 2019; Published: 7 August 2019

**Abstract:** Quantum behavior of two oscillator modes, with mutually balanced gain and loss and coupled via linear coupling (including energy conserving as well as energy non-conserving terms) and nonlinear cross-Kerr effect, is investigated. Stationary states are found and their stability analysis is given. Exceptional points for PT -symmetric cases are identified. Quantum dynamics treated by the model of linear operator corrections to a classical solution reveals nonclassical properties of individual modes (squeezing) as well as their entanglement.

**Keywords:** nonlinearly coupled oscillators; PT symmetry; cross-Kerr nonlinearity; stability analysis; quantum properties

#### **1. Introduction**

PT -symmetric systems, which contain gain and loss in mutual balance, have been extensively analyzed in various configurations and from different points of view since the pioneering work by Bender and Boettcher occurred [1–3]. The simplest system is composed of two linearly coupled oscillator modes, one exhibiting gain and the other loss [4]. In real physical applications, there occur additional nonlinear Kerr-type terms in both oscillator modes. They originate in physical models of mode amplification and attenuation typically realized via two-level atoms [5]. These models were developed and extensively discussed in the semiclassical and quantum theories of lasers [6]. Stationary states then occur in such systems due to this nonlinearity. The simplest model of two coupled nonlinear oscillator modes has been generalized to include more oscillators in various configurations. The obtained models were applied in many areas of physics including optical coupled structures [7–10], optical waveguides [11,12], coupled optical micro-resonators [13–15], optical lattices [16–19], opto-mechanical systems [20,21], etc.

Recently, attention has been devoted to the consistent quantum description of PT -symmetric systems. To guarantee the validity of commutation relations among the field operators during the evolution, the fluctuating quantum Langevin forces with specific properties have to be considered in the system [22–25]. As a consequence, the noise in the system constantly increases during the evolution both owing to the amplification and attenuation [24]. Despite this, quantum PT -symmetric systems exhibit interesting and appealing features, such as enhancement of interactions around and at exceptional points (EPs) [26] or quantum Zeno effect [27]. The enhancement of nonlinear interactions then opens the door for the generation of nonclassical light (squeezing) and entangled states [28–31].

Here, we investigate the behavior of a specific form of the model of two coupled oscillator modes with amplification and attenuation that includes only the cross-Kerr nonlinear term. In addition, linear coupling of both modes through *χ*(2) parametric interaction that does not conserve energy, is considered [32,33]. It adds or removes photons simultaneously in both modes. This coupling leads to

nonclassical properties of the fields and the occurrence of entanglement between the modes [22,32,34], together with the cross-Kerr nonlinear coupling [35,36]. The cross-Kerr coupling is known to play a significant role in quantum non-demolition measurement [37], generation of the states defined in finite-dimensional Hilbert spaces [38] and generation of maximally entangled Bell-type states [39,40]). In general, the cross-Kerr coupling appearing in the so-called Kerr couplers considerably changes their quantum properties [41–43]. The cross-Kerr nonlinearity can even enhance the usual Kerr effect, e.g., when squeezing effects are analyzed [44].

We show that continuous sets of stationary states occur in the model and we analyze their stability. Then, in the framework of the model of quantum superposition of signal and noise [22], we address squeezed-state generation and generation of entangled states around the stationary states. We note that if one of the oscillator modes in the analyzed model attains an additional Kerr nonlinear term, only the trivial stationary states exist. On the other hand, if the standard Kerr nonlinear terms are attributed to both oscillator modes, the system behavior considerably changes and only discrete stationary states are found [45]. We note that related systems were analyzed from the point of view of squeezed-state generation in [29] (without parametric interaction) and [31] (no Kerr terms, parametric interaction in individual modes) and quantum-noise generation [25] (without parametric interaction). In addition, the work [46], where breaking of the oscillatory regime in a classical two-mode PT -symmetric system with the Kerr nonlinearity due to larger modes intensities is reported, is worth mentioning.

The paper is organized as follows. In Section 2, the analyzed system is defined and the corresponding Heisenberg equations are given. Stationary states and their stability are investigated in Section 3. Nonclassical properties of the evolving states are discussed in Section 4. Section 5 presents conclusions.

#### **2. Quantum Hamiltonian and Dynamical Equations**

Introducing annihilation (*a*ˆ*j*) and creation (*a*ˆ† *<sup>j</sup>* ) operators of photons for oscillator modes 1 (*j* = 1) and 2 (*j* = 2) with identical frequencies *ω*, the considered system is described by the following interaction Hamiltonian *H*ˆ [22]:

$$\hat{H} = -i\gamma\_1 \mathfrak{d}\_1^\dagger \mathfrak{d}\_1 - i\gamma\_2 \mathfrak{d}\_2^\dagger \mathfrak{d}\_2 + \left[\epsilon \mathfrak{d}\_1^\dagger \mathfrak{d}\_2 + \kappa \mathfrak{d}\_1 \mathfrak{d}\_2 + \text{h.c.}\right] + \beta\_\mathfrak{c} \mathfrak{d}\_1^\dagger \mathfrak{d}\_2^\dagger \mathfrak{d}\_1 \mathfrak{d}\_2. \tag{1}$$

We assume that mode 1 is attenuated *γ*<sup>1</sup> ≥ 0 and mode 2 is amplified *γ*<sup>2</sup> ≤ 0. Transfer of energy in the system is described by linear coupling constants  and *κ*. Whereas the coupling constant *-* characterizes transfer of energy between the modes, the constant *κ* quantifies energy inserted and removed to/from both modes in the same amount in the *χ*(2) parametric process. The nonlinear coupling constant *β<sup>c</sup>* characterizes in Equation (1) the cross-Kerr nonlinear term that is responsible for the occurrence of stationary states. Both *χ*(2) term of parametric interaction and cross-Kerr term occur together in nonlinear photonic structures [33] (waveguides, nonlinear fibers). Symbol h.c. replaces the Hermitian conjugated terms.

The Hamiltonian *<sup>H</sup>*<sup>ˆ</sup> in Equation (1) attains its PT -symmetric form provided that the constants *<sup>γ</sup>*1, *γ*2, *-*, *κ* and *β<sup>c</sup>* are real and

$$
\gamma\_1 = -\gamma\_2 \equiv \gamma \ge 0. \tag{2}
$$

Moreover, to allow for simple physical interpretation, PT -symmetric Hamiltonians are usually applied for the range of parameters in which their linear parts are endowed with real eigenvalues. For the Hamiltonian *H*ˆ in Equation (1), this occurs provided that

$$
\kappa^2 - \kappa^2 - \gamma^2 \ge 0. \tag{3}
$$

Points in the space of parameters for which equality in Equation (3) holds identify systems with specific properties. They are called exceptional points (EPs) [1]. Without the loss of generality, we further assume *-*> 0 and *κ* ≥ 0.

Applying the canonical commutation relations for field operators [5] we obtain the Heisenberg equations from the Hamiltonian *H*ˆ in Equation (1),

$$\begin{aligned} \frac{d\hat{\mathfrak{a}}\_1}{dt} &= -\gamma\_1 \mathfrak{a}\_1 - i\varepsilon \mathfrak{a}\_2 - i\kappa \mathfrak{a}\_2^\dagger - i\beta\_c \mathfrak{a}\_2^\dagger \mathfrak{a}\_2 \mathfrak{a}\_1 + \hat{l}\_1, \\ \frac{d\hat{\mathfrak{a}}\_2}{dt} &= -\gamma\_2 \mathfrak{a}\_2 - i\varepsilon \mathfrak{a}\_1 - i\kappa \mathfrak{a}\_1^\dagger - i\beta\_c \mathfrak{a}\_1^\dagger \mathfrak{a}\_1 \mathfrak{a}\_2 + \hat{l}\_2. \end{aligned} \tag{4}$$

and the Hermitian-conjugated ones. The fluctuating Langevin operator forces ˆ *l*<sup>1</sup> and ˆ *l*<sup>2</sup> are introduced in Equation (4) in relation to attenuation in mode 1 and amplification in mode 2, respectively. Their properties [22,23,25],

$$
\begin{aligned}
\langle \hat{l}\_1^\dagger(t) \hat{l}\_1(t') \rangle &= 0, & \langle \hat{l}\_1(t) \hat{l}\_1^\dagger(t') \rangle &= 2\gamma\_1 \delta(t - t'), \\
\langle \hat{l}\_2^\dagger(t) \hat{l}\_2(t') \rangle &= -2\gamma\_2 \delta(t - t'), & \langle \hat{l}\_2(t) \hat{l}\_2^\dagger(t') \rangle &= 0,
\end{aligned}
\tag{5}$$

guarantee validity of the commutation relations for the field operators during the evolution. In mode 1, they express the fluctuation-dissipation theorem [47,48] according to which any dissipation of the energy from a system has to be accompanied by back-action from the environment. In analogy, in mode 2, the Langevin forces represent a part of the 'fluctuation-amplification theorem' that occurs as a consequence of consistent adding the energy into the system [6]. Symbol *δ* means the Dirac function.

#### **3. Stationary States and Their Stability**

First, we address the Heisenberg equations in Equation (4) in their 'classical' noiseless limit, i.e., we write them for the coherent states |*αj* with complex amplitudes *α<sup>j</sup>* = *<sup>j</sup>* exp(*iϕj*), *j* = 1, 2:

$$\frac{d\varrho\_1}{dt} = -\gamma\_1 \varrho\_1 + \left[\epsilon \sin(\varrho) - \kappa \sin(\psi)\right] \varrho\_2. \tag{6}$$

$$\frac{d\varrho\_2}{dt} = -\gamma\_2 \varrho\_2 - \left[\epsilon \sin(\varrho) + \kappa \sin(\psi)\right] \varrho\_1. \tag{7}$$

$$\frac{d\varrho\_1}{dt} = -\left[\epsilon\cos(\varrho) + \kappa\cos(\psi)\right]\frac{\varrho\_2}{\varrho\_1} - \beta\_\varepsilon \varrho\_{2'}^2 \tag{8}$$

$$\frac{d\varrho\_2}{dt} = -\left[\epsilon\cos(\varrho) + \kappa\cos(\psi)\right]\frac{\varrho\_1}{\varrho\_2} - \beta\_\varepsilon\varrho\_1^2. \tag{9}$$

In Equations (6) and (7), we suitably replace the phases *ϕ*<sup>1</sup> and *ϕ*<sup>2</sup> by their sum *ψ* = *ϕ*<sup>2</sup> + *ϕ*<sup>1</sup> and difference *ϕ* = *ϕ*<sup>2</sup> − *ϕ*1.

To reveal the stationary complex amplitudes *α*<sup>1</sup> and *α*2, we set the time derivatives in Equations (6)–(9) to zero. Before analyzing the obtained equations in detail, we note that there exist the trivial stationary states with st <sup>1</sup> = st <sup>2</sup> = 0 and arbitrary values of phases *<sup>ϕ</sup>*st <sup>1</sup> and *<sup>ϕ</sup>*st 2 . Under the stationary conditions, Equations (8) and (9) for the phases *ϕ*<sup>1</sup> and *ϕ*<sup>2</sup> are dependent and, e.g., Equation (8) gives us:

$$
\varrho\_1^{\rm st} \varrho\_2^{\rm st} = - (\mathbf{c}\_{\rm c} + \mathbf{c}\_{\rm \kappa}) / \beta\_{\rm c}. \tag{10}
$$

To simplify the notation, we use in Equation (10) and below the following functions that substitute the phases *ψ* and *ϕ*:

$$\mathbf{s}\_{\mathbf{c}} = \boldsymbol{\varepsilon}\sin(\boldsymbol{\varphi}^{\mathrm{st}}), \ \mathbf{c}\_{\mathbf{c}} = \boldsymbol{\varepsilon}\cos(\boldsymbol{\varphi}^{\mathrm{st}}), \ \mathbf{s}\_{\mathbf{x}} = \boldsymbol{\kappa}\sin(\boldsymbol{\psi}^{\mathrm{st}}), \ \mathbf{c}\_{\mathbf{x}} = \boldsymbol{\kappa}\cos(\boldsymbol{\psi}^{\mathrm{st}}).\tag{11}$$

Furthermore, the coupled Equations (6) and (7) considered to be functions of st <sup>1</sup> and st <sup>2</sup> have a nontrivial solution provided that their determinant is zero:

$$
\mathbf{s}\_{\mathbf{c}}^2 - \mathbf{s}\_{\mathbf{c}}^2 + \gamma\_1 \gamma\_2 = 0. \tag{12}
$$

Equation (6) when combined with Equation (10) gives us the stationary solution for amplitudes:

$$
\rho\_{1,2}^{\rm st} = \sqrt{\frac{\gamma\_{2,1}}{\beta\_{\rm c}} \frac{\mathbf{c\_{c}} + \mathbf{c\_{x}}}{\pm \mathbf{s\_{c}} + \mathbf{s\_{x}}}}.\tag{13}
$$

According to Equations (12) and (13), one phase, e.g., *ψ*st in the stationary solution is free. The other phase, *ϕ*st, has to fulfill Equation (12) that admits in general four solutions. The amplitudes st 1,2 determined by Equation (13) have to be real and also Equation (10) has to give nonnegative st <sup>1</sup> st 2 .

To address stability of the stationary solution, we derive the linearized equations for deviations *δ*1, *δ*2, *δϕ* and *δψ* from their corresponding stationary values st <sup>1</sup> , st <sup>2</sup> , *<sup>ϕ</sup>*st, and *<sup>ψ</sup>*st:

$$
\frac{d}{dt} \begin{bmatrix} \delta \varrho\_1 \\ \delta \varrho\_2 \\ \delta \varrho \\ \delta \Psi \end{bmatrix} = \begin{bmatrix} -\gamma\_1 & s\_{\mathfrak{c}} - s\_{\mathfrak{k}} & c\_{\mathfrak{c}} \varrho\_2^{\mathfrak{k}} & -c\_{\mathfrak{k}} \varrho\_2^{\mathfrak{k}} \\ -s\_{\mathfrak{c}} - s\_{\mathfrak{k}} & -\gamma\_2 & -c\_{\mathfrak{c}} \varrho\_1^{\mathfrak{k}} & -c\_{\mathfrak{k}} \varrho\_1^{\mathfrak{k}} \\ G^+ & H^+ & I^+ \mathfrak{s}\_{\mathfrak{c}} & I^+ \mathfrak{s}\_{\mathfrak{k}} \\ G^- & H^- & I^- \mathfrak{s}\_{\mathfrak{c}} & I^- \mathfrak{s}\_{\mathfrak{k}} \end{bmatrix} \begin{bmatrix} \delta \varrho\_1 \\ \delta \varrho\_2 \\ \delta \varrho \\ \delta \Psi \end{bmatrix}. \tag{14}
$$

The parameters *G*±, *H*± and *I*± are given as follows:

$$\begin{array}{rcl} G^{\pm} &=& \frac{\beta\_{\mathfrak{c}} \mathfrak{c}\_{1}^{\mathfrak{st}}}{\mathbf{s}\_{\mathfrak{c}} - \mathbf{s}\_{\mathfrak{c}}} \left[ \mathbf{s}\_{\mathfrak{c}} \left( \mp \gamma\_{1} / \gamma\_{2} - 1 \right) + \mathbf{s}\_{\mathfrak{c}} \left( \mp \gamma\_{1} / \gamma\_{2} + 1 \right) \right], \\ H^{\pm} &=& \frac{\beta\_{\mathfrak{c}} \mathfrak{c}\_{2}^{\mathfrak{st}}}{\mathbf{s}\_{\mathfrak{c}} + \mathbf{s}\_{\mathfrak{c}}} \left[ \mathbf{s}\_{\mathfrak{c}} \left( \gamma\_{2} / \gamma\_{1} \pm 1 \right) + \mathbf{s}\_{\mathfrak{c}} \left( -\gamma\_{2} / \gamma\_{1} \pm 1 \right) \right], \\ I^{\pm} &=& \left( \mathbf{s}\_{\mathfrak{c}} - \mathbf{s}\_{\mathfrak{c}} \right) / \gamma\_{1} \pm \left( \mathbf{s}\_{\mathfrak{c}} + \mathbf{s}\_{\mathfrak{c}} \right) / \gamma\_{2}. \end{array} \tag{15}$$

For the PT -symmetric case, eigenvalues *ν* of the dynamical matrix from Equation (14) are given as roots of the following polynomial:

$$\nu \left( \nu^{\natural} + b\nu + c \right) = 0,$$

$$b = 4(\mathbf{c}\_{\mathbf{c}} + \mathbf{c}\_{\mathbf{k}})(\mathbf{c}\_{\mathbf{c}}\mathbf{s}\_{\mathbf{k}}^2 + \mathbf{c}\_{\mathbf{k}}\mathbf{s}\_{\mathbf{c}}^2)/\gamma^2, \ c = -8\mathbf{s}\_{\mathbf{c}}\mathbf{s}\_{\mathbf{k}}(\mathbf{c}\_{\mathbf{c}} + \mathbf{c}\_{\mathbf{k}})^2/\gamma. \tag{16}$$

The eigenvalue *ν*<sup>1</sup> = 0 is related to the freedom in determining, e.g., the phase *ψ*st of a stationary state. Provided that c*-* + c*<sup>κ</sup>* = 0, we have *ν*1−<sup>4</sup> = 0 and Equation (13) gives us the trivial stationary state st <sup>1</sup> = st <sup>2</sup> = 0 lying on the border of stability.

Assuming PT -symmetry and special case without *<sup>χ</sup>*(2) interaction (*<sup>κ</sup>* <sup>=</sup> 0), we have s*<sup>κ</sup>* <sup>=</sup> <sup>c</sup>*<sup>κ</sup>* <sup>=</sup> <sup>0</sup> and the phase *ψ*st is arbitrary. On the other hand, one solution for the remaining parameters of the stationary state is derived from Equations (12) and (13) for *β<sup>c</sup>* > 0 in the following implicit form:

$$\mathbf{c}\_{\mathfrak{c}} = \gamma,\ \mathbf{c}\_{\mathfrak{c}} = -\sqrt{\epsilon^2 - \gamma^2},\ \ q\_1^{\mathfrak{st}} = \varrho\_2^{\mathfrak{st}} = \sqrt[4]{\epsilon^2 - \gamma^2} / \sqrt{\beta\_{\mathfrak{c}}}.\tag{17}$$

Similarly, we reveal one stationary solution for *β<sup>c</sup>* < 0:

$$\mathbf{s}\_{\mathfrak{c}} = \gamma, \; \mathbf{c}\_{\mathfrak{c}} = \sqrt{\epsilon^2 - \gamma^2}, \; \varrho\_1^{\mathfrak{st}} = \varrho\_2^{\mathfrak{st}} = \sqrt[4]{\epsilon^2 - \gamma^2} / \sqrt{-\beta\_{\mathfrak{c}}}.\tag{18}$$

Eigenvalues of the dynamical matrix in Equation (14) are obtained for both solutions in Equations (17) and (18) as *ν*1−<sup>4</sup> = 0, i.e., the states are at the border of stability.

For nonzero *κ*, we first address the stationary states in EPs (PT -symmetric case) for specific values of the phase *ψ*st. The trivial stationary states with st <sup>1</sup> = st <sup>2</sup> = 0 and zero eigenvalues *ν*1−<sup>4</sup> = 0 in the stability analysis are found under the conditions summarized in Table 1.

**Table 1.** Parameters of the trivial stationary states st <sup>1</sup> <sup>=</sup> st <sup>2</sup> = 0 with zero eigenvalues *<sup>ν</sup>*1−<sup>4</sup> = 0 for specific values of the phase *ψ*st.


The nontrivial stationary solution with st <sup>1</sup> = st <sup>2</sup> = 2*κ*/*β<sup>c</sup>* in EPs is revealed for *<sup>ψ</sup>*st = *<sup>π</sup>* and *β<sup>c</sup>* > 0:

$$\mathbf{s}\_{\mathbf{k}} = \mathbf{0}, \ \mathbf{c}\_{\mathbf{k}} = -\mathbf{c}, \ \mathbf{s}\_{\mathbf{c}} = \gamma, \ \mathbf{c}\_{\mathbf{c}} = -\mathbf{c}. \tag{19}$$

On the other hand, we have for *β<sup>c</sup>* < 0 and *ψ*st = 0 the stationary solution with st <sup>1</sup> = st <sup>2</sup> = 2*κ*/(−*βc*) in EPs provided that

$$\mathbf{s}\_{\mathbf{x}} = \mathbf{0}, \ \mathbf{c}\_{\mathbf{x}} = \kappa, \ \mathbf{s}\_{\mathbf{c}} = \gamma, \ \mathbf{c}\_{\mathbf{c}} = \kappa. \tag{20}$$

Both solutions in Equations (19) and (20) have the same eigenvalues *ν*1,2 = 0 and *ν*3,4 = ±*i*2 <sup>√</sup>2*<sup>κ</sup>* in the stability analysis, i.e., no amplification of amplitude fluctuations occur.

For an arbitrary phase *ψ*st, Equation (12) admits four possible values for the phase *ϕ*st:

$$q\mathbf{q}\_1^{\rm st} = q\mathbf{q}\_{\rm base}^{\rm st}, \quad q\mathbf{q}\_2^{\rm st} = \boldsymbol{\pi} - q\mathbf{q}\_{\rm base}^{\rm st}, \quad q\mathbf{q}\_3^{\rm st} = \boldsymbol{\pi} + q\mathbf{q}\_{\rm base}^{\rm st}, \quad q\mathbf{q}\_4^{\rm st} = 2\boldsymbol{\pi} - q\mathbf{q}\_{\rm base}^{\rm st} \tag{21}$$

where *ϕ*st base = arcsin( *κ*<sup>2</sup> sin2(*ψ*st) + *γ*2/*-*). However, only some of them lead to real and nonnegative amplitudes st <sup>1</sup> and st <sup>2</sup> in Equation (13) and nonnegative expression in Equation (10). According to Equation (12), |s*-*|≥|s*κ*|. Equation (12) can also be recast into the form suitable for the discussion:

$$
\mathfrak{c}\_{\mathfrak{e}}^{2} - \mathfrak{c}\_{\mathfrak{x}}^{2} = \mathfrak{e}^{2} - \mathfrak{x}^{2} - \gamma^{2}.\tag{22}
$$

Considering Equation (22) for *-*<sup>2</sup> <sup>−</sup> *<sup>κ</sup>*<sup>2</sup> <sup>−</sup> *<sup>γ</sup>*<sup>2</sup> <sup>≥</sup> 0 we have <sup>|</sup>c*-*|≥|c*κ*|. If *β<sup>c</sup>* > 0 [*β<sup>c</sup>* < 0], Equation (10) requires c*-* ≤ 0 [c*-* <sup>≥</sup> 0] and this admits the phases *<sup>ϕ</sup>*st <sup>2</sup> and *<sup>ϕ</sup>*st <sup>3</sup> [*ϕ*st <sup>1</sup> and *<sup>ϕ</sup>*st 4 ]. The expression (13) for st <sup>1</sup> requires s*-* <sup>≥</sup> 0 independently of the sign of *<sup>β</sup><sup>c</sup>* and so only the phases *<sup>ϕ</sup>*st 1 and *ϕ*st <sup>2</sup> can be considered. Both conditions are fulfilled only for the phase *<sup>ϕ</sup>*st <sup>2</sup> [*ϕ*st <sup>1</sup> ] for *β<sup>c</sup>* > 0 [*β<sup>c</sup>* < 0].

On the other hand, we have |c*-*|≤|c*κ*| for *-*<sup>2</sup> <sup>−</sup> *<sup>κ</sup>*<sup>2</sup> <sup>−</sup> *<sup>γ</sup>*<sup>2</sup> <sup>≤</sup> 0. In this case, c*<sup>κ</sup>* <sup>≤</sup> 0 [c*<sup>κ</sup>* <sup>≥</sup> 0] is needed in Equation (10) for *<sup>β</sup><sup>c</sup>* <sup>&</sup>gt; 0 [*β<sup>c</sup>* <sup>&</sup>lt; 0] and so *<sup>ψ</sup>*st <sup>∈</sup>*π*/2, 3*π*/2 [*ψ*st <sup>∈</sup>0, *<sup>π</sup>*/2∪3*π*/2, 2*π*]. Similarly as above, s*-* <sup>≥</sup> 0 guarantees nonnegative expression (13) for st <sup>1</sup> for arbitrary *βc*, i.e., only the phases *ϕ*st <sup>1</sup> and *<sup>ϕ</sup>*st <sup>2</sup> are allowed. According to both conditions, nontrivial stationary states are expected for the phases *ϕ*st <sup>1</sup> and *<sup>ϕ</sup>*st <sup>2</sup> in the interval of phase *<sup>ψ</sup>*st <sup>∈</sup>*π*/2, 3*π*/2 [*ψ*st <sup>∈</sup>0, *<sup>π</sup>*/2∪3*π*/2, 2*π*] for *β<sup>c</sup>* > 0 [*β<sup>c</sup>* < 0].

The EPs occurring at the border of two above discussed regions need special attention. According to Equation (22), we have |c*-*| = |c*κ*| at the EPs. The analysis reveals that on the top of the stationary states characterized in the above two paragraphs, the trivial stationary states with st <sup>1</sup> = st <sup>2</sup> = 0 exist for the phase *ϕ*st <sup>3</sup> in the interval *<sup>ψ</sup>*st <sup>∈</sup>0, *<sup>π</sup>*/2∪3*π*/2, 2*π* and for the phase *<sup>ϕ</sup>*st <sup>4</sup> in the interval *<sup>ψ</sup>*st <sup>∈</sup>*π*/2, 3*π*/2 independently of the sign of *<sup>β</sup>c*.

The above general conclusions are further illustrated in the graphs in Figure 1 where the stationary states and their stability are analyzed in the plane (*γ*/*-*, *ψ*st) for the case *β<sup>c</sup>* > 0. In Figure 1, we characterize the stationary states by intensities st2 <sup>1</sup> and st2 <sup>2</sup> of modes 1 and 2, respectively. We judge the stationary states according to the maximal values of imaginary and real parts of the complex frequencies *ν*1−4. A positive (negative) imaginary part means amplification (attenuation) of amplitude fluctuations around the stationary state (we note the inverse notation for signs for *ν* and *γ*). A nonzero real part then indicates oscillations in the evolution of amplitude fluctuations. We have *κ*/*-*= 0.5 for

the graphs drawn in Figure 1 and so the EP occurs for *γ*EP/*-* <sup>=</sup> <sup>√</sup>3/2 <sup>≈</sup> 0.87. The stationary solutions for the phase *ϕ*st <sup>1</sup> exist only in the area with exponential increase of amplitudes (*γ* ≥ *γ*EP), they are unstable and amplitude fluctuations oscillate. Only at the EP, the amplitudes st <sup>1</sup> and st <sup>2</sup> are zero and the state is at the border of stability. On the other hand, there exist stationary states for the phase *ϕ*st <sup>2</sup> in the oscillatory regime of amplitude evolution (*γ* < *γ*EP). According to the graph in Figure 1b amplitude fluctuations around these stationary states oscillate. They are amplified except for the line *ψ*st = *π* that lies at the border of stability (see Figure 1d). According to the graphs in Figure 1e–h, the pattern of intensity st2 <sup>1</sup> of mode 1 is a mirror image of the pattern of intensity st2 <sup>2</sup> of mode 2 with respect to the plane *<sup>ψ</sup>*st <sup>=</sup> *<sup>π</sup>* (compare Equation (13) for *<sup>γ</sup>*<sup>2</sup> <sup>=</sup> <sup>−</sup>*γ*<sup>1</sup> and <sup>±</sup>s*κ*).

The analysis of the graphs for the case with *β<sup>c</sup>* < 0 drawn under the conditions of the graphs in Figure 1 reveals similarity provided that we replace *ϕ*st <sup>1</sup> by *<sup>ϕ</sup>*st <sup>2</sup> , shift the phase *<sup>ψ</sup>*st by *<sup>π</sup>* and exchange mode amplitudes st <sup>1</sup> and st <sup>2</sup> . This similarity is illustrated in the graphs in Figure 2 where the stationary intensities and stability parameters are drawn at the EP for both cases. Considering *β<sup>c</sup>* > 0 [*β<sup>c</sup>* < 0] and following the graphs in Figure 2, there exist only the trivial stationary states with st <sup>1</sup> = st <sup>2</sup> = 0 at the border of stability for *<sup>ψ</sup>*st <sup>∈</sup>0, *<sup>π</sup>*/2∪3*π*/2, 2*π* [*ψ*st <sup>∈</sup>*π*/2, 3*π*/2] and *<sup>ϕ</sup>*st <sup>2</sup> [*ϕ*st <sup>1</sup> ]. For the remaining phases *ψ*st nontrivial stationary states are found. These states are unstable except for the state with *ψ*st = *π* [*ψ*st = 0] and *ϕ*st <sup>2</sup> [*ϕ*st <sup>1</sup> ] that is at the border of stability and amplitude fluctuations around this state oscillate. Parameters for this state are given in Equation (19) (Equation (20)).

**Figure 1.** Maximal values of the real and imaginary parts of four complex frequencies *ν* in the stability analysis expressed via functions Ω*<sup>ν</sup>* = log[1 + |Re(*ν*)|] (**a**,**b**) and Γ*<sup>ν</sup>* = sign[Im(*ν*)] log[1 + |Im(*ν*)|] (**c**,**d**), respectively, and intensities st2 <sup>1</sup> (**e**,**f**) and st2 <sup>2</sup> (**g**,**h**) of modes 1 and 2, respectively, as they depend on dimensionless attenuation/amplification parameter *γ*/ and phase *ψ*st for stationary states with *ϕ*st <sup>1</sup> (**a**,**c**,**e**,**g**) and *<sup>ϕ</sup>*st <sup>2</sup> (**b**,**d**,**f**,**h**) defined in Equation (21); symbol sign gives the sign of the argument, log means the decimal logarithm, Re (Im) stands for the real (imaginary) part of the argument, *κ*/*-* = 0.5, and *βc*/*-*= 0.1.

**Figure 2.** Maximal values of the real and imaginary parts of four complex frequencies *ν* in the stability analysis expressed by functions <sup>Ω</sup>*<sup>ν</sup>* (**a**,**d**) and <sup>Γ</sup>*<sup>ν</sup>* (**b**,**e**), respectively, and intensities st2 of mode 1 (∗) and 2 () (**c**,**f**) as they depend on phase *<sup>ψ</sup>*st for *<sup>ϕ</sup>*st <sup>2</sup> (see Equation (21)) and *βc*/*-* = 0.1 (**a**–**c**) and *ϕ*st <sup>1</sup> and *βc*/*-* <sup>=</sup> <sup>−</sup>0.1 (**d**–**f**); EP condition *<sup>γ</sup>* <sup>=</sup> <sup>√</sup>*-*<sup>2</sup> − *<sup>κ</sup>*<sup>2</sup> is assumed. Functions <sup>Ω</sup>*<sup>ν</sup>* and <sup>Γ</sup>*<sup>ν</sup>* as well as the other parameters are given in the caption to Figure 1.

#### **4. Quantum Properties of the Evolving States**

We analyze the properties of states evolving from the stationary states determined above and compare them with those characterizing the states originating from non-stationary states. For an initial stationary state, the evolution is described by the Heisenberg equations in Equation (4) linearized around the initial complex amplitudes *α*st <sup>1</sup> and *<sup>α</sup>*st <sup>2</sup> (*a*ˆ*<sup>j</sup>* = *<sup>α</sup>*st *<sup>j</sup>* + *δa*ˆ*j*, *j* = 1, 2),

$$\begin{array}{rcl} \frac{d\delta\hat{a}\_{1}}{dt} &=& -\left(\gamma\_{1} + i\beta\_{c}|a\_{2}^{\mathrm{st}}|^{2}\right)\delta\mathbb{i}\_{1} - i\left(\boldsymbol{\varepsilon} + \beta\_{c}a\_{1}^{\mathrm{st}}a\_{2}^{\mathrm{st}}\right)\delta\mathbb{i}\_{2} - i\left(\boldsymbol{\kappa} + \beta\_{c}a\_{1}^{\mathrm{st}}a\_{2}^{\mathrm{st}}\right)\delta\mathbb{i}\_{2}^{\dagger} + \hat{l}\_{1} \\ \frac{d\delta\hat{a}\_{2}}{dt} &=& -i\left(\boldsymbol{\varepsilon} + \beta\_{c}a\_{1}^{\mathrm{st}}a\_{2}^{\mathrm{st}}\right)\delta\mathbb{i}\_{1} - i\left(\boldsymbol{\kappa} + \beta\_{c}a\_{1}^{\mathrm{st}}a\_{2}^{\mathrm{st}}\right)\delta\mathbb{i}\_{1}^{\dagger} - \left(\gamma\_{2} + i\beta\_{c}|a\_{1}^{\mathrm{st}}|^{2}\right)\delta\mathbb{i}\_{2} + \hat{l}\_{2}, \end{array} \tag{23}$$

and the Hermitian-conjugated ones. When an initial non-stationary state is assumed, we numerically solve the classical nonlinear Equations (6)–(9) and linearize the Heisenberg equations around the evolving complex amplitudes *α*1(*t*) and *α*2(*t*) [45]. In both cases the solution can be expressed in the following general form

$$
\begin{bmatrix}
\delta\boldsymbol{\delta}\_{1}(t) \\
\delta\boldsymbol{\hat{a}}\_{2}(t) \\
\end{bmatrix} = \mathbf{U}(t) \begin{bmatrix}
\delta\boldsymbol{\delta}\_{1}(0) \\
\delta\boldsymbol{\hat{a}}\_{2}(0) \\
\end{bmatrix} + \mathbf{V}(t) \begin{bmatrix}
\delta\boldsymbol{\delta}\_{1}^{\dagger}(0) \\
\delta\boldsymbol{\hat{a}}\_{2}^{\dagger}(0) \\
\end{bmatrix} + \begin{bmatrix}
\boldsymbol{\hat{f}}\_{1}(t) \\
\boldsymbol{\hat{f}}\_{2}(t) \\
\end{bmatrix},
$$

$$
\begin{bmatrix}
\boldsymbol{\hat{f}}\_{1}(t) \\
\boldsymbol{\hat{f}}\_{2}(t) \\
\end{bmatrix} = \int\_{0}^{t} dt' \mathbf{U}(t - t') \begin{bmatrix}
\boldsymbol{\hat{f}}\_{1}(t') \\
\boldsymbol{\hat{f}}\_{2}(t') \\
\end{bmatrix} + \int\_{0}^{t} dt' \mathbf{V}(t - t') \begin{bmatrix}
\boldsymbol{\hat{f}}\_{1}^{\dagger}(t') \\
\boldsymbol{\hat{f}}\_{2}^{\dagger}(t') \\
\end{bmatrix},
\tag{24}
$$

in which the matrices **U**(*t*) and **V**(*t*) and correlation functions of the fluctuating operator forces ˆ *fj*(*t*), *j* = 1, 2, are determined numerically in general [45] and analytically under specific conditions [35,45].

We consider the initial vacuum state for the evolution of operator amplitude corrections. In this case, the evolving states remain Gaussian and so the following six correlation functions characterize them completely [22,35]:

$$
\langle \delta \mathfrak{d}\_j^\dagger(t) \delta \mathfrak{d}\_j(t) \rangle,\quad \langle [\delta \mathfrak{d}\_j(t)]^2 \rangle,\quad j = 1,2,\quad \langle \delta \mathfrak{d}\_1^\dagger(t) \delta \mathfrak{d}\_2(t) \rangle,\quad \langle \delta \mathfrak{d}\_1(t) \delta \mathfrak{d}\_2(t) \rangle.\tag{25}
$$

They are easily determined from Equation (24). All quantities characterizing the evolving states can then be expressed in terms of the correlation functions (25). For example, the principal squeeze variance of mode *j* is obtained as [49]

$$\lambda\_j = 1 + 2\left[ \langle \delta \hat{a}\_j^\dagger(t) \delta \hat{a}\_j(t) \rangle - |\langle [\delta \hat{a}\_j(t)]^2 \rangle| \right], \quad j = 1, 2. \tag{26}$$

Determination of the covariance matrix in the symmetric operator ordering then allows to reach the logarithmic negativity *EN* that is a suitable quantifier of the entanglement between the modes (for details, see [50,51]).

In Figure 3, we compare the state evolution around a stationary state with that occurring around a non-stationary state. As a stationary state, we consider the state given in Equation (19). The analyzed non-stationary state evolves from the state that differs from that in Equation (19) in the phase *ψ*init = 0:

$$\mathbf{s}\_{\mathbf{x}}^{\text{init}} = 0, \quad \mathbf{c}\_{\mathbf{x}}^{\text{init}} = \mathbf{x}, \; \mathbf{s}\_{\mathbf{c}}^{\text{init}} = -\gamma, \; \mathbf{c}\_{\mathbf{c}}^{\text{init}} = \mathbf{x}. \tag{27}$$

For the stationary state that lies at the border of stability, both intensities <sup>2</sup> <sup>1</sup>(*t*) and <sup>2</sup> <sup>2</sup>(*t*) increase during the evolution and the fluctuating forces give dominant contribution to this increase (compare solid and dashed curves in Figure 3a). Contrary to this, for the initial non-stationary state the intensity 2 <sup>1</sup>(*t*) of attenuated mode 1 first considerably decreases whereas the intensity <sup>2</sup> <sup>2</sup>(*t*) of amplified mode 2 increases constantly. According to the curves in Figure 3d, the relative contribution of fluctuating forces to the dynamics of intensities is small. In both cases, only the attenuated mode 1 exhibits squeezing (for the principal squeeze variances *λ*1,2, see Figure 3b,e) and both modes are entangled (for the logarithmic negativity *EN*, see Figure 3c,f) for a limited time period. It is worth noting that both squeezing and entanglement are stronger for the initial non-stationary state. The comparison of curves in Figure 3b,d for the principal squeeze variances *λ*1,2 and in Figure 3c,e for the logarithmic negativity *EN* drawn with and without the inclusion of fluctuating forces clearly documents substantial role of these forces in consistent description of PT -symmetric quantum systems.

**Figure 3.** Intensities <sup>2</sup> (**a**,**d**) and principal squeeze variances *<sup>λ</sup>* of mode 1 (∗) and 2 () (**b**,**e**) and logarithmic negativity *EN* (**c**,**f**) as they evolve along dimensionless time *t* for initial stationary (non-stationary) state with *ψ*st = *π* [*ψ*init = 0], 1,2 = 2*κ*/*β<sup>c</sup>* and *ϕ*st <sup>2</sup> (**a**–**c**) [*ϕ*init <sup>4</sup> (**c**–**e**)] given in Equation (19) (Equation (27)), assuming *<sup>γ</sup>* <sup>=</sup> <sup>√</sup>*-*<sup>2</sup> − *<sup>κ</sup>*2, *<sup>κ</sup>*/*-* = 0.5, and *βc*/*-* = 0.1. The initial vacuum state is assumed, evolution is treated with [without] fluctuating forces (black solid [red dashed] curves).

#### **5. Conclusions**

Two oscillator modes with balanced attenuation and amplification were considered to be mutually coupled via the usual linear coupling, *χ*(2) parametric process and cross-Kerr nonlinearity. Nontrivial stationary states that occur owing to the cross-Kerr nonlinearity were identified and their stability was determined. The stationary states typically form one-parameter systems. There occur only unstable stationary states and states lying at the border of stability (zero imaginary parts of frequencies in the stability analysis). The solution of linearized operator equations for mode amplitudes around these stationary states revealed nonclassical properties of the evolving states (single-mode squeezing, entanglement). Initial non-stationary states seem to be more suitable for nonclassical-state generation than the stationary ones at the border of stability. Substantial role of the fluctuating Langevin operator forces in consistent description of the system was demonstrated.

**Author Contributions:** Conceptualization, J.P.J. and A.L.; Investigation, J.P.J. and A.L.; Software, J.P.J.; Writing—original draft, J.P.J.; Writing—review & editing, J.P.J. and A.L.

**Funding:** J.P. acknowledges the support by the GA CR project 18-22102S. A.L. gratefully acknowledges the ˇ support from the project IGA\_PrF\_2019\_008 of Palacký University.

**Acknowledgments:** The authors thank J. K. Kalaga and W. Leo ´nski for discussions about PT-symmetry and Kerr effect.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
