Revisiting the Optical PT-Symmetric Dimer

Abstract

Optics has proved a fertile ground for the experimental simulation of quantum mechanics. Most recently, optical realizations of $\mathcal{PT}$-symmetric quantum mechanics have been shown, both theoretically and experimentally, opening the door to international efforts aiming at the design of practical optical devices exploiting this symmetry. Here, we focus on the optical $\mathcal{PT}$-symmetric dimer, a two-waveguide coupler where the materials show symmetric effective gain and loss, and provide a review of the linear and nonlinear optical realizations from a symmetry-based point of view. We go beyond a simple review of the literature and show that the dimer is just the smallest of a class of planar N-waveguide couplers that are the optical realization of the Lorentz group in 2 + 1 dimensions. Furthermore, we provide a formulation to describe light propagation through waveguide couplers described by non-Hermitian mode coupling matrices based on a non-Hermitian generalization of the Ehrenfest theorem.

1. Introduction

The desire to create an optical directional coupler, a device composed of parallel optical waveguides close enough that leaked energy is transferred between them, led to the exploration of waveguide creation in semi-conductors via proton bombardment [1]. At the time, power losses played an interesting role, and the nascent mode coupling theory [2] allowed the theoretical description of linear loses in such devices [3],

$$\begin{array}{c}\hfill -i\frac{d}{dz}{E}_n(z)=\frac{i\alpha }{2}{E}_n(z)+K[E_{n-1}(z)+E_{n+1}(z)]\end{array}$$

(1)

where the real numbers α and K are the effective linear loss, identical in all implanted waveguides, and the effective waveguide coupling strength, also identical for the whole system, in that order. This, to the best of our knowledge, was the first theoretical description of an experimental N-waveguide coupler including losses in the form of a Schrödinger-like equation involving a non-Hermitian Hamiltonian. Almost twenty years later, the desire to create an intensity-dependent switch working at low power levels took another team of researchers to explore twin core nonlinear couplers with gain and loss. Again, coupled mode theory allowed the description of such devices [4],

$$\begin{array}{ccc}\hfill -i\frac{d}{dz}{E}_1(z)& =& \left({\beta }_1+i{\alpha }_1\right){\mathcal{E}}_1+\Delta {\beta }_1\left({|{\mathcal{E}}_1|}^2\right){\mathcal{E}}_1+K{\mathcal{E}}_2\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill -i\frac{d}{dz}{E}_2(z)& =& \left({\beta }_2+i{\alpha }_2\right){\mathcal{E}}_2+\Delta {\beta }_2\left({|{\mathcal{E}}_2|}^2\right){\mathcal{E}}_2+K{\mathcal{E}}_1\hfill \end{array}$$

(3)

where we have kept the notation used before and introduce the effective real part of the refractive index, ${\beta }_j$, and the real function $\Delta {\beta }_j({|{\mathcal{E}}_j|}^2)=\pm {\kappa }_j{|{\mathcal{E}}_j|}^2$ that describes an effective Kerr nonlinearity-induced change in the refractive index of the j-th core, positive for self-focusing and negative for self-defocusing materials. At this point in history, there existed experimental and theoretical work describing an optical dimer where the waveguides present effective loss and gain and a nonlinearity, but a little something was missing. A couple of years later, a theory exploring a particular type of non-Hermitian Hamiltonians with real spectra was brought forward in quantum mechanics [5]. These Hamiltonians were invariant under space-time reflection, received the name of $\mathcal{PT}$-symmetric and opened new avenues of research in quantum mechanics, as well as other areas of physics and mathematics, cf. [6,7] and the references therein.

In optics, it took a few years more to propose two seminal ideas. The first one regarded single elements and showed, in particular, that an optical planar slab waveguide composed of two media with linear gain and loss can be described by a Schrödinger-like equation under dynamics dictated by a $\mathcal{P}\mathcal{T}$-symmetric Hamiltonian, where the optical refractive index played the role of quantum-like potential and propagation distance that of time [8]. The second one dealt with composite systems, where a mature mode coupling theory produced a theory of coupled optical $𝒫$$𝒯$-symmetric structures [9]. In the ten years following those first proposals for an optical realization of $\mathcal{P}\mathcal{T}$-symmetry, work has been reported on slab waveguides [10,11], Bragg scatterers [12,13,14,15,16,17], as well as linear [18,19,20,21,22,23,24,25,26] and nonlinear [27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45] coupled waveguides, to mention just a few. Research in this field is slowly getting to information technologies applications with recent proposals of all-optical $\mathcal{PT}$-symmetric logic gates [46] and amplitude-to-phase converters [47].

Here, we will provide a review of the optical $\mathcal{PT}$-symmetric dimer. First, in Section 2, we will introduce a two-waveguide coupler where component waveguides show effective complex refractive indices with an identical real part. Starting from this device, we will recover the effective mode coupling differential equation set for the linear $\mathcal{PT}$-symmetric dimer, which describes a nonunitary optical device showing symmetric effective loss and gain feasible for passive and active optical realizations. Then, we will recover the dispersion relation for the dimer that shows three regimes, one with real eigenvalues, the $\mathcal{PT}$-symmetric regime, another with fully-degenerate eigenvalues equal to zero, the fully-degenerate regime, and a third one with purely imaginary eigenvalues, the broken symmetry regime. We will construct an analytic propagator that will show asymmetric amplifying oscillator, power amplification and exponential amplification behaviors in each of these regimes. We will also show that it is possible to uncouple the mode coupling differential equation set of the $\mathcal{PT}$-symmetric dimer. The resulting second order differential equations and boundary conditions for the field amplitudes propagating through each waveguides take the form of nonlocal oscillators with positive potential, a free particle traveling through a nonlocal medium and nonlocal oscillators with inverted potential, in each of the regimes. In the final part of Section 2, we will bring forward the renormalized field approach that helps us cast the linear $\mathcal{PT}$-symmetric dimer as a nonlinear dimer with imaginary Kerr nonlinearity, either in a self- or a cross-modulation scheme, and allowing us to realize an asymptotic behavior that depends just on the gain to coupling strength ratio of the device. In Section 3, we will discuss the linear $\mathcal{PT}$-symmetric dimer when both waveguides show the same effective self-focusing Kerr nonlinearity. We will show the stable nonlinear modes of the device, discuss its dynamics in terms of the passive Kerr two-waveguide coupler that allows for coherent and localized oscillations between the waveguide field modes and show that the inclusion of symmetric gain and loss breaks these dynamics, producing localization in the gain waveguide without showing an asymptotic behavior. In Section 4, we will extend the linear $\mathcal{PT}$-symmetric dimer to planar N-waveguide couplers using finite dimensional matrix representations of a complexified version of $SU(2)$. We have previously shown [48] that the $\mathcal{PT}$-symmetric dimer and its extensions to planar N-waveguide couplers possess an $SO(2,1)$ symmetry realized in a finite dimensional non-unitary irreducible representation. This representation is accomplished through complexification of $SU(2)$, $\{{\mathbb{J}}_x,{\mathbb{J}}_y,{\mathbb{J}}_z\}\to \{{\mathbb{J}}_x,i{\mathbb{J}}_y,i{\mathbb{J}}_z\}\equiv \{{\mathbb{K}}_z,{\mathbb{K}}_x,{\mathbb{K}}_y\}$, and it allows us to provide the dispersion relation and a closed form analytic propagator, which have the same regimes and dynamics found for the dimer. We will show that the renormalized field approach provides us with an asymptotic behavior that is independent of the initial field distribution and depends just on the waveguide number and the effective gain to coupling ratio. Then, in Section 5, we will introduce a modified version of the Ehrenfest theorem suitable for non-Hermitian Hamiltonians and show how it can help us define the dynamics of an $N^2$-dimensional generalized Stokes vector for the planar N-waveguide couplers discussed in Section 4. In Section 6, we will go back to the dimer, but consider the propagation of quantum fields. In the quantum regime, spontaneous generation and absorption of electromagnetic radiation should be considered when using media with linear gain or loss. We will show the solution for the quantum linear $\mathcal{PT}$-symmetric dimer and discuss the generation of light from a vacuum due to spontaneous processes in the absence of fields impinging at the device. Finally, we will produce a brief summary and discuss future avenues regarding non-Hermitian optical systems.

2. Linear $\mathcal{PT}$-Symmetric Dimer

A two-waveguide coupler can be described by the following differential equation system via mode coupling theory,

$$\begin{array}{c}\hfill -i{\partial }_z\left(\begin{array}{c}{E}_1(z)\\ E_2(z)\end{array}\right)=\left(\begin{array}{cc}{n}_1& g\\ g& n_2\end{array}\right)\left(\begin{array}{c}{E}_1(z)\\ E_2(z)\end{array}\right)\end{array}$$

(4)

where the complex numbers $E_j$ are the field amplitudes at each waveguide, the parameters $n_j$ and g, which in the most general case can be complex, are the effective refractive indices and waveguide coupling (the latter will be assumed real from now on, in that order) and the operator ${\partial }_z$ stands for the derivative with respect to the propagation distance z. Note that using field amplitudes of the form:

$$E_j(z)=e^{i{n}_{+}z}{\mathcal{E}}_j(z)$$

(5)

with $n_{\pm }=(n_1\pm n_2)/2$, reduces the system to an effective mode coupling matrix,

$$\begin{array}{c}\hfill -i{\partial }_z\left(\begin{array}{c}{\mathcal{E}}_1(z)\\ {\mathcal{E}}_2(z)\end{array}\right)=\left(\begin{array}{cc}{n}_{-}& g\\ g& -n_{-}\end{array}\right)\left(\begin{array}{c}{\mathcal{E}}_1(z)\\ {\mathcal{E}}_2(z)\end{array}\right)\end{array}$$

(6)

The mode coupling matrix has eigenvalues:

$$\begin{array}{c}\hfill {\epsilon }_{\pm }=\pm \sqrt{g^2+n_{-}^2}\end{array}$$

(7)

The eigenvalues are real for the case of the identical real part of the effective refractive indices, $\Re (n_{-})=0$, and the imaginary part less than the value of the effective coupling, $\Im (n_{-})<g$. The eigenvalues degenerate in the case $\Re (n_{-})=0$ and $\Im (n_{-})=g$. They are purely imaginary for $\Re (n_{-})=0$ and $\Im (n_{-})>g$ and, finally, complex elsewhere. This general non-Hermitian dimer has a rich structure that deserves further attention, but right now, we are interested in just the $\mathcal{PT}$-symmetric case.

In quantum mechanics, $𝒫$$𝒯$-symmetry refers to space-time reflection symmetry [6]. In discrete optical couplers, we can consider waveguide permutation and propagation inversion as equivalent to space and time reflection, respectively. Then, in order to recover the standard linear $\mathcal{P}\mathcal{T}$-symmetric dimer, we need to work with waveguides that have the same effective refractive indices, $\Re (n_1)=\Re (n_2)$, such that $n_{-}$ is purely imaginary, and we can write a differential set [48],

$$\begin{array}{c}\hfill -i{\partial }_{\zeta }\left(\begin{array}{c}{\mathcal{E}}_1(\zeta )\\ {\mathcal{E}}_2(\zeta )\end{array}\right)=\left(\begin{array}{cc}i\gamma & 1\\ 1& -i\gamma \end{array}\right)\left(\begin{array}{c}{\mathcal{E}}_1(\zeta )\\ {\mathcal{E}}_2(\zeta )\end{array}\right)\end{array}$$

(8)

where the permutation of the waveguides, $\pm i\gamma \to \mp i\gamma$, and propagation reversal, $\zeta \to -\zeta$, leaves the system invariant. Note that we have used the effective coupling parameter to scale the propagation distance, $\zeta =gz$, such that we deal with a single parameter given by the effective refractive index to coupling ratio, $\gamma =\Im (n_{-})/g$ with $\Re (n_{-})=0$. This mathematical model is equivalent to considering waveguides with effective pure linear loss and gain, $i\gamma$ and $-i\gamma$, in that order. In optics, the linear $\mathcal{PT}$-symmetric dimer has been experimentally demonstrated in passive lossy waveguides [49,50] (Figure 1a,b), as well as active, pumped waveguides [51] (Figure 1c) and pumped whispering-gallery mode microcavities [52,53] (Figure 1d) with linear gain. The experimental demonstration of $\mathcal{PT}$-symmetric devices is not limited to optical resonators; they have also been realized with operational amplifiers in electronics [54].

2.1. Quantum Mechanics, Linear Algebra Approach

Let us try to address propagation in the linear $\mathcal{PT}$-symmetric dimer. First, note that we can cast its coupled differential equation system, Equation (8), in a vector form,

$$\begin{array}{c}\hfill -i{\partial }_{\zeta }|\mathcal{E}(\zeta )\rangle =\mathbb{H}|\mathcal{E}(\zeta )\rangle \end{array}$$

(9)

similar to Schrödinger equation. Now, the coupling matrix takes the role of a Hamiltonian,

$$\begin{array}{ccc}\hfill \mathbb{H}& =& i\gamma {\sigma }_z+{\sigma }_x\hfill \end{array}$$

(10)

where the operators ${\sigma }_j$ with $j=x,y,z$ are Pauli matrices. We will keep the standard sign convention used in classical optics, although this Hamiltonian may confuse those working in quantum mechanics, as it implies that a parallel spin-field configuration is the most energetic configuration. In any given case, it is straightforward to find the propagator for this ζ-independent Schrödinger-like equation [48],

$$\begin{array}{}(11)& \hfill \mathbb{U}(\zeta )& =& e^{i\mathbb{H}\zeta }\hfill (12)& & =& \cos (\Omega \zeta )\mathrm{𝟙}+\frac{i}{\Omega }\sin (\Omega \zeta )\mathbb{H}\hfill \end{array}$$

such that the propagated fields through the device are given in terms of the initial field configuration,

$$\begin{array}{c}\hfill |\mathcal{E}(\zeta )\rangle =\mathbb{U}(\zeta )|\mathcal{E}(0)\rangle \end{array}$$

(13)

with the dispersion relation given by:

$$\begin{array}{c}\hfill \Omega =\sqrt{1-{\gamma }^2}\end{array}$$

(14)

which can be real for $\gamma <1$, zero for $\gamma =1$ or purely imaginary for $\gamma >1$. Note that the eigenvalues of the coupling matrix are given by $\pm \Omega$, and they become fully degenerate at $\gamma =1$. Note that this degeneracy is different from that in Hermitian systems; here, the actual dimension of the system becomes one at this point known in the theory of linear operators as the Kato exceptional point [55]. Figure 2 shows the behavior of the coupling matrix eigenvalues as a function of the gain to coupling ratio as they go from purely real, Figure 2a, degenerate to zero, Figure 2b, and become purely imaginary, Figure 2c.

It is straightforward to show that in the $\mathcal{PT}$-symmetric regime, where the eigenvalues of the mode coupling matrix are real, the device behaves as an asymmetric oscillator with amplification,

$$\begin{array}{ccc}\hfill \mathbb{U}(\zeta )& =& \cos (\Omega \zeta )\mathrm{𝟙}+\frac{i}{\Omega }\sin (\Omega \zeta )\mathbb{H},\gamma <1\hfill \end{array}$$

(15)

Once the $\mathcal{P}\mathcal{T}$-symmetry is broken, we have two distinct cases, the fully-degenerate one where both eigenvalues are zero and the device shows amplification ruled by a power law,

$$\begin{array}{ccc}\hfill \mathbb{U}(\zeta )& =& \mathrm{𝟙}+i\zeta \mathbb{H},\gamma =1\hfill \end{array}$$

(16)

and the case of purely imaginary eigenvalues, where the amplification is exponential,

$$\begin{array}{ccc}\hfill \mathbb{U}(\zeta )& =& \cosh (|\Omega |\zeta )\mathrm{𝟙}+\frac{i}{|\Omega |}\sinh (|\Omega |\zeta )\mathbb{H},\gamma >1\hfill \end{array}$$

(17)

Figure 3 shows the absolute field amplitude propagating through a $\mathcal{PT}$-symmetric dimer when light impinges just at the first waveguide in a device with parameters in the regime with real eigenvalues, Figure 3a, fully degenerate eigenvalues, Figure 3b, and imaginary eigenvalues, Figure 3c. Now, while we find this algebraic approach short, elegant and elucidating, it is not the only available method to infer the properties of the $\mathcal{PT}$-symmetric dimer.

2.2. Nonlocal Oscillator, Partial Differential Equation Approach

It is also possible to decouple the dimer differential equation set, Equation (8), into two second order differential equations that have the same form of the one-dimensional wave equation [48],

$$\begin{array}{c}\hfill {\partial }_{\zeta }^2{\mathcal{E}}_k(\zeta )+{\Omega }^2{\mathcal{E}}_k(\zeta )=0,k=1,2\end{array}$$

(18)

The boundary conditions are given by the initial fields impinging the device, ${\mathcal{E}}_k(0)$, and their first derivatives,

$$\begin{array}{ccc}\hfill {\partial }_{\zeta }{\mathcal{E}}_1(0)& =& -\gamma {\mathcal{E}}_1(0)+i{\mathcal{E}}_2(0)\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill {\partial }_{\zeta }{\mathcal{E}}_2(0)& =& \gamma {\mathcal{E}}_2(0)+i{\mathcal{E}}_1(0)\hfill \end{array}$$

(20)

obtained from the dimer mode coupling equation, Equation (8). Note that the latter takes us far from the standard one-dimensional wave equation, because it refers to nonlocal media involving linear loss and gain. Nevertheless, we can follow the standard partial differential equation approach.

In the $\mathcal{PT}$-symmetric case, the characteristic equation is positive, ${\Omega }^2>0$, and we can write the second order differential equation as that of a standard oscillator,

$$\begin{array}{c}\hfill {\partial }_{\zeta }^2{\mathcal{E}}_k(\zeta )+{\Omega }^2{\mathcal{E}}_k(\zeta )=0\end{array}$$

(21)

thus, we can use its well-known solution plus our particular boundary conditions to obtain the propagated fields,

$$\begin{array}{ccc}\hfill {\mathcal{E}}_1(\zeta )& =& {\mathcal{E}}_1(0)\cos (|\Omega |\zeta )-\frac{1}{|\Omega |}[\gamma {\mathcal{E}}_1(0)-i{\mathcal{E}}_2(0)]\sin (|\Omega |\zeta )\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill {\mathcal{E}}_2(\zeta )& =& {\mathcal{E}}_2(0)\cos (|\Omega |\zeta )+\frac{1}{|\Omega |}[\gamma {\mathcal{E}}_2(0)+i{\mathcal{E}}_1(0)]\sin (|\Omega |\zeta )\hfill \end{array}$$

(23)

These fields allow us to describe the dimer as an asymmetric periodic oscillator with amplification. It is not a harmonic oscillator due to the boundary condition on the first derivatives. In the fully-degenerate case, where the characteristic equation is equal to zero, ${\Omega }^2=0$,

$$\begin{array}{c}\hfill {\partial }_{\zeta }^2{\mathcal{E}}_k(\zeta )=0\end{array}$$

(24)

we can think of light propagating through the dimer as a free particle through some nonlocal media with linear unitary gain and loss, which yields amplification following a power law,

$$\begin{array}{ccc}\hfill {\mathcal{E}}_1\left(\zeta \right)& =& {\mathcal{E}}_1\left(0\right)-\left[{\mathcal{E}}_1\left(0\right)-i{\mathcal{E}}_2\left(0\right)\right]\zeta \hfill \end{array}$$

(25)

$$\begin{array}{ccc}\hfill {\mathcal{E}}_2\left(\zeta \right)& =& {\mathcal{E}}_2\left(0\right)+\left[{\mathcal{E}}_2\left(0\right)+i{\mathcal{E}}_1\left(0\right)\right]\zeta \hfill \end{array}$$

(26)

Finally, in the broken symmetry case, the characteristic equation is negative, ${\Omega }^2<0$, and we can write the second order differential equation as an inverted oscillator,

$$\begin{array}{c}\hfill {\partial }_{\zeta }^2{\mathcal{E}}_k\left(\zeta \right)-{|\Omega |}^2{\mathcal{E}}_k\left(\zeta \right)=0\end{array}$$

(27)

that provides us with a device that amplifies initial fields following an exponential law,

$$\begin{array}{ccc}\hfill {\mathcal{E}}_1\left(\zeta \right)& =& {\mathcal{E}}_1\left(0\right)\cosh \left(|\Omega |\zeta \right)-\frac{1}{|\Omega |}\left[\gamma {\mathcal{E}}_1\left(0\right)-i{\mathcal{E}}_2\left(0\right)\right]\sinh \left(|\Omega |\zeta \right)\hfill \end{array}$$

(28)

$$\begin{array}{ccc}\hfill {\mathcal{E}}_2\left(\zeta \right)& =& {\mathcal{E}}_2\left(0\right)\cosh \left(|\Omega |\zeta \right)+\frac{1}{|\Omega |}\left[\gamma {\mathcal{E}}_2\left(0\right)+i{\mathcal{E}}_1\left(0\right)\right]\sinh \left(|\Omega |\zeta \right)\hfill \end{array}$$

(29)

All of these solutions are just the explicit form of the propagated field found earlier by a purely algebraic approach, Equations (13)–(17), and tap into the well-known one-dimensional wave equation with the difference that an effective nonlocal active medium is provided by the first derivative boundary conditions.

2.3. Nonlinear Oscillator, Renormalized Fields Approach

So far, we have managed to provide an algebraic propagator and to describe the field behavior in the three possible regimes of the linear $\mathcal{PT}$-symmetric dimer. Now, we can bring forward a complementary view that can give us asymptotic information of the broken symmetry phases. Let us define instantaneous renormalized fields [48],

$$\begin{array}{c}\hfill {\tilde{\mathcal{E}}}_k\left(\zeta \right)=\frac{{\mathcal{E}}_k\left(\zeta \right)}{\sqrt{{|{\mathcal{E}}_1\left(\zeta \right)|}^2+{|{\mathcal{E}}_2\left(\zeta \right)|}^2}}\end{array}$$

(30)

such that the total renormalized field intensity at each propagation distance is always the unit, ${\sum }_{k=1}^2{|{\tilde{\mathcal{E}}}_k\left(\zeta \right)|}^2=1$. In this picture, it is easier to realize that light intensity through the $\mathcal{PT}$-symmetric device behaves like a non-harmonic oscillator, Figure 4a. Furthermore, this allows us to conduct asymptotic analysis in the broken symmetry phases. In the fully-degenerate phase, $\gamma =1$, such that $\Omega =0$, it is possible to calculate the asymptotic behavior of the fields’ intensity as the scaled propagation distance goes to infinity and to find out that the renormalized optical power is balanced in both waveguides, independent of the initial field distribution,

$$\begin{array}{c}\hfill \underset{\zeta \to \infty }{\lim }{|{\tilde{\mathcal{E}}}_k\left(\zeta \right)|}^2=\frac{1}{2},\gamma =1,\Omega =0\end{array}$$

(31)

as shown in Figure 4b. In the broken symmetry regime, where the eigenvalues are purely imaginary, the asymptotic intensity distribution depends on the effective gain to coupling ratio, and we can include the previous result,

$$\begin{array}{ccc}\hfill \underset{\zeta \to \infty }{\lim }{|{\tilde{\mathcal{E}}}_1\left(\zeta \right)|}^2& =& \frac{1}{2\gamma }{\left(\gamma +\sqrt{{\gamma }^2-1}\right)}^{-1}\hfill \end{array}$$

(32)

$$\begin{array}{ccc}\hfill \underset{\zeta \to \infty }{\lim }{|{\tilde{\mathcal{E}}}_2\left(\zeta \right)|}^2& =& \frac{1}{2\gamma }\left(\gamma +\sqrt{{\gamma }^2-1}\right),\gamma \ge 1\hfill \end{array}$$

(33)

In other words, for a device long enough, in the broken symmetry region, the input field distribution has no effect on the output field intensity distribution. The latter is governed only by the gain to coupling ratio of the device, γ (Figure 4c).

Furthermore, this approach also allows us to show that the linear $\mathcal{PT}$-symmetric dimmer is the irreducible form of two equivalent reducible nonlinear models, one with linear loss and gain, as well as an imaginary nonlinearity [48],

$$\begin{array}{ccc}\hfill -i{\partial }_{\zeta }{\tilde{\mathcal{E}}}_1(\zeta )& =& {\tilde{\mathcal{E}}}_2(\zeta )+2i\gamma \left[1-{|{\tilde{\mathcal{E}}}_1(\zeta )|}^2\right]{\tilde{\mathcal{E}}}_1(\zeta )\hfill \end{array}$$

(34)

$$\begin{array}{ccc}\hfill -i{\partial }_{\zeta }{\tilde{\mathcal{E}}}_2(z)& =& {\tilde{\mathcal{E}}}_1(\zeta )-2i\gamma \left[1-{|{\tilde{\mathcal{E}}}_2(\zeta )|}^2\right]{\tilde{\mathcal{E}}}_2(\zeta )\hfill \end{array}$$

(35)

that has a form reminiscent of van der Pol oscillators [56] in the sense that there is a change of sign dependent on the intensity for the self-modulation of the waveguide fields. Actual $\mathcal{PT}$-symmetric-coupled van der Pol oscillators have been implemented in the laboratory, as well as their optical realizations [57]. The other form corresponds to an imaginary cross-nonlinearity,

$$\begin{array}{ccc}\hfill -i{\partial }_{\zeta }{\tilde{\mathcal{E}}}_1(\zeta )& =& {\tilde{\mathcal{E}}}_2(\zeta )+2i\gamma {|{\tilde{\mathcal{E}}}_2(\zeta )|}^2{\tilde{\mathcal{E}}}_1(\zeta )\hfill \end{array}$$

(36)

$$\begin{array}{ccc}\hfill -i{\partial }_{\zeta }{\tilde{\mathcal{E}}}_2(\zeta )& =& {\tilde{\mathcal{E}}}_1(\zeta )-2i\gamma {|{\tilde{\mathcal{E}}}_1(\zeta )|}^2{\tilde{\mathcal{E}}}_2(\zeta )\hfill \end{array}$$

(37)

We can follow a standard approach from non-Hermitian Bose–Hubbard dimers [58,59,60,61] used to deal with optical nonlinear $\mathcal{PT}$-symmetric dimers [28], and introduce a Stokes vector, $\overrightarrow{\tilde{\mathcal{S}}}=({\tilde{\mathcal{S}}}_x,{\tilde{\mathcal{S}}}_y,{\tilde{\mathcal{S}}}_z)$, with components given by,

$$\begin{array}{c}\hfill {\tilde{\mathcal{S}}}_j(\zeta )=\langle \tilde{\mathcal{E}}(\zeta )|{\sigma }_j|\tilde{\mathcal{E}}(\zeta )\rangle ,j=x,y,z\end{array}$$

(38)

Again, the matrices ${\sigma }_j$ are Pauli matrices; the notation $|\tilde{\mathcal{E}}(\zeta )\rangle$ is a column vector containing the renormalized field amplitudes equivalent to that defined in Section 2 for the field amplitudes; and the new conjugate transpose notation $\langle \tilde{\mathcal{E}}(\zeta )|={(|{\tilde{\mathcal{E}}}^{*}(\zeta )\rangle )}^T$ is a row vector with the conjugate renormalized field amplitudes as components. In terms of the renormalized field amplitudes,

$$\begin{array}{ccc}\hfill {\tilde{\mathcal{S}}}_x(\zeta )& =& {\tilde{\mathcal{E}}}_1^{*}(\zeta ){\tilde{\mathcal{E}}}_2(\zeta )+{\tilde{\mathcal{E}}}_1(\zeta ){\tilde{\mathcal{E}}}_2^{*}(\zeta )\hfill \end{array}$$

(39)

$$\begin{array}{ccc}\hfill {\tilde{\mathcal{S}}}_y(\zeta )& =& -i[{\tilde{\mathcal{E}}}_1^{*}(\zeta ){\tilde{\mathcal{E}}}_2(\zeta )-{\tilde{\mathcal{E}}}_1(\zeta ){\tilde{\mathcal{E}}}_2^{*}(\zeta )]\hfill \end{array}$$

(40)

$$\begin{array}{ccc}\hfill {\tilde{\mathcal{S}}}_z(\zeta )& =& {|{\tilde{\mathcal{E}}}_1(\zeta )|}^2-{|{\tilde{\mathcal{E}}}_2(\zeta )|}^2\hfill \end{array}$$

(41)

In the case at hand, where the nonlinearity is reducible, the Stokes vector norm is the unit and a constant of motion,

$$\begin{array}{c}\hfill \tilde{\mathcal{S}}(\zeta )=\sqrt{{\tilde{\mathcal{S}}}_x^2(\zeta )+{\tilde{\mathcal{S}}}_y^2(\zeta )+{\tilde{\mathcal{S}}}_z^2(\zeta )}={|{\tilde{\mathcal{E}}}_1(\zeta )|}^2+{|{\tilde{\mathcal{E}}}_2(\zeta )|}^2=1\end{array}$$

(42)

The Stokes vector components verify the following set of coupled differential equations:

$$\begin{array}{ccc}\hfill {\partial }_{\zeta }\tilde{\mathcal{S}}(\zeta )& =& 0\hfill \end{array}$$

(43)

$$\begin{array}{ccc}\hfill {\partial }_{\zeta }{\tilde{\mathcal{S}}}_x(\zeta )& =& 2\gamma {\tilde{\mathcal{S}}}_x(\zeta ){\tilde{\mathcal{S}}}_z(\zeta )\hfill \end{array}$$

(44)

$$\begin{array}{ccc}\hfill {\partial }_{\zeta }{\tilde{\mathcal{S}}}_y(\zeta )& =& 2{\tilde{\mathcal{S}}}_z(\zeta )[1+\gamma {\tilde{\mathcal{S}}}_y(\zeta )]\hfill \end{array}$$

(45)

$$\begin{array}{ccc}\hfill {\partial }_{\zeta }{\tilde{\mathcal{S}}}_z(\zeta )& =& -2\{{\tilde{\mathcal{S}}}_y(\zeta )+\gamma [1-{\tilde{\mathcal{S}}}_z^2(\zeta )]\}\hfill \end{array}$$

(46)

Note that these equations have a set of stable points only outside the $\mathcal{PT}$-symmetric regime, $\gamma \ge 1$, given by ${\tilde{\mathcal{S}}}_s(\zeta )=(0,-{\gamma }^{-1},-{\gamma }^{-1}\sqrt{{\gamma }^2-1})$. It is also possible to describe the asymptotic behavior of the renormalized Stokes vector components,

$$\begin{array}{ccc}\hfill \underset{\zeta \to \infty }{\lim }{\tilde{\mathcal{S}}}_x(\zeta )& =& 0\hfill \end{array}$$

(47)

$$\begin{array}{ccc}\hfill \underset{\zeta \to \infty }{\lim }{\tilde{\mathcal{S}}}_y(\zeta )& =& -\frac{1}{\gamma }\hfill \end{array}$$

(48)

$$\begin{array}{ccc}\hfill \underset{\zeta \to \infty }{\lim }{\tilde{\mathcal{S}}}_z(\zeta )& =& -\frac{\sqrt{{\gamma }^2-1}}{\gamma },\gamma \ge 1\hfill \end{array}$$

(49)

The Stokes vector approach allows us to visualize the field propagation as a trajectory on a unit sphere. Any initial condition, in the fully-degenerate and broken symmetry regimes, will converge asymptotically to the same stable point on the sphere of the differential equation, ${\lim }_{\zeta \to \infty }\overrightarrow{\tilde{S}}(\zeta )={\tilde{\mathcal{S}}}_s(\zeta )=(0,-{\gamma }^{-1},-{\gamma }^{-1}\sqrt{{\gamma }^2-1})$, in a dimer described by the gain to coupling ratio γ. Figure 5 shows the Stokes vector propagation related to the examples given in Figure 4 and an additional initial condition set to show that the asymptotic behavior is independent of the initial conditions. This asymptotical behavior suggest the use of this device as an unidirectional variable amplitude coupler.

3. Nonlinear $\mathcal{PT}$-Symmetric Dimer

Now, we will discuss a variation of the $\mathcal{PT}$-symmetric dimer where the waveguides have an additional Kerr nonlinearity,

$$\begin{array}{ccc}\hfill -i{\partial }_{\zeta }{\mathcal{E}}_1(\zeta )& =& \left[i\gamma +\kappa {|{\mathcal{E}}_1(\zeta )|}^2\right]{\mathcal{E}}_1(\zeta )+{\mathcal{E}}_2(\zeta )\hfill \end{array}$$

(50)

$$\begin{array}{ccc}\hfill -i{\partial }_{\zeta }{\mathcal{E}}_2(\zeta )& =& \left[-i\gamma +\kappa {|{\mathcal{E}}_2(\zeta )|}^2\right]{\mathcal{E}}_2(\zeta )+{\mathcal{E}}_1(\zeta )\hfill \end{array}$$

(51)

For the sake of simplicity, we will consider the effective Kerr nonlinearity to coupling strength ratio, κ, to be equal in both waveguides. Stability in this system has been formally discussed in the case of $\kappa =1$, and it was found that the zero equilibrium state is neutrally stable in the $\mathcal{PT}$-symmetric regime, $\gamma <1$, and that the total light intensity at the waveguides is bounded from above by the initial intensity amplified by an exponential gain, ${|{\mathcal{E}}_1(\zeta )|}^2+{|{\mathcal{E}}_2(\zeta )|}^2\le ({|{\mathcal{E}}_1(0)|}^2+{|{\mathcal{E}}_2(0)|}^2)e^{2\gamma \zeta }$ [35]. In order to deal with the dynamics, it is possible to introduce a Stokes vector defined, now, in terms of the field amplitudes [28],

$$\begin{array}{ccc}\hfill {\mathcal{S}}_x(\zeta )& =& {\mathcal{E}}_1^{*}(\zeta ){\mathcal{E}}_2(\zeta )+{\mathcal{E}}_1(\zeta ){\mathcal{E}}_2^{*}(\zeta )\hfill \end{array}$$

(52)

$$\begin{array}{ccc}\hfill {\mathcal{S}}_y(\zeta )& =& -i\left[{\mathcal{E}}_1^{*}(\zeta ){\mathcal{E}}_2(\zeta )-{\mathcal{E}}_1(\zeta ){\mathcal{E}}_2^{*}(\zeta )\right]\hfill \end{array}$$

(53)

$$\begin{array}{ccc}\hfill {\mathcal{S}}_z(\zeta )& =& {|{\mathcal{E}}_1(\zeta )|}^2-{|{\mathcal{E}}_2(\zeta )|}^2\hfill \end{array}$$

(54)

Thus, the norm of the redefined Stokes vector, the total intensity at the waveguides, is no longer a constant of motion,

$$\begin{array}{c}\hfill \mathcal{S}(\zeta )=\sqrt{{\mathcal{S}}_x^2(\zeta )+{\mathcal{S}}_y^2(\zeta )+{\mathcal{S}}_z^2(\zeta )}={|{\mathcal{E}}_1(\zeta )|}^2+{|{\mathcal{E}}_2(\zeta )|}^2\end{array}$$

(55)

and its dynamics,

$$\begin{array}{ccc}\hfill {\partial }_{\zeta }\mathcal{S}\left(\zeta \right)& =& -2\gamma {\mathcal{S}}_z\left(\zeta \right)\hfill \end{array}$$

(56)

$$\begin{array}{ccc}\hfill {\partial }_{\zeta }{\mathcal{S}}_x\left(\zeta \right)& =& \kappa {\mathcal{S}}_y\left(\zeta \right){\mathcal{S}}_z\left(\zeta \right)\hfill \end{array}$$

(57)

$$\begin{array}{ccc}\hfill {\partial }_{\zeta }{\mathcal{S}}_y\left(\zeta \right)& =& {\mathcal{S}}_z\left(\zeta \right)\left[2-\kappa {\mathcal{S}}_x\left(\zeta \right)\right]\hfill \end{array}$$

(58)

$$\begin{array}{ccc}\hfill {\partial }_{\zeta }{\mathcal{S}}_z\left(\zeta \right)& =& -2\left[{\mathcal{S}}_y\left(\zeta \right)+\gamma \mathcal{S}\left(\zeta \right)\right]\hfill \end{array}$$

(59)

will not be restricted to the unit sphere. These dynamics have been shown to be an optical simulation of a relativistic massless particle of negative charge in a pseudo-electromagnetic field [28]. Note that there is a set of stable points, ${\partial }_{\zeta }{\overrightarrow{S}}_s(\zeta )=0$, for the effective gain to coupling ratio $\gamma <1$,

$$\begin{array}{c}\hfill {\overrightarrow{S}}_s(\zeta )=\left\{\left(\frac{2}{\kappa },\pm \frac{2\gamma }{\kappa \sqrt{1-{\gamma }^2}},0\right),\left(\pm \frac{1}{\gamma }\sqrt{(1-{\gamma }^2)}|S_y(\zeta )|,{\mathcal{S}}_y(\zeta ),0\right)\right\}\end{array}$$

(60)

All of these stable nonlinear modes are waveguide fields of the form ${\mathcal{E}}_j=\mathcal{A}{e}^{i{\varphi }_j}$ with $\mathcal{A}\ge 0$ and phase difference constricted by the relation,

$$\begin{array}{c}\hfill \tan \left({\varphi }_1-{\varphi }_2\right)=\pm \frac{\gamma }{\sqrt{1-{\gamma }^2}}\end{array}$$

(61)

There are no stable points outside the $\mathcal{PT}$-symmetric regime, $\gamma \ge 1$.

In order to create intuition, let us start with the passive self-focusing two-waveguide coupler, $\gamma =0$, that has two constants of motion in the form of the Stokes vector norm, $S(\zeta )=1$, such that ${\partial }_{\zeta }S(\zeta )=0$, and the Hamiltonian-like quantity,

$$\begin{array}{c}\hfill \mathcal{H}=\frac{1}{2}\kappa {\mathcal{S}}_z^2(\zeta )+2{\mathcal{S}}_x(\zeta )\end{array}$$

(62)

Note that the system is integrable, and the conservation of the Stokes vector norm allows for the parametrization ${\mathcal{S}}_x(\zeta )=\sqrt{1-{\mathcal{S}}_z^2}\cos \varphi (\zeta )$ and ${\mathcal{S}}_y(\zeta )=\sqrt{1-{\mathcal{S}}_z^2}\sin \varphi (\zeta )$, such that we can write a Hamilton–Jacobi model,

$$\begin{array}{c}\hfill \mathcal{H}=\frac{1}{2}\kappa {\mathcal{S}}_z^2(\zeta )+2\sqrt{1-{\mathcal{S}}_z^2(\zeta )}\cos \varphi (\zeta ),{\partial }_{\zeta }{\mathcal{S}}_z(\zeta )=\frac{\partial \mathcal{H}}{\partial \varphi },{\partial }_{\zeta }\varphi (\zeta )=-\frac{\partial \mathcal{H}}{\partial {\mathcal{S}}_z}\end{array}$$

(63)

equivalent to that of a nonrigid pendulum or a Bose–Josephson junction [62,63]. This particular configuration allows for Rabi oscillations (Figure 6a) below the critical effective Kerr nonlinearity to coupling ratio, $\kappa =2$ (Figure 6b); above this critical value, the system can show both Rabi and Josephson oscillations (Figure 6c). In other words, the initial field amplitudes either coherently oscillate between the waveguides or localize at the waveguide where they were originally prepared, depending on both the initial field distribution and the effective Kerr nonlinearity of the device.

If we now include the symmetric gain and loss, $\gamma \ne 0$, the system is still integrable, as we can write two constants of motion [28],

$$\begin{array}{ccc}\hfill {\mathcal{S}}_C\left(\zeta \right)& =& \sqrt{{\left[\kappa {\mathcal{S}}_x\left(\zeta \right)-2\right]}^2+{\left[\kappa {\mathcal{S}}_y\left(\zeta \right)\right]}^2}\hfill \end{array}$$

(64)

$$\begin{array}{ccc}\hfill {\mathcal{S}}_J\left(\zeta \right)& =& S\left(\zeta \right)+\frac{2\gamma }{\kappa }\arctan \left[\frac{\kappa {\mathcal{S}}_x\left(\zeta \right)-2}{\kappa {\mathcal{S}}_y\left(\zeta \right)}\right]\hfill \end{array}$$

(65)

In the literature, it has been found that the system is stable in the interval $0\le \gamma <1$ for effective nonlinearity $\kappa =1$ [35], and numerical arguments have been given in the most general case [28]. Figure 7 shows how the dynamics of the passive nonlinear dimer are affected by the addition of a small effective gain to coupling ratio to the system. Below the critical nonlinearity for the passive system (Figure 7a), we can still find the coherent oscillation behavior of the linear $\mathcal{PT}$-symmetric dimer, but the former stable point is no longer a fixed point of the system. As the nonlinearity increases, we can see how device parameters and initial conditions start having an effect on the dynamics (Figure 7b), until a point where it is possible to have unstable light localization at the waveguide where light originally impinged (Figure 7c). Here, the constant of motion plays an important role as an accuracy test for the process of numerically solving the coupled nonlinear system.

Now, for a given effective Kerr nonlinearity to coupling ratio, $\kappa =1$, we can approximate a critical effective gain to coupling ratio, ${\gamma }_c\simeq {\pi }^{-1}$, where a change of dynamics occur. First, the reciprocity condition regarding the exchange of waveguides for a given initial field distribution breaks as we get closer to the approximate critical gain to coupling ratio (note the shift to the left in Figure 8a,b), then, after we cross the critical value, the field intensity at the gain waveguide gets localized and experiences an exponential gain, while the field in the lossy waveguide diminishes independently of the initial field amplitude distribution (Figure 8c). Note, the dynamics above the critical gain to coupling ratio do not tend to the fixed point of the system as in the linear device; the fields do not seem to show a constant behavior in the asymptotic limit. This localization with amplification in the gain waveguide suggests the use of these devices as optical diodes [28].

4. Linear $\mathcal{PT}$-Symmetric Planar $\mathit{N}$-Waveguide Coupler

As we said before, the linear dimer is the device with the smallest dimension showing $\mathcal{PT}$-symmetry. Its algebraic structure, Equation (10), and the fact that Pauli matrices are the two by two scaled matrix representation of the $SU(2)$ group in dimension two, ${\sigma }_k=2J_k$, suggest that a larger class of N-waveguide devices can be constructed with higher dimensional matrix representations of $SU(2)$. Therefore, we can construct a mode coupling equation set describing a device of N waveguides following Equation (8),

$$\begin{array}{c}\hfill i{\partial }_z|\mathcal{E}(z)\rangle =2\left(i{n}_{-}{\mathbb{J}}_z+g{\mathbb{J}}_x\right)|\mathcal{E}(z)\rangle \end{array}$$

(66)

Now, we choose to introduce a different scaled propagation, $\xi =2\zeta =2gz$, to recover a single parameter Schrödinger-like equation with the form,

$$\begin{array}{c}\hfill i{\partial }_{\xi }|\mathcal{E}(\xi )\rangle =\mathbb{H}|\mathcal{E}(\xi )\rangle ,\mathbb{H}=i\gamma {\mathbb{J}}_z+{\mathbb{J}}_x\end{array}$$

(67)

where the effective refractive index to coupling ratio is the same as in the standard dimer, $\gamma =n_{-}/g$. This mode coupling matrix with $i\gamma \in \mathbb{R}$ describes planar N-waveguide couplers with identical real effective refractive indices and underlying $SU(2)$ symmetry that show harmonic oscillator behavior [64]. These devices have been used to produce perfect state transfer in both the classical [65,66,67,68] and quantum regimes [69]. In order to deal with a linear $\mathcal{PT}$-symmetric N-waveguide coupler, we must consider pure imaginary effective refractive indices, $\gamma \in \mathbb{R}$, and thus, the realization of higher finite-dimensional non-unitary representation of $SO(2,1)$, implemented as a complexified version of $SU(2)$, is involved. In matrix form, the group generators have the following elements,

$$\begin{array}{ccc}\hfill {\left[{\mathbb{J}}_x\right]}_{m,n}& =& \frac{1}{2}\left[{\delta }_{m-1,n}\sqrt{n\left(2j-n+1\right)}+{\delta }_{m+1,n}\sqrt{m\left(2j-m+1\right)}\right]\hfill \end{array}$$

(68)

$$\begin{array}{ccc}\hfill {\left[{\mathbb{J}}_y\right]}_{m,n}& =& \frac{i}{2}\left[{\delta }_{m-1,n}\sqrt{n\left(2j-n+1\right)}-{\delta }_{m+1,n}\sqrt{m\left(2j-m+1\right)}\right]\hfill \end{array}$$

(69)

$$\begin{array}{ccc}\hfill {\left[{\mathbb{J}}_z\right]}_{m,n}& =& {\delta }_{m,n}\left(j-m+1\right)m,n=1,2,\dots ,N\hfill \end{array}$$

(70)

with the Kronecker delta given by ${\delta }_{m,n}$ and the Bargmann parameter by $j=(N-1)/2$. These matrices fulfill the commutation relation $\left[{\mathbb{J}}_i,{\mathbb{J}}_j\right]=i{ϵ}_{ijk}{\mathbb{J}}_z$, where ${ϵ}_{ijk}$ is the Levi–Civita symbol, and commute with the Casimir operator ${\mathbb{J}}^2={\mathbb{J}}_x^2+{\mathbb{J}}_y^2+{\mathbb{J}}_z^2=j(j+1)\mathrm{𝟙}$, $\left[{\mathbb{J}}_j,{\mathbb{J}}^2\right]=0$. In the standard differential form, this is equivalent to the coupled mode set,

$$-i{\partial }_{\xi }{\mathcal{E}}_k(\xi )=\frac{i\gamma }{2}\left(N-2k+1\right){\mathcal{E}}_k(\xi )+\frac{1}{2}\sqrt{(k-1)(N-k-1)}{\mathcal{E}}_{k-1}(\xi )+\frac{1}{2}\sqrt{k(N-k)}{\mathcal{E}}_{k+1}(\xi )$$

(71)

Following the vector notation, we can construct a field vector as,

$$\begin{array}{c}\hfill |\mathcal{E}(\xi )\rangle ={\displaystyle \sum _{k=1}^N}{\mathcal{E}}_k(\xi )|j,j-k+1\rangle \end{array}$$

(72)

where we can define the n-th element of the standard basis as:

$$\begin{array}{c}\hfill {\left[|j,m\rangle \right]}_n={\delta }_{j-m+1,n},m=-j,-j+1,\dots ,j-1,j,n=1,2,\dots ,N\end{array}$$

(73)

such that we can define more helpful generators with their corresponding actions,

$$\begin{array}{ccc}\hfill J_z|j,m\rangle & =& m|j,m\rangle \hfill \end{array}$$

(74)

$$\begin{array}{ccc}\hfill J_{\pm }|j,m\rangle & =& \sqrt{(m+1)(2j+1\mp m)}|j,m\pm 1\rangle \hfill \end{array}$$

(75)

where we have defined the ladder operators ${\mathbb{J}}_{\pm }={\mathbb{J}}_x\pm i{\mathbb{J}}_y$ that fulfill $\left[{\mathbb{J}}_z,{\mathbb{J}}_{\pm }\right]=\pm {\mathbb{J}}_{\pm }$. Following the Gilmore–Perelomov approach for $SU(2)$ [48,64], we can find the n-th eigenvalue of the mode coupling matrix,

$$\begin{array}{c}\hfill {\Omega }_n=(j-n+1)\Omega ,n=1,\dots ,N\end{array}$$

(76)

and obtain the same structure found for the dimer. All of the eigenvalues will be real numbers for $\gamma <1$ (Figure 9a), completely degenerate and equal to zero ${\Omega }_m=0$ for $\gamma =1$ (Figure 9b) and imaginary for $\gamma >1$ (Figure 9c). This so-called collapse of the eigenvalues is a direct consequence of the underlying symmetry [70].

Furthermore, we can provide the propagation matrix elements,

$$\begin{array}{c}{\left[\mathbb{U}\left(\xi \right)\right]}_{m,n}={\Omega }^{-2j}\sqrt{\left(\begin{array}{c}2j\\ m-1\end{array}\right)\left(\begin{array}{c}2j\\ n-1\end{array}\right)}{\left(\Omega \cos \frac{\Omega }{2}\xi -\gamma \sin \frac{\Omega }{2}\xi \right)}^{2(j+1)-m-n}\hfill \\ \times {\left(i\sin \frac{\Omega }{2}\xi \right)}^{m+n-2}{}_2{F}_1\left(1-m,1-n,-2j,{\Omega }^2{\csc }^2\frac{\Omega }{2}\xi \right),\gamma <1\hfill \end{array}$$

(77)

$$\begin{array}{c}{\left[\mathbb{U}\left(\xi \right)\right]}_{m,n}={(-2)}^{-2j}\sqrt{\left(\begin{array}{c}2j\\ m-1\end{array}\right)\left(\begin{array}{c}2j\\ n-1\end{array}\right)}{\left(\xi -2\right)}^{2(j+1)-m-n}{\left(i\xi \right)}^{m+n-2}\hfill \\ \times {}_2{F}_1\left(1-m,1-n,-2j,\frac{4}{{\xi }^2}\right),\gamma =1\hfill \end{array}$$

(78)

$$\begin{array}{c}{\left[\mathbb{U}\left(\xi \right)\right]}_{m,n}={\left(|\Omega |\right)}^{-2j}\sqrt{\left(\begin{array}{c}2j\\ m-1\end{array}\right)\left(\begin{array}{c}2j\\ n-1\end{array}\right)}{\left(|\Omega |\cosh \frac{|\Omega |}{2}\xi -\gamma \sinh \frac{|\Omega |}{2}\xi \right)}^{2(j+1)-m-n}\hfill \\ \times {\left(i\sinh \frac{|\Omega |}{2}\xi \right)}^{m+n-2}{}_2{F}_1\left(1-m,1-n,-{2j,|\Omega |}^2{\text{csch}}^2\frac{|\Omega |}{2}\xi \right),\gamma >1\hfill \end{array}$$

(79)

where the notations $\left(\begin{array}{c}a\\ b\end{array}\right)$ and ${}_{2}F_1(a,b,c,z)$ stand for the binomial coefficient and Gauss hypergeometric function, in that order. Again, we will have three distinct propagation behaviors as demonstrated for the dimer. These behaviors are simpler to visualize if we define renormalized field amplitudes,

$$\begin{array}{c}\hfill {\tilde{\mathcal{E}}}_k\left(\xi \right)=\frac{{\mathcal{E}}_k\left(\xi \right)}{\sqrt{{\sum }_{p=1}^N{|{\mathcal{E}}_p\left(\xi \right)|}^2}}\end{array}$$

(80)

Now, we can see periodic amplified oscillations in the $\mathcal{PT}$-symmetric regime, $\gamma <1$ (Figure 10a), amplification following a power law in the fully-degenerate regime, $\gamma =1$ (Figure 10b), and asymmetric amplification following an exponential law in the broken symmetry regime, $\gamma >1$ (Figure 10c).

Furthermore, we can derive the set of reducible coupled nonlinear equations for the renormalized field amplitudes,

$$\begin{array}{l}-i\partial {\tilde{\mathcal{E}}}_k\left(\xi \right)=\frac{i\gamma }{2}\left[\left(N-2k+1\right)+\sum _{p=1}^N\left(N-2p+1\right){|{\tilde{\mathcal{E}}}_p\left(\xi \right)|}^2\right]{\tilde{\mathcal{E}}}_k\left(\xi \right)+\\ \phantom{-i\partial {\tilde{\mathcal{E}}}_k\left(\xi \right)=}+\frac{1}{2}\sqrt{(k-1)(N-k-1)}{\tilde{\mathcal{E}}}_{k-1}\left(\xi \right)+\frac{1}{2}\sqrt{k(N-k)}{\tilde{\mathcal{E}}}_{k+1}\left(\xi \right),k=1,2,\dots ,N\end{array}$$

(81)

It is cumbersome, but possible to show that the asymptotic response of an N-waveguide coupler in the fully-degenerate and broken symmetry regimes is independent of the input field amplitude distribution,

$$\begin{array}{ccc}\hfill \underset{\xi \to \infty }{\lim }{|{\tilde{\mathcal{E}}}_k\left(\xi \right)|}^2& =& \frac{1}{{\left(2\gamma \right)}^{2j}}\left(\begin{array}{c}2j\\ k-1\end{array}\right){\left(\gamma +\sqrt{{\gamma }^2-1}\right)}^{2\left(k-j-1\right)},\gamma \ge 1\hfill \end{array}$$

(82)

which can be seen in Figure 10b for $\gamma =1$ and Figure 10c for $\gamma >1$. As expected from the mathematical description, in the fully-degenerate regime (Figure 10b), the extremal waveguides, those with major effective losses and gain, will transmit fields with smaller amplitudes than those in the central waveguides, because the asymptotic intensity distribution follows the binomial coefficient. In the broken symmetry regime (Figure 10c), the field intensity correlates with the strength of the gain or loss; the most intense field will travel through the waveguide with the larger effective gain and the less intense through the one with the larger effective loss.

The equivalent three-dimensional formulation of the Stokes vector for renormalized fields,

$$\begin{array}{c}\hfill {\tilde{\mathcal{J}}}_k\left(\xi \right)=\langle \tilde{\mathcal{E}}\left(\xi \right)|{\mathbb{J}}_k|\tilde{\mathcal{E}}\left(\xi \right)\rangle ,k=x,y,z\end{array}$$

(83)

yields components of the following form,

$$\begin{array}{ccc}\hfill {\tilde{\mathcal{J}}}_x\left(\xi \right)& =& \frac{1}{2}{\displaystyle \sum _{k=1}^{2j+1}}\sqrt{k(2j-k+1)}[{\tilde{\mathcal{E}}}_{k+1}^{*}{\tilde{\mathcal{E}}}_k+{\tilde{\mathcal{E}}}_k^{*}{\tilde{\mathcal{E}}}_{k+1}]\hfill \end{array}$$

(84)

$$\begin{array}{ccc}\hfill {\tilde{\mathcal{J}}}_y\left(\xi \right)& =& \frac{i}{2}{\displaystyle \sum _{k=1}^{2j+1}}\sqrt{k(2j-k+1)}[{\tilde{\mathcal{E}}}_{k+1}^{*}{\tilde{\mathcal{E}}}_k-{\tilde{\mathcal{E}}}_k^{*}{\tilde{\mathcal{E}}}_{k+1}]\hfill \end{array}$$

(85)

$$\begin{array}{ccc}\hfill {\tilde{\mathcal{J}}}_z\left(\xi \right)& =& {\displaystyle \sum _{k=1}^{2j+1}}\left(j-k+1\right){|{\tilde{\mathcal{E}}}_k\left(\xi \right)|}^2\hfill \end{array}$$

(86)

Here, the conserved variable is the Casimir operator,

$$\begin{array}{ccc}\hfill \tilde{\mathcal{C}}\left(\xi \right)& =& \langle \tilde{\mathcal{E}}\left(\xi \right)|\left[{\mathbb{J}}_x^2+{\mathbb{J}}_y^2+{\mathbb{J}}_z^2\right]|\tilde{\mathcal{E}}\left(\xi \right)\rangle =j(j+1)\hfill \end{array}$$

(87)

and it is important to emphasize that the norm of this three-dimensional Stokes vector is no longer a constant of motion,

$$\begin{array}{c}\hfill \tilde{\mathcal{J}}\left(\xi \right)=\sqrt{{\tilde{\mathcal{J}}}_x^2\left(\xi \right)+{\tilde{\mathcal{J}}}_y^2\left(\xi \right)+{\tilde{\mathcal{J}}}_z^2\left(\xi \right)}\end{array}$$

(88)

The reason behind this is that a complex vector of dimension N with unit norm can be represented as a point on the surface of a ball of unit radius in dimension $N^2-1$. For example, the renormalized fields through a two-waveguide coupler, two-dimensional complex vector of unit norm, can be represented on the surface of a three-dimensional ball; in other words, a two-dimensional sphere, where the Stokes vector norm is a constant of motion,

$$\begin{array}{c}\hfill \tilde{\mathcal{J}}\left(\xi \right)=j\left[{|{\tilde{\mathcal{E}}}_1\left(\xi \right)|}^2+{|{\tilde{\mathcal{E}}}_2\left(\xi \right)|}^2\right]=\frac{1}{2},j=\frac{1}{2}\end{array}$$

(89)

In general, we should use the surface of a $(N^2-1)$-ball of unit radius in order to describe properly the field amplitudes propagating through an N-waveguide coupler. This does not make it simpler to visualize the dynamics, so we favor a projection from $(N^2-1)$-dimensional to three-dimensional space with the price of loosing the unit norm for all cases, but $N=2$, where we can write ${\tilde{\mathcal{J}}}_k(\xi /2)=2{\tilde{\mathcal{S}}}_k(\zeta )$ with $j=1/2$. Note that this projection also allows us to derive an asymptotic expression for the z-component of the Stokes vector via Equation (82), and heuristically propose the rest,

$$\begin{array}{ccc}\hfill \underset{\xi \to \infty }{\lim }{\tilde{\mathcal{J}}}_x\left(\xi \right)& =& 0\hfill \end{array}$$

(90)

$$\begin{array}{ccc}\hfill \underset{\xi \to \infty }{\lim }{\tilde{\mathcal{J}}}_y\left(\xi \right)& =& -\frac{1}{\gamma }j\hfill \end{array}$$

(91)

$$\begin{array}{ccc}\hfill \underset{\xi \to \infty }{\lim }{\tilde{\mathcal{J}}}_z\left(\xi \right)& =& -\frac{\sqrt{{\gamma }^2-1}}{\gamma }j,\hfill \end{array}$$

(92)

$$\begin{array}{ccc}\hfill \underset{\xi \to \infty }{\lim }\tilde{\mathcal{J}}\left(\xi \right)& =& j\hfill \end{array}$$

(93)

this was confirmed numerically over a random sample of initial states and gain to coupling ratios outside the $\mathcal{PT}$-symmetric regime, $j\in \left[1/2,10\right]$ and $\gamma \in \left[1,3\right]$. Figure 11 shows the propagation of the renormalized Stokes vector in a six-waveguide coupler, $j=5/2$, with parameters in the $\mathcal{PT}$-symmetric, $\gamma =0.5$ (Figure 11a), fully-degenerate, $\gamma =1$ (Figure 11b), and broken symmetry, $\gamma =1.5$ (Figure 11c), regimes for light impinging the first waveguide of the coupler in black. Furthermore, the renormalized Stokes vector propagation for an initial field amplitude distribution corresponding to the eigenstate of $J_x$ with eigenvalue $-j$ is plotted in red to show the asymptotic behavior outside the $\mathcal{PT}$-symmetric regime.

5. Non-Hermitian Ehrenfest Theorem and Generalized Stokes Vector

So far, we have studied the propagation of classical light through a class of $\mathcal{PT}$-symmetric devices with underlying $SO(2,1)$ symmetry that includes the linear and nonlinear $\mathcal{PT}$-symmetric dimer. We have seen that propagation on these linear and nonlinear devices is more involved than in the passive case, $\gamma =0$, but the linear algebra approach has helped us produce closed form propagators for the lineal case, and the renormalized fields approach and its Stokes vector representation have allowed us to find stationary states and to visualize the propagation dynamics. Here, we shall try to understand the underlying reasons for this more complex propagation behavior.

As we mentioned earlier, the propagation of light through tight-binding N-waveguide couplers can be modeled by mode coupling theory in a form similar to the Schrödinger equation, Equation (9). Thus, finding a propagator, Equation (13), provides us with the information of the complex field amplitudes as they propagate through each waveguide. In quantum mechanics, we can also ask about the propagation of the expectation value, which in the optical picture translates to the following,

$$\begin{array}{c}\hfill \mathcal{A}\left(z\right)=\langle \mathcal{E}\left(z\right)\left|\mathbb{A}\right(z\left)\right|\mathcal{E}\left(z\right)\rangle \end{array}$$

(94)

for the operator $\mathbb{A}(z)$, which can vary with propagation. This is exactly what has been done when studying the propagation of the Stokes vector defined as the expectation value of the $SU(2)$ generators, Equations (38) and (83), with the peculiarity that we used renormalized field amplitudes instead of just the field amplitudes in the linear cases. The Ehrenfest theorem relates the variation with the propagation of the mean value with the dynamics of the model,

$$\begin{array}{c}\hfill \frac{d}{dz}\mathcal{A}(z)=〈\frac{d}{dz}\mathbb{A}(z)〉\end{array}$$

(95)

In standard Hermitian quantum mechanics, the Heisenberg equation provides the equation of motion for the operator, but in our non-Hermitian model, we need to go beyond this. Any non-Hermitian operator $\mathbb{A}(z)$ can be decomposed,

$$\begin{array}{c}\hfill \mathbb{A}(z)={\mathbb{A}}_H(z)+{\mathbb{A}}_S(z)\end{array}$$

(96)

into Hermitian, ${\mathbb{A}}_H^{†}(z)={\mathbb{A}}_H(z)$, and skew-Hermitian, ${\mathbb{A}}_S^{†}(z)=-{\mathbb{A}}_S(z)$, parts. Thus, we can define a commutator between non-Hermitian operators,

$$\begin{array}{c}\hfill {\left[\mathbb{A}(z),\mathbb{B}(z)\right]}_{NH}=\mathbb{A}(z)\mathbb{B}(z)-{\mathbb{B}}^{†}(z){\mathbb{A}}^{†}(z)\end{array}$$

(97)

such that we can write a Heisenberg-like equation of motion ruled by a non-Hermitian Hamiltonian,

$$\begin{array}{c}\hfill \frac{d}{dz}\mathbb{A}(z)=i{\left[\mathbb{A}(z),\mathbb{H}(z)\right]}_{NH}+\frac{\partial }{\partial z}\mathbb{A}(z)\end{array}$$

(98)

and recover a non-Hermitian generalization of the Ehrenfest theorem for tight binding non-Hermitian waveguide couplers,

$$\begin{array}{c}\hfill \frac{d}{dz}\mathcal{A}(z)=i〈{\left[\mathbb{A}(z),\mathbb{H}(z)\right]}_{NH}〉+〈\frac{\partial }{\partial z}\mathbb{A}(z)〉\end{array}$$

(99)

Note that equivalent modified Heisenberg equations of motion have been derived following alternate methods in non-Hermitian quantum mechanics [71]. Let us consider, as an example, the N-waveguide coupler of the last section, Equation (67), where the effective mode-coupling matrix acting as the Hamiltonian has a Hermitian part, ${\mathbb{H}}_H={\mathbb{J}}_x$, and a skew-Hermitian part, ${\mathbb{H}}_S=i\gamma {\mathbb{J}}_z$. Here, the Ehrenfest theorem can be simplified to the expression,

$$\begin{array}{c}\hfill {\partial }_{\xi }\mathcal{A}(\xi )=i〈\left[\mathbb{A}(\xi ),{\mathbb{H}}_H\right]〉+i〈\left\{\mathbb{A}(\xi ),{\mathbb{H}}_S\right\}〉+〈{\partial }_{\xi }\mathbb{A}(\xi )〉\end{array}$$

(100)

where the standard commutator, $\left[\mathbb{A},\mathbb{B}\right]=\mathbb{A}\mathbb{B}-\mathbb{B}\mathbb{A}$, and anti-commutator, $\left\{\mathbb{A},\mathbb{B}\right\}=\mathbb{A}\mathbb{B}+\mathbb{B}\mathbb{A}$, have been used. We can see that the anti-commutator term in this expression will be proportional to the gain to coupling strength ratio, γ, and this is the culprit behind the more complex behavior of our general class of N-waveguide lattices with underlying complexified $SU(2)$ symmetry. Note that we can also use this result to derive the conserved quantities of the model by solving ${\partial }_{\xi }\mathcal{A}(\xi )=i〈\left[\mathbb{A}(\xi ),{\mathbb{H}}_H\right]〉+i〈\left\{\mathbb{A}(\xi ),{\mathbb{H}}_S\right\}〉+〈{\partial }_{\xi }\mathbb{A}(\xi )〉=0$.

As a practical example, let us derive the equations of motion for the Stokes vector for the linear $\mathcal{PT}$-symmetric dimer, ${\mathcal{S}}_k=\langle {\sigma }_k\rangle$ with $k=0,x,y,z$, where ${\sigma }_0$ is the identity matrix, given by Equations (52)–(55). Here, the Hermitian and skew-Hermitian parts of the mode coupling matrix are ${\mathbb{H}}_H={\sigma }_x$ and ${\mathbb{H}}_S=i\gamma {\sigma }_z$. Now, any linear Hermitian operator for this system can be written as the linear superposition of the matrices ${\sigma }_k$,

$$\begin{array}{c}\hfill \mathbb{A}(\zeta )=\sum _{k=0,x,y,z}{a}_k^{(\mathbb{A})}(\zeta ){\sigma }_{\mu }\end{array}$$

(101)

and the propagation of its expectation values, according to the non-Hermitian Ehrenfest theorem, are given by the following expression,

$${\partial }_{\zeta }\mathcal{A}(\zeta )=2\left[a_y(\zeta ){\mathcal{S}}_z(\zeta )-a_z(\zeta ){\mathcal{S}}_y(\zeta )\right]-2\gamma \left[a_0(\zeta ){\mathcal{S}}_z(\zeta )+a_z(\zeta ){\mathcal{S}}_0(\zeta )\right]+\sum _{k=0,x,y,z}{\mathcal{S}}_k(\zeta ){\partial }_{\zeta }{a}_k(\zeta )$$

(102)

Thus, noting that for the Stokes vectors, the coefficients are constant, $a_k^{({\sigma }_l)}(\zeta )={\delta }_{l,k}$, we can write the evolution for the components of the Stokes vector without field renormalization,

$$\begin{array}{ccc}\hfill {\partial }_{\zeta }{\mathcal{S}}_0(\zeta )& =& -2\gamma {\mathcal{S}}_z(\zeta )\hfill \end{array}$$

(103)

$$\begin{array}{ccc}\hfill {\partial }_{\zeta }{\mathcal{S}}_x(\zeta )& =& 0\hfill \end{array}$$

(104)

$$\begin{array}{ccc}\hfill {\partial }_{\zeta }{\mathcal{S}}_y(\zeta )& =& 2{\mathcal{S}}_z(\zeta )\hfill \end{array}$$

(105)

$$\begin{array}{ccc}\hfill {\partial }_{\zeta }{\mathcal{S}}_z(\zeta )& =& -2\left[{\mathcal{S}}_y(\zeta )+2\gamma {\mathcal{S}}_0(\zeta )\right]\hfill \end{array}$$

(106)

which are in complete agreement with what we obtain from the nonlinear $\mathcal{PT}$-symmetric dimer, Equations (56)–(59), if we kill the effective nonlinearity to coupling strength ratio, $\kappa =0$. Note that the total intensity, ${\mathcal{S}}_0(\zeta )$, is not conserved as expected from non-Hermitian dynamics. Note that, in this case, the total intensity, ${\mathcal{S}}_0(\zeta )$, coincides with the norm of the Stokes vector, $\mathcal{S}(\zeta )=\sqrt{{\mathcal{S}}_x^2(\zeta )+{\mathcal{S}}_y^2(\zeta )+{\mathcal{S}}_z^2(\zeta )}$, and we can recover the renormalized Stokes vector dynamics, Equations (43)–(46), if we define a renormalized Stokes vector, ${\tilde{\mathcal{S}}}_k={\mathcal{S}}_k/{\mathcal{S}}_0$, and use the equations of motion found here. If we were to find a constant of motion, ${\mathbb{S}}_c(\zeta )$, then its components should satisfy,

$$\begin{array}{ccc}\hfill {\partial }_{\zeta }{a}_0^{({\mathbb{S}}_c)}(\zeta )& =& 2\gamma a_z^{{\mathbb{S}}_c}(\zeta )\hfill \end{array}$$

(107)

$$\begin{array}{ccc}\hfill {\partial }_{\zeta }{a}_x^{({\mathbb{S}}_c)}(\zeta )& =& 0\hfill \end{array}$$

(108)

$$\begin{array}{ccc}\hfill {\partial }_{\zeta }{a}_y^{({\mathbb{S}}_c)}(\zeta )& =& 2a_z^{({\mathbb{S}}_c)}(\zeta )\hfill \end{array}$$

(109)

$$\begin{array}{ccc}\hfill {\partial }_{\zeta }{a}_z^{({\mathbb{S}}_c)}(\zeta )& =& 2\left[\gamma a_0^{({\mathbb{S}}_c)}(\zeta )-a_y^{({\mathbb{S}}_C)}(\zeta )\right]\hfill \end{array}$$

(110)

A particular solution to this set of equation is $a_k^{({\mathbb{S}}_c)}(\zeta )={\delta }_{k,x}$ in agreement with Equation (105)

For the general case of the planar N-waveguide coupler, the situation is far more complex as we are dealing with square matrices of dimension N. In order to construct any given Hermitian operator of this dimension, we need a basis with a total of $N^2$ matrices; these are provided by the standard unitary group of degree N, $SU(N)$, plus the identity. This way, we will work with a set of $N^2$ operators where the first four elements are the representation of $SU(2)$ in dimension N plus the unity, ${\mathbb{J}}_k$ with $k=0,x,y,z$, that form $U(N)$,

$$\begin{array}{c}\hfill \mathbb{A}(\xi )={\displaystyle \sum _{k=0}^{N^2-1}}{a}_k^{(\mathbb{A})}(\xi ){\mathbb{J}}_k\end{array}$$

(111)

where we have just implicitly made the change ${\mathbb{J}}_x={\mathbb{J}}_1$, ${\mathbb{J}}_y={\mathbb{J}}_2$, ${\mathbb{J}}_z={\mathbb{J}}_3$. Thus, if we define a generalized Stokes vector for the planar N-waveguide coupler, it will have dimension $N^2$, and the zeroth component will be the total intensity in the system,

$$\begin{array}{c}\hfill {\mathcal{J}}_0(\xi )={\displaystyle \sum _{k=1}^N}{|{\mathcal{E}}_k(\xi )|}^2\end{array}$$

(112)

but in this case, the zeroth component of the generalized Stokes vector, ${\mathcal{J}}_0(\zeta )$, is still the total intensity, but does not coincide with the norm of the generalized Stokes vector, $\mathcal{J}=\sqrt{{\sum }_k{|{\mathcal{J}}_k|}^2}$. Thus, a graphical representation on the sphere will just be a projection of the propagation dynamics occurring on a $(N^2-1)$-dimensional hypersphere as mentioned before. Note that for passive devices, $\gamma =0$, the propagation equations for the Stokes vectors do not involve any other functions, and the propagation dynamics is restricted to the subgroup $SU(2)$ of $U(N)$, recovering the results of [64].

6. Quantum $𝒫$$𝒯$-Symmetric Dimer

Let us turn our attention now to the propagation of nonclassical light. In the quantum regime, it is possible to describe two-waveguide couplers with the following effective Hamiltonian [72],

$$\begin{array}{c}\hfill \widehat{H}=n_1{\widehat{a}}_1^{†}{\widehat{a}}_1+n_2{\widehat{a}}_2^{†}{\widehat{a}}_2+g\left({\widehat{a}}_1^{†}{\widehat{a}}_2+{\widehat{a}}_1{\widehat{a}}_2^{†}\right)\end{array}$$

(113)

where we have kept the notation for the effective refractive indices and evanescent coupling strength, $n_j$ with $j=1,2$ and g, in that order. At this point, we can use the Schwinger two-boson representation of $SU(2)$ [73],

$$\begin{array}{c}\hfill J_x=\frac{1}{2}\left({\widehat{a}}_1^{†}{\widehat{a}}_2+{\widehat{a}}_1{\widehat{a}}_2^{†}\right),J_y=-\frac{i}{2}\left({\widehat{a}}_1^{†}{\widehat{a}}_2-{\widehat{a}}_1{\widehat{a}}_2^{†}\right),J_z=\frac{1}{2}\left({\widehat{a}}_1^{†}{\widehat{a}}_1-{\widehat{a}}_2^{†}{\widehat{a}}_2\right)\end{array}$$

(114)

to write an effective Hamiltonian with underlying $SU(2)$ symmetry,

$$\begin{array}{c}\hfill H=\omega J_z+J_x,\omega =\frac{1}{g}\left(n_1-n_2\right)\end{array}$$

(115)

that answers to the effective Schrödinger equation,

$$-i{\partial }_{\xi }|\mathcal{E}(\xi )\rangle =H|\mathcal{E}(\xi )\rangle ,\xi =2gz$$

(116)

In the single photon regime,

$$\begin{array}{c}\hfill |\mathcal{E}(\xi )\rangle ={\mathcal{E}}_1(\xi )|1,0\rangle +{\mathcal{E}}_2(\xi )|0,1\rangle \end{array}$$

(117)

we recover the differential equation set describing the standard two-waveguide coupler,

$$\begin{array}{c}\hfill -i{\partial }_{\xi }\left(\begin{array}{c}{\mathcal{E}}_1\left(\xi \right)\\ {\mathcal{E}}_2\left(\xi \right)\end{array}\right)=\left(\begin{array}{cc}\omega & 1\\ 1& -\omega \end{array}\right)\left(\begin{array}{c}{\mathcal{E}}_1\left(\xi \right)\\ {\mathcal{E}}_2\left(\xi \right)\end{array}\right)\end{array}$$

(118)

This approach suffices for the analysis of ideal dimers without gain or loses, $\omega \in \mathbb{R}$, where the total photon number,

$$\begin{array}{c}\hfill \widehat{n}={\widehat{a}}_1^{†}{\widehat{a}}_1+{\widehat{a}}_2^{†}{\widehat{a}}_2\end{array}$$

(119)

of the initial state, $\langle \widehat{n}\rangle =\langle \mathcal{E}(0)|\widehat{n}|\mathcal{E}(0)\rangle$, determines the dimension of the $SU(2)$ representation to be used, $N=\langle \widehat{n}\rangle +1=2j+1$, with the eigenbasis of $J_z$ given by the following,

$$\begin{array}{c}\hfill |j,m\rangle =|n_1,n_2\rangle ,j=\frac{n_1+n_2}{2},m=\frac{n_1-n_2}{2}\end{array}$$

(120)

such that we can use the results proposed for classical waveguides’ couplers with underlying $SU(2)$ symmetry [64] to calculate relevant quantities, like the mean photon number at each waveguide.

The inclusion of linear loses and gain is not a trivial matter, and it is simpler to discuss in the Heisenberg picture [74],

$$\begin{array}{c}\hfill \frac{d}{dz}\widehat{O}(z)=-i\left[\widehat{H},\widehat{O}(z)\right]+{\partial }_z\widehat{O}(z)\end{array}$$

(121)

where we have accounted for the change from time to distance propagation. Let us go straight to the $\mathcal{PT}$-symmetric dimer, with the identical real part of the refractive index and moving into a rotating frame,

$$\begin{array}{c}\hfill \left(\begin{array}{c}{\widehat{a}}_1\left(z\right)\\ {\widehat{a}}_2\left(z\right)\end{array}\right)=e^{i{n}_{+}z}\left(\begin{array}{c}{\widehat{o}}_1\left(z\right)\\ {\widehat{o}}_2\left(z\right)\end{array}\right)\end{array}$$

(122)

such that, again, we can define a scaled propagation, $\zeta =gz$, and include spontaneous processes arising from the quantum description of materials with linear loss or gain processes [74],

$$\begin{array}{c}\hfill \frac{d}{d\zeta }\left(\begin{array}{c}{\widehat{o}}_1\left(\zeta \right)\\ {\widehat{o}}_2\left(\zeta \right)\end{array}\right)=i\mathbb{H}\left(\begin{array}{c}{\widehat{o}}_1\left(\zeta \right)\\ {\widehat{o}}_2\left(\zeta \right)\end{array}\right)+\mathrm{𝟙}\left(\begin{array}{c}{\widehat{f}}_1\left(\zeta \right)\\ {\widehat{f}}_2\left(\zeta \right)\end{array}\right),\mathbb{H}=\left(\begin{array}{cc}i\gamma & 1\\ 1& -i\gamma \end{array}\right)\end{array}$$

(123)

where the first term in the right-hand side is related to propagation through the quantum two-waveguide coupler with linear loss and gain, $\gamma \in \mathbb{R}$, and the second term describes the Gaussian random processes of emission and absorption, a result arising from the linear materials in an equivalent treatment to that used in the quantum description of the laser [75],

$$\begin{array}{ccc}\hfill \langle {\widehat{f}}_1^{†}\left(\zeta \right){\widehat{f}}_1\left({\zeta }^{\prime }\right)\rangle & =& 0,\langle {\widehat{f}}_1\left(\zeta \right){\widehat{f}}_1^{†}\left({\zeta }^{\prime }\right)\rangle =2\gamma \delta (\zeta -{\zeta }^{\prime }),\hfill \end{array}$$

(124)

$$\begin{array}{ccc}\hfill \langle {\widehat{f}}_2^{†}\left(\zeta \right){\widehat{f}}_2\left({\zeta }^{\prime }\right)\rangle & =& 2\gamma \delta (\zeta -{\zeta }^{\prime }),\langle {\widehat{f}}_2\left(\zeta \right){\widehat{f}}_2^{†}\left({\zeta }^{\prime }\right)\rangle =0\hfill \end{array}$$

(125)

The formal solution for this differential equation yields the propagation of the annihilation operators,

$$\begin{array}{c}\hfill \left(\begin{array}{c}{\widehat{a}}_1\left(\zeta \right)\\ {\widehat{a}}_2\left(\zeta \right)\end{array}\right)=e^{i\mathbb{H}\zeta }\left(\begin{array}{c}{\widehat{a}}_1\left(0\right)\\ {\widehat{a}}_2\left(0\right)\end{array}\right)+{\int }_0^{\zeta }{e}^{i\mathbb{H}\left(\zeta -t\right)}\left(\begin{array}{c}{\widehat{f}}_1\left(t\right)\\ {\widehat{f}}_2\left(t\right)\end{array}\right)dt,\end{array}$$

(126)

where we have obviated the common phase factor $e^{i{n}_{+}\zeta /g}$ that does not play any important role. Note that we can use the propagator we already found for the classical dimer, $\mathbb{U}=e^{i\mathbb{H}\zeta }$ in Equation (11), for the first term in the right-hand side.

Now, in order to realize the effect of processes induced by the linear materials, let us focus on spontaneous generation in the absence of fields in both waveguides. In the classical case, there will be no light at all propagating through the waveguides, but in the quantum case, even with an initial vacuum state, we can calculate the spontaneous generation at each waveguide [74],

$$\begin{array}{ccc}\hfill S_1& =& 2\gamma {\int }_0^{\zeta }{\left|{\left[e^{i\mathbb{H}t}\right]}_{12}\right|}^{2}dt,\hfill \end{array}$$

(127)

$$\begin{array}{ccc}\hfill S_2& =& 2\gamma {\int }_0^{\zeta }{\left|{\left[e^{i\mathbb{H}t}\right]}_{22}\right|}^{2}dt.\hfill \end{array}$$

(128)

It shows in the symmetric regime, $\gamma <1$,

$$\begin{array}{ccc}\hfill S_1& =& \frac{1}{4{\Omega }^3}\left[2\Omega \zeta -\sin \left(2\Omega \zeta \right)\right],\hfill \end{array}$$

(129)

$$\begin{array}{ccc}\hfill S_2& =& \frac{1}{4{\Omega }^3}[2\Omega \zeta +(1-2{\gamma }^2)\sin \left(2\Omega \zeta \right)+4\gamma \Omega {\sin }^2\left(2\Omega \zeta \right)],\hfill \end{array}$$

(130)

a linear increase with a periodic modulation (Figure 12a). In the fully-degenerate case, $\gamma =1$, the spontaneous generation,

$$\begin{array}{ccc}\hfill S_1& =& \frac{1}{3}{\zeta }^3,\hfill \end{array}$$

(131)

$$\begin{array}{ccc}\hfill S_2& =& \frac{\zeta }{3}(3+3\zeta +{\zeta }^2),\hfill \end{array}$$

(132)

follows a cubic polynomial (Figure 12b), and in the broken symmetry regime, $\gamma >1$,

$$\begin{array}{ccc}\hfill S_1& =& \frac{1}{{4|\Omega |}^3}\left[-2|\Omega |\zeta +\sinh \left(2|\Omega |\zeta \right)\right],\hfill \end{array}$$

(133)

$$\begin{array}{ccc}\hfill S_2& =& \frac{1}{{4|\Omega |}^3}[-2|\Omega |\zeta -(1-2{\gamma }^2)\sinh \left(2|\Omega |\zeta \right)+4\gamma |\Omega |{\sinh }^2\left(2|\Omega |\zeta \right)],\hfill \end{array}$$

(134)

it shows exponential amplification (Figure 12c). Further discussion regarding the effect of spontaneous processes on the propagation of diverse nonclassical fields through a linear $\mathcal{PT}$-symmetric dimer can be found in Ref. [74].

7. Conclusions

We have presented a review of the $\mathcal{PT}$-symmetric dimer in its linear, nonlinear and quantum versions and show that it belongs to a symmetry class with underlying $SO(2,1)$ symmetry, realized as a complexification of the $SU(2)$ group, that allows the description of N-waveguide couplers. We have aimed to present a coherent narrative of the different approaches to the optical $\mathcal{PT}$-symmetric dimer and relate them to the underlying symmetry of the model. In doing this, we introduce the idea of using a non-Hermitian version of the Ehrenfest theorem to approach the propagation dynamics of waveguide couplers described by non-Hermitian mode coupling matrices.

The field is young, and there still exist fundamental open questions on the subject, such as the analytic determination of critical effective nonlinearity to coupling ratios for the Kerr nonlinear $\mathcal{PT}$-symmetric dimer; the need for a deeper understanding of the non-unitary finite dimensional representations of $SO(2,1)$, realized without resorting to the complexified $SU(2)$ representations; the generalization to propagation-dependent photonic systems together with the possible applications; just to mention a few that we hope to address in future work.