Jamming anomaly in PT-symmetric systems

We consider a spatially extended system consisting of two coupled elements, where the energy is gained in the first element and dissipated in the second one. To make stationary processes possible, the structure is designed so that the gain exactly balances the loss. The energy is assumed to be gained and dissipated locally, that is, there are no flows to or from infinity. Mathematically, this system is modelled by a parity-time (${ \mathcal P }{ \mathcal T }$) symmetric Schrödinger equation, linear or non-linear, with a two-centre localised complex potential and vanishing boundary conditions at infinity.

This bicentric structure is prototypical for a variety of physical settings, such as the electromagnetic field between the active and passive coupled parallel waveguides [1, 2] and microtoroid resonators [3], or the pumped and lossy atomic cells [4]. A similar gain–loss dimer is formed by the Bose–Einstein condensate trapped in two communicating wells, with one well leaking and the other one being loaded with atoms [5]. One more context is provided by the quasi-one-dimensional polaritonic waveguide. Here, the bicentric gain–loss structure can be modelled on the transistor switch realised experimentally in [6]. The localized injection of polaritons at one point and enhanced absorption at another point of the waveguide can be exploited as a means of the current control. A ${ \mathcal P }{ \mathcal T }$-symmetric system consisting of two coupled mechanical resonators, a damped nonlinear and a driven linear one, also deserves to be mentioned [7]. The system was proposed as a concept for an on-chip microscale phonon diode. (More ${ \mathcal P }{ \mathcal T }$-symmetric bicentric structures are discussed in a recent review [8] where one can find all the pertinent references.)

Stationary processes in the bicentric structure are sustained by the energy flux from the active to the passive element. We study the control of the flux—and therefore, control of the total gain and loss rate in the corresponding elements—by means of the variations of the gain–loss coefficient of the system. Intuitively, the increase of this gain and loss amplitude is expected to intensify the flux. The aim of the current paper is to show that the anomalous behaviour is also possible; namely, the flux may drop as a result of the gain–loss amplitude increase. For reasons explained below, we are referring to this phenomenon as the jamming anomaly.

The outline of the paper is as follows. After formulating the anomaly mathematically (section 2), we exemplify the general concept with two exactly solvable models. First, we demonstrate the anomalous behaviour in the linear Schrödinger equation with the ${ \mathcal P }{ \mathcal T }$-symmetric double-delta well potential (section 3) and in its non-linear counterpart (section 4). To show that the double well is not imperative for the anomaly to occur, we then consider a potential with a single well (the ${ \mathcal P }{ \mathcal T }$-symmetric Scarff potential, section 5).

Both potentials considered reveal jamming near their exceptional points or just before the corresponding eigenvalue immerses in the continuous spectrum. In section 6, we argue that any ${ \mathcal P }{ \mathcal T }$-symmetric potential with an exceptional point in its spectrum should exhibit the flux anomaly in the gain–loss parameter region adjacent to that point. To test whether the proximity to an exceptional point (or to the edge of the continuum) is not only sufficient but necessary for the anomaly to occur, we examine an exactly-solvable example of a potential with an entirely discrete spectrum free from exceptional points (section 7). The corresponding flux shows no anomalous behaviour.

The conclusions of this study are summarised in section 8.

2.1. Flux across the gain–loss interface

The system is described by the nonlinear Schrödinger equation,

$$\begin{eqnarray}&&{\rm{i}}{{\rm{\Psi }}}_{t}+{{\rm{\Psi }}}_{{xx}}-V(x){\rm{\Psi }}+g{\rm{\Psi }}| {\rm{\Psi }}{| }^{2}=0,\end{eqnarray} \tag{ 1 }$$

or its linear counterpart defined by setting g = 0 in (1). Here Ψ is either a dimensionless amplitude of the electric field or an order parameter—depending on whether one aims for the optical or boson-condensate interpretation of equation (1). In the context of atomic and polaritonic condensates, t denotes a properly scaled time. In the paraxial optics application, t measures the propagation distance while x is the coordinate transversal to the beam.

In equation (1), $V(x)$ is the ${ \mathcal P }{ \mathcal T }$-symmetric potential:

$$\begin{eqnarray*}&&V(x)=U(x)+{\rm{i}}{W}(x),\end{eqnarray*}$$

with $U(-x)=U(x)$ and $W(-x)=-W(x)$. A product of the nonhermitian quantum mechanics [9, 10], the ${ \mathcal P }{ \mathcal T }$-symmetric potentials are becoming increasingly relevant in optics and other applied disciplines. (See, e.g. the recent reviews [2, 8].)

For simplicity of analysis, we assume that x lies on the infinite line, $-\infty \lt x\lt \infty$. The right half of the x-line is characterised by loss and the left half by gain:

$$\begin{eqnarray*}&&W(x)\gt 0\ {\rm{if}}\ x\lt 0;\qquad W(x)\lt 0\ {\rm{if}}\ x\gt 0.\end{eqnarray*}$$

The infinite line can serve as an approximation for a ring-shaped structure (see figure 1), provided $U(x)$ and $W(x)$ are rapidly decaying functions with the decay ranges much shorter than the length of the ring. We consider localised solutions of equation (1)—that is, solutions with the asymptotic behaviour ${\rm{\Psi }}(x,t)\to 0$ as $x\to \pm \infty$.

Zoom In Zoom Out Reset image size

The integrals ${\int }_{-\infty }^{0}| {\rm{\Psi }}{| }^{2}{\rm{d}}{x}$ and ${\int }_{0}^{\infty }| {\rm{\Psi }}{| }^{2}{\rm{d}}{x}$ give the beam power (or the number of particles) captured in the regions with gain and loss, respectively. The power (the number) in the dissipative region varies according to

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}{\int }_{0}^{\infty }| {\rm{\Psi }}{| }^{2}{\rm{d}}{x}=J(t)+2{\int }_{0}^{\infty }| {\rm{\Psi }}{| }^{2}W(x){\rm{d}}{x},\end{eqnarray} \tag{ 2 }$$

where we have defined the flux across the gain–loss interface:

$$\begin{eqnarray}&&{\left.J={\rm{i}}({{\rm{\Psi }}}_{x}^{* }{\rm{\Psi }}-{{\rm{\Psi }}}^{* }{{\rm{\Psi }}}_{x})\right|}_{x=0}.\end{eqnarray} \tag{ 3 }$$

Equation (2) implies that two sources of the power variation in the dissipative domain are (a) the energy flux from the region with gain and (b) the losses suffered between x = 0 and $\infty$. Likewise, the variation of the power between x = 0 and $-\infty$, is given by the gain in that region less the flux through the origin.

In this paper, we focus on stationary regimes, ${\rm{\Psi }}(x,t)=\psi (x){{\rm{e}}}^{{\rm{i}}{\kappa }^{2}t}$. Here the real ${\kappa }^{2}$ represents the propagation constant in optics and chemical potential in the context of boson condensates. The spatial part of the separable solution Ψ satisfies

$$\begin{eqnarray}&&-{\psi }_{{xx}}+V(x)\psi -g\psi | \psi {| }^{2}=-{\kappa }^{2}\psi .\end{eqnarray} \tag{ 4 }$$

The boundary conditions

$$\begin{eqnarray}&&\psi (x)\to 0\quad {\rm{as}}\ x\to \pm \infty \end{eqnarray} \tag{ 5 }$$

require that ${\kappa }^{2}$ be taken positive.

In the stationary regime, the power density $| {\rm{\Psi }}{| }^{2}$ and the flux (3) are time-independent: $| {\rm{\Psi }}{| }^{2}=| \psi (x){| }^{2}$ and

$$\begin{eqnarray}&&{J={\rm{i}}({\psi }_{x}^{* }\psi -{\psi }^{* }{\psi }_{x})| }_{x=0}.\end{eqnarray} \tag{ 6 }$$

In this case equation (2) gives

$$\begin{eqnarray}&&J=-2{\int }_{0}^{\infty }| \psi {| }^{2}W(x)\;{\rm{d}}{x}\gt 0.\end{eqnarray} \tag{ 7 }$$

Assume that the rate of gain and loss is controlled by a single parameter $\gamma \gt 0;$ for instance, let

$$\begin{eqnarray*}&&W=\gamma { \mathcal W }(x),\end{eqnarray*}$$

where ${ \mathcal W }(x)$ is γ-independent. We will be referring to γ as the gain–loss amplitude.

In the linear case (g = 0), the constant $-{\kappa }^{2}$ arises as an eigenvalue of the Schrödinger operator (4). For the given γ, there is a sequence of eigenvalues $-{\kappa }_{n}^{2},n=0,1,2,\ldots$. When $\gamma =0$, the operator (4) is hermitian and all eigenvalues are real. The eigenvalues remain on the real line even when γ is made small nonzero; this is a fascinating consequence of the ${ \mathcal P }{ \mathcal T }$ symmetry [9, 10]. As γ is increased further, then, at some point ${\gamma }_{c}$ (commonly referred to as the exceptional point [11, 12]), one or more pairs of eigenvalues may merge and become complex. In the region above the exceptional point, the ${ \mathcal P }{ \mathcal T }$ symmetry is said to be spontaneously broken [3, 9, 10].

The present paper deals with the unbroken region, $\gamma \lt {\gamma }_{c}$. Here all eigenvalues $-{\kappa }_{n}^{2}$ are real, and each of the corresponding eigenfunctions ${\psi }_{n}(x)$ can be brought to the ${ \mathcal P }{ \mathcal T }$-symmetric form by a constant phase shift [10, 13]. That is to say, there is no loss of generality in assuming that

$$\begin{eqnarray}&&{\psi }_{n}(-x)={\psi }_{n}^{* }(x).\end{eqnarray} \tag{ 8 }$$

This will be our routine assumption in what follows.

The question that we concern ourselves with is how the interfacial flux J associated with the normalised eigenfunction ${\psi }_{n}(x)$, responds to the variation of γ.

A similar question can be posed for the nonlinear Schrödinger equation. Unlike a generic dissipative system, a typical ${ \mathcal P }{ \mathcal T }$-symmetric potential supports a one-parameter family of localised modes for every fixed γ in some interval [5, 8, 14–16]. The total power carried by the localised mode,

$$\begin{eqnarray*}&&P={\int }_{-\infty }^{\infty }| \psi {| }^{2}{\rm{d}}{x},\end{eqnarray*}$$

is a function of γ and κ, where κ is the parameter of the family. Setting the power to a certain fixed value, e.g. setting $P(\gamma ,\kappa )=1$, makes κ an (implicit) function of γ: $\kappa ={\kappa }_{P}(\gamma )$. Our aim is to find out how the flux J associated with the nonlinear mode with the parameter $\kappa ={\kappa }_{P}(\gamma )$ changes as γ is varied.

2.2. Toy model: the ${\boldsymbol{ \mathcal Q }}{\boldsymbol{ \mathcal T }}$-symmetric dimer

At first glance, the increase of the gain–loss amplitude γ should boost the flux: the greater the power growth rate in the active region and the faster the power loss in the dissipative domain, the more energy flows across the gain–loss interface.

To illustrate this intuitively appealing idea, assume that the active and lossy elements of our ${ \mathcal P }{ \mathcal T }$-symmetric doublet are point-like objects—that is, they do not have any internal structure. In this case the ${ \mathcal P }{ \mathcal T }$-symmetric system can be described by a two-component vector equation known as the nonlinear Schrödinger dimer [1, 2, 17]:

$$\begin{eqnarray*}{\rm{i}}{\dot{{\rm{\Psi }}}}_{1}+{{\rm{\Psi }}}_{2}+g| {{\rm{\Psi }}}_{1}{| }^{2}{{\rm{\Psi }}}_{1} & = & {\rm{i}}\gamma {{\rm{\Psi }}}_{1},\\ {\rm{i}}{\dot{{\rm{\Psi }}}}_{2}+{{\rm{\Psi }}}_{1}+g| {{\rm{\Psi }}}_{2}{| }^{2}{{\rm{\Psi }}}_{2} & = & -{\rm{i}}\gamma {{\rm{\Psi }}}_{2}.\end{eqnarray*}$$

Here ${{\rm{\Psi }}}_{1}$ and ${{\rm{\Psi }}}_{2}$ are complex amplitudes of the active and lossy modes, respectively. The powers carried by the two modes are $| {{\rm{\Psi }}}_{1}{| }^{2}$ and $| {{\rm{\Psi }}}_{2}{| }^{2}$, and the interfacial flux is given by

$$\begin{eqnarray*}&&J={\rm{i}}({{\rm{\Psi }}}_{1}{{\rm{\Psi }}}_{2}^{* }-{{\rm{\Psi }}}_{2}{{\rm{\Psi }}}_{1}^{* }).\end{eqnarray*}$$

Stationary solutions of the dimer have the form ${{\rm{\Psi }}}_{1}(t)={{\rm{e}}}^{i{\kappa }^{2}t}{\psi }_{1}$ and ${{\rm{\Psi }}}_{2}(t)={{\rm{e}}}^{i{\kappa }^{2}t}{\psi }_{2}$, where ${\psi }_{1}$ and ${\psi }_{2}$ are complex coefficients determined by initial conditions, and ${\kappa }^{2}$ is a real propagation constant. The stationary mode powers ${P}_{1}=| {\psi }_{1}{| }^{2}$ and ${P}_{2}=| {\psi }_{2}{| }^{2}$ are equal: ${P}_{\mathrm{1,2}}=P$. The corresponding flux is

$$\begin{eqnarray}&&J=2\gamma P.\end{eqnarray} \tag{ 9 }$$

The modes' common power is related to κ and γ by one of the following two equations:

$$\begin{eqnarray}&&P=\displaystyle \frac{1}{g}({\kappa }^{2}\pm \sqrt{1-{\gamma }^{2}}).\end{eqnarray} \tag{ 10 }$$

In order to determine the effect of the gain–loss amplitude variation on the flux (9), we consider the initial conditions ${\psi }_{\mathrm{1,2}}$ with the γ-independent mode powers: ${P}_{1}={P}_{2}$ ($=P$). In the linear case (g = 0) this simply means that the vector $({\psi }_{1},{\psi }_{2})$ is normalised; for example, we can let P = 1. In the nonlinear situation ($g\ne 0$) we can also choose P = 1; this choice selects one value of the propagation constant out of the family in (10): ${\kappa }^{2}=g\mp \sqrt{1-{\gamma }^{2}}$.

Differentiating (9) with respect to γ and keeping in mind that ${\rm{d}}{P}/{\rm{d}}\gamma =0$, gives

$$\begin{eqnarray*}&&\displaystyle \frac{{\rm{d}}{J}}{{\rm{d}}\gamma }=2P.\end{eqnarray*}$$

This is always positive—in agreement with our intuition.

2.3. Jamming anomaly in the ${\boldsymbol{ \mathcal Q }}{\boldsymbol{ \mathcal T }}$-symmetric Schrödinger equation

Returning to our spatially extended system, equation (7) yields

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{J}}{{\rm{d}}\gamma }=-2{\int }_{0}^{\infty }| \psi {| }^{2}\;{ \mathcal W }\;{\rm{d}}{x}-2\gamma {\int }_{0}^{\infty }\displaystyle \frac{\partial | \psi {| }^{2}}{\partial \gamma }{ \mathcal W }\;{\rm{d}}{x}.\end{eqnarray} \tag{ 11 }$$

The first term in the right-hand side of (11) is positive. The second term can be negative, but if γ is small, the sum shall nevertheless be positive, ${\rm{d}}{J}/{\rm{d}}\gamma \gt 0$, —as in the structureless example of section 2.2.

Note that the initial conditions giving rise to stationary solutions of (1) cannot be chosen arbitrarily—in particular, one is not free to choose $\partial | \psi {| }^{2}/\partial \gamma =0$. The dependence of $\psi (x)$ on γ is controlled by the nonlinear boundary-value problem (4). Therefore if γ is large enough, the second term in (11) may become dominant and so the anomalous regimes with ${\rm{d}}{J}/{\rm{d}}\gamma \lt 0$ cannot be ruled out a priori. Below, we produce several systems that do display this anomalous behaviour.

Before proceeding to the analysis of the specific examples, it is appropriate to mention that the flux from the active to the lossy region bears some analogy with the traffic flow. Let γ be the traffic density on a highway and J the traffic flow through it. When the traffic is building up, then, as long as J stays under the road capacity value, the greater the density, the greater the traffic flow. However once the capacity of the highway has been reached, any further increase of γ results in the drop of J. The traffic becomes congested.

With this analogy in mind, the regime where the inequality ${\rm{d}}{J}/{\rm{d}}\gamma \lt 0$ holds will be referred to as the jamming anomaly. (A more detailed discussion of the jamming metaphor appears in section 8.)

Our first choice is the ${ \mathcal P }{ \mathcal T }$-symmetric double-delta well potential

$$\begin{eqnarray}&&V(x)=-\left[\delta \left(x-\displaystyle \frac{L}{2}\right)+\delta \left(x+\displaystyle \frac{L}{2}\right)\right]+{\rm{i}}\gamma \left[\delta \left(x+\displaystyle \frac{L}{2}\right)-\delta \left(x-\displaystyle \frac{L}{2}\right)\right],\end{eqnarray} \tag{ 12 }$$

where $L\gt 0$ is the distance between the wells. The well on the left (the one at $x=-L/2$) gains and the one on the right ($x=L/2$) loses power.

The linear (g = 0) Schrödinger equation (4) with the potential (12) has one or two discrete eigenvalues. The corresponding eigenfunctions are given by the common expression [5, 18]:

$$\begin{eqnarray}\psi (x)=\left\{\begin{array}{ll}\tfrac{{{\rm{e}}}^{{\rm{i}}\phi }+{{\rm{e}}}^{\kappa L-{\rm{i}}\phi }}{2\sqrt{N}}{{\rm{e}}}^{\kappa x}, & x\leqslant -\tfrac{L}{2};\\ \tfrac{\mathrm{cosh}(\kappa x+{\rm{i}}\phi )}{\sqrt{N}}, & -\tfrac{L}{2}\leqslant x\leqslant \tfrac{L}{2};\\ \tfrac{{{\rm{e}}}^{-{\rm{i}}\phi }+{{\rm{e}}}^{\kappa L+i\phi }}{2\sqrt{N}}{{\rm{e}}}^{-\kappa x}, & x\geqslant \tfrac{L}{2}.\end{array}\right.\end{eqnarray} \tag{ 13 }$$

Here κ is a root of the transcendental equation

$$\begin{eqnarray}&&{{\rm{e}}}^{-2\kappa L}=\displaystyle \frac{{\gamma }^{2}+{(2\kappa -1)}^{2}}{{\gamma }^{2}+1},\end{eqnarray} \tag{ 14 }$$

and the constant angle is defined by

$$\begin{eqnarray*}&&{{\rm{e}}}^{2{\rm{i}}\phi }=\displaystyle \frac{2\kappa -1+{\rm{i}}\gamma }{1-{\rm{i}}\gamma }{{\rm{e}}}^{\kappa L}.\end{eqnarray*}$$

Letting

$$\begin{eqnarray*}&&N=\displaystyle \frac{\mathrm{cos}2\phi }{2}\left(L+\displaystyle \frac{1}{\kappa }\right)+\displaystyle \frac{{{\rm{e}}}^{\kappa L}}{2\kappa },\end{eqnarray*}$$

the eigenfunction (13) is normalised to unity:

$$\begin{eqnarray}&&{\int }_{-\infty }^{\infty }| \psi {| }^{2}\;{\rm{d}}{x}=1.\end{eqnarray} \tag{ 15 }$$

The analysis of the secular equation (14) can be carried out without any use of a computer [19]. The number of eigenvalues and the behaviour of the corresponding branches depend on the distance between the wells.

When $L\lt 1$, there is a single positive branch of eigenvalues $\kappa (\gamma )$ which decays, monotonically, as γ is increased from zero to ${\gamma }_{0}$ (figure 2, leftmost panel in the top row). As γ reaches ${\gamma }_{0}$,

$$\begin{eqnarray*}&&{\gamma }_{0}=\sqrt{2/L-1},\end{eqnarray*}$$

the exponent κ falls to zero and the eigenvalue $-{\kappa }^{2}$ immerses in the continuous spectrum.

Zoom In Zoom Out Reset image size

When L is taken between 1 and 2, then, depending on the interval of γ values, the function $\kappa (\gamma )$ has one or two branches (figure 2, second panel in the top row). As γ is raised from 0 to ${\gamma }_{c}$, one branch descends, monotonically, from ${\kappa }^{(b)}$ to ${\kappa }_{c}$. In the interval ${\gamma }_{0}\lt \gamma \lt {\gamma }_{c}$, there also is a monotonically growing branch. On this branch, κ increases from 0 to ${\kappa }_{c}$ as γ changes from ${\gamma }_{0}$ to ${\gamma }_{c}$. As γ reaches ${\gamma }_{c}$, the two eigenvalues merge and become complex.

Finally, when $L\geqslant 2$, the monotonically increasing and decreasing branches of $\kappa (\gamma )$ exist over the same interval $0\leqslant \gamma \leqslant {\gamma }_{c}$ (figure 2, the last two panels in the top row). As γ grows from 0 to ${\gamma }_{c}$, one branch of κ grows from ${\kappa }^{(a)}$ to ${\kappa }_{c}$ whereas the other one descends from ${\kappa }^{(b)}$ to ${\kappa }_{c}$.

The interfacial flux (3) associated with the eigenfunction (13), is

$$\begin{eqnarray}&&J=\displaystyle \frac{\kappa }{N}\mathrm{sin}2\phi .\end{eqnarray} \tag{ 16 }$$

The gain–loss coefficient can be expressed as an explicit function of the eigenvalue:

$$\begin{eqnarray*}&&\gamma =\sqrt{\displaystyle \frac{4\kappa (1-\kappa )}{1-{{\rm{e}}}^{-2\kappa L}}-1}.\end{eqnarray*}$$

Therefore the N and the $\mathrm{sin}2\phi$ factors in (16) can also be explicitly parametrised by κ. As a result, the dependence $J(\gamma )$ can be represented by parametric equations $\gamma =\gamma (\kappa ),J=J(\kappa )$. This dependence is shown in the panels making up the bottom row in figure 2. Highlighted are the intervals of γ where ${\rm{d}}{J}/{\rm{d}}\gamma \lt 0$.

The variation of J with γ is governed by

$$\begin{eqnarray*}&&\displaystyle \frac{{\rm{d}}{J}}{{\rm{d}}\gamma }=\displaystyle \frac{{\rm{d}}{J}/{\rm{d}}\kappa }{{\rm{d}}\gamma /{\rm{d}}\kappa }.\end{eqnarray*}$$

It is instructive to consider the situation where the eigenvalue $-{\kappa }^{2}$ approaches the continuous spectrum: $\kappa \to 0$. In this limit, we have

$$\begin{eqnarray*}&&\displaystyle \frac{{\rm{d}}{J}}{{\rm{d}}\kappa }=2{\gamma }_{0}+O(\kappa ),\quad \displaystyle \frac{{\rm{d}}\gamma }{{\rm{d}}\kappa }=\displaystyle \frac{L-1}{{\gamma }_{0}L}+O(\kappa ).\end{eqnarray*}$$

In case $L\lt 1$, the behaviour of the flux is anomalous: ${\rm{d}}{J}/{\rm{d}}\gamma \lt 0$. (See figure 2, leftmost panel in the bottom row.)

To explore jamming in the presence of nonlinearity, we let $g\ne 0$ in the stationary equation (4) with the same two-delta well potential (12). In what follows, the inter-well distance L and the nonlinearity strength g are assumed to be fixed.

Fixing a particular value of γ gives rise to a branch of solutions with a continuously varied κ, rather than a countable set of ${\kappa }_{n}$ as in the linear case. Since we are interested in the variation of the interfacial flux J due to factors other than the increase of the total power $P=\int | \psi {| }^{2}{\rm{d}}{x}$, we impose the condition

$$\begin{eqnarray}&&\int | \psi {| }^{2}{\rm{d}}{x}=1.\end{eqnarray} \tag{ 17 }$$

The condition (17) was imposed in the linear case as well—recall equation (15). There, it did not affect the eigenvalues κ and was only necessary for the gauging of the flux J. In contrast, when $g\ne 0$, condition (17) selects a particular κ from a continuous family, and therefore represents a nontrivial constraint. The nonlinear Schrödinger equation (4) and this constraint form a system of two simultaneous equations. We note, in passing, that condition (17) has a clear physical meaning in the context of the two-trap boson-gas interpretation of equation (4) [5].

The system (4) + (17) can be solved exactly—in the sense that its solution can be reduced to a system of two transcendental equations [19]. Instead of finding the 'eigenvalue' κ and 'eigenfunction' ψ for the given γ, it is convenient to assume that the parameter κ is given—and solve for ψ and γ. Having determined $J(\kappa )$ and $\gamma (\kappa )$, one can readily plot $J(\gamma )$ as a parametric curve.

The nonlinear mode is found [19] to be

$$\begin{eqnarray}&&\psi (x)=\sqrt{\displaystyle \frac{2}{g}}\kappa \;\varphi (x),\end{eqnarray} \tag{ 18a }$$

where

$$\begin{eqnarray}\varphi =\left\{\begin{array}{lr}\mathrm{sech}(\kappa x+\mu ){{\rm{e}}}^{-{\rm{i}}\chi }, & x\leqslant -L/2\\ \sqrt{\rho (x)}{{\rm{e}}}^{{\rm{i}}\theta (x)}, & -L/2\leqslant x\leqslant L/2\\ \mathrm{sech}(\kappa x-\mu ){{\rm{e}}}^{{\rm{i}}\chi }, & x\geqslant L/2\end{array}\right..\end{eqnarray} \tag{ 18b }$$

Here χ is a constant phase defined by $\chi =\theta (L/2)$.

In the middle range, $-L/2\leqslant x\leqslant L/2$, the modulus-squared and phase of are given by

$$\begin{eqnarray}&&\rho (x)=(\alpha -\beta ){{\rm{cn}}}^{2}(K-\sqrt{2\alpha +\beta -1}\kappa x,k)+\beta \end{eqnarray} \tag{ 19 }$$

and

$$\begin{eqnarray*}&&\theta (x)=\kappa \sqrt{\alpha \beta (\alpha +\beta -1)}{\int }_{0}^{x}\displaystyle \frac{{\rm{d}}{s}}{\rho (s)},\end{eqnarray*}$$

respectively. Here α and β are two constants satisfying $\alpha \geqslant \beta \geqslant 0$ and $\alpha +\beta \gt 1$.

In equation (19), ${\rm{cn}}$ is the Jacobi elliptic cosine and K is the complete elliptic integral of the first kind. The elliptic modulus of the Jacobi function ${\rm{cn}}(z,k)$ and the argument of $K(k)$, is expressible in terms of α and β:

$$\begin{eqnarray*}&&k=\sqrt{\displaystyle \frac{\alpha -\beta }{2\alpha +\beta -1}}.\end{eqnarray*}$$

The parameter μ in (18b) is also expressible via these constants:

$$\begin{eqnarray}&&{\mathrm{sech}}^{2}\left(\mu -\displaystyle \frac{\kappa L}{2}\right)=\beta +(\alpha -\beta ){\mathrm{cn}}^{2}(K-y,k).\end{eqnarray} \tag{ 20 }$$

Here we have denoted

$$\begin{eqnarray*}&&y=\displaystyle \frac{L\kappa }{2}\sqrt{2\alpha +\beta -1}.\end{eqnarray*}$$

The constants α and β are determined as two components of a vector root of the following system of two transcendental equations:

$$\begin{eqnarray}&&{\zeta }^{2}+\beta +(\alpha -\beta ){\mathrm{cn}}^{2}(K-y,k)-1=0,\end{eqnarray} \tag{ 21a }$$

$$\begin{eqnarray}&&{(1+\kappa \zeta )}^{2}-{\kappa }^{2}{S}^{2}=0.\end{eqnarray} \tag{ 21b }$$

In (21a)–(21b), S² and ζ are functions of α and β:

$$\begin{eqnarray*}{S}^{2} & = & \displaystyle \frac{(\alpha +\beta )(\alpha +\beta -1)-\alpha \beta }{\beta +(\alpha -\beta ){\mathrm{cn}}^{2}(K-y,k)}-\displaystyle \frac{\alpha \beta (\alpha +\beta -1)}{{[\beta +(\alpha -\beta ){\mathrm{cn}}^{2}(K-y,k)]}^{2}}\\ & & +1-\beta -(\alpha -\beta ){\mathrm{cn}}^{2}(K-y,k)\end{eqnarray*}$$

and

$$\begin{eqnarray*}&&\zeta =\displaystyle \frac{\alpha +\beta -1}{\sqrt{2\alpha +\beta -1}}y-1+\displaystyle \frac{g}{4\kappa }-\sqrt{2\alpha +\beta -1}\left\{E\left(\displaystyle \frac{\pi }{2},k\right)-E[\mathrm{am}(K-y),k]\right\},\end{eqnarray*}$$

where $E[\mathrm{am}(z),k]$ is the incomplete elliptic integral of the second kind,

$$\begin{eqnarray*}&&E[\mathrm{am}(z),k]={\int }_{0}^{\mathrm{am}(z)}\sqrt{1-{k}^{2}{\mathrm{sin}}^{2}\theta }\;{\rm{d}}\theta ={\int }_{0}^{z}{\mathrm{dn}}^{2}(w,k){\rm{d}}{w},\end{eqnarray*}$$

and $E\left(\tfrac{\pi }{2},k\right)$ is the complete elliptic integral:

$$\begin{eqnarray*}&&E\left(\displaystyle \frac{\pi }{2},k\right)={\int }_{0}^{\pi /2}\sqrt{1-{k}^{2}{\mathrm{sin}}^{2}\theta }\;{\rm{d}}\theta .\end{eqnarray*}$$

Note that the transcendental equations (21a)–(21b) involve an explicit dependence on κ but not on γ. Hence $\alpha =\alpha (\kappa ),\beta =\beta (\kappa )$. Once a root of (21a)–(21b) has been determined, the gain–loss coefficient is expressible as

$$\begin{eqnarray}&&\gamma (\kappa )=\displaystyle \frac{\sqrt{\alpha \beta (\alpha +\beta -1)}\kappa }{\beta +(\alpha -\beta ){\mathrm{cn}}^{2}(K-y,k)}.\end{eqnarray} \tag{ 22 }$$

Substituting the function (18a)–(18b) in (6) gives the interfacial flux J, also as a function of κ:

$$\begin{eqnarray}&&J(\kappa )=2\kappa \sqrt{\alpha \beta (\alpha +\beta -1)}.\end{eqnarray} \tag{ 23 }$$

Solving the transcendental system (21a)–(21b) for a range of κ-values, we obtain the $J(\gamma )$ dependence in the parametric form: $\gamma =\gamma (\kappa )$, $J=J(\kappa )$. This is plotted in figure 3, for a representative pair of L and g. The dependence of the flux on the gain–loss coefficient features a segment with a negative slope indicating the jamming anomaly.

Zoom In Zoom Out Reset image size

The vital property of a nonlinear mode is its stability. To classify the stability of the stationary solution (18a)–(18b), we consider a solution to equation (1) of the form

$$\begin{eqnarray*}&&{\rm{\Psi }}(x,t)={{\rm{e}}}^{{\rm{i}}{\kappa }^{2}t}\{\psi (x)+\epsilon \quad {{\rm{e}}}^{\lambda t}[p(x)+{\rm{i}}{q}(x)]\},\end{eqnarray*}$$

with p and q real. Linearising in small gives an eigenvalue problem

$$\begin{eqnarray}&&{ \mathcal H }\vec{p}=\lambda { \mathcal J }\quad \vec{p},\end{eqnarray} \tag{ 24 }$$

where ${ \mathcal H }$ and ${ \mathcal J }$ are 2 × 2 matrices, and $\vec{p}$ is a two-component vector-column:

$$\begin{eqnarray}{ \mathcal H } & = & -I\displaystyle \frac{{{\rm{d}}}^{2}}{{{\rm{d}}{x}}^{2}}+\left(\begin{array}{lr}{\kappa }^{2}+U & -W\\ W & {\kappa }^{2}+U\end{array}\right)-g\left(\begin{array}{lr}3{{ \mathcal R }}^{2}+{{ \mathcal I }}^{2} & -2{ \mathcal R }{ \mathcal I }\\ 2{ \mathcal R }{ \mathcal I } & 3{{ \mathcal I }}^{2}+{{ \mathcal R }}^{2}\end{array}\right),\\ { \mathcal J } & = & \left(\begin{array}{lr}0 & -1\\ 1 & 0\end{array}\right),\quad I=\left(\begin{array}{lr}1 & 0\\ 0 & 1\end{array}\right),\quad \vec{p}=\left(\begin{array}{c}p\\ q\end{array}\right).\end{eqnarray} \tag{ 25 }$$

In (25), ${ \mathcal R }$ and ${ \mathcal I }$ are the real and imaginary part of $\psi (x)$:

$$\begin{eqnarray*}&&\psi ={ \mathcal R }(x)+{\rm{i}}{ \mathcal I }(x).\end{eqnarray*}$$

The stationary solution ${\rm{\Psi }}={{\rm{e}}}^{{\rm{i}}{\kappa }^{2}t}\psi (x)$ is classified as unstable if the operator ${{ \mathcal J }}^{-1}{ \mathcal H }$ has (at least one) eigenvalue λ with a positive real part.

The eigenvalue problem (24) was solved numerically. The resulting stability and instability domains are demarcated in figure 3. The central conclusion of this computer analysis is that the stationary nonlinear mode (18a)–(18b) is stable in a sizeable interval of the gain–loss amplitudes.

Does one need to have two potential wells in order to observe the jamming effect? The aim of this section is to demonstrate the same phenomenon in a single-well potential—yet with the bicentric distribution of gain and loss.

This complex profile is known as the ${ \mathcal P }{ \mathcal T }$-symmetric Scarff II potential [20–22]:

$$\begin{eqnarray}&&V(x)=-q\;{{\rm{sech}}}^{2}\;x-{\rm{i}}\gamma \;{\rm{sech}}\;x\mathrm{tanh}x.\end{eqnarray} \tag{ 26 }$$

This time the potential well is centred at the origin, while the gain and loss are continuously distributed over the entire left (gain) and right (loss) semiaxis.

A sufficiently deep potential well (26) can support an arbitrarily large number of bound states. For simplicity, we focus on the two lowest eigenvalues pertaining to the single-hump eigenfunctions. By the direct substitution one can verify that the eigenvalues are given by

$$\begin{eqnarray*}{\kappa }_{1} & = & \displaystyle \frac{\sqrt{{\gamma }_{c}+\gamma }+\sqrt{{\gamma }_{c}-\gamma }}{2}-\displaystyle \frac{1}{2},\\ {\kappa }_{2} & = & \displaystyle \frac{\sqrt{{\gamma }_{c}+\gamma }-\sqrt{{\gamma }_{c}-\gamma }}{2}-\displaystyle \frac{1}{2},\end{eqnarray*}$$

where

$$\begin{eqnarray*}&&{\gamma }_{c}=q+\displaystyle \frac{1}{4},\end{eqnarray*}$$

and we have assumed that $\gamma \leqslant {\gamma }_{c}$. The corresponding eigenfunctions are

$$\begin{eqnarray}{\psi }_{1} & = & {\pi }^{-1/4}\sqrt{\displaystyle \frac{{\rm{\Gamma }}({\kappa }_{1}+1/2)}{{\rm{\Gamma }}({\kappa }_{1})}}{({\rm{sech}}x)}^{{\kappa }_{1}}\mathrm{exp}\left\{{\rm{i}}\displaystyle \frac{2{\kappa }_{2}+1}{2}\mathrm{arctan}(\mathrm{sinh}x)\right\},\\ {\psi }_{2} & = & {\pi }^{-1/4}\sqrt{\displaystyle \frac{{\rm{\Gamma }}({\kappa }_{2}+1/2)}{{\rm{\Gamma }}({\kappa }_{2})}}{({\rm{sech}}x)}^{{\kappa }_{2}}\mathrm{exp}\left\{{\rm{i}}\displaystyle \frac{2{\kappa }_{1}+1}{2}\mathrm{arctan}(\mathrm{sinh}x)\right\},\end{eqnarray} \tag{ 27 }$$

where ${\rm{\Gamma }}(s)$ is the gamma function. (For the complete list of eigenvalues and eigenfunctions, see [22].)

When $q\lt \tfrac{1}{4}$, there is a single eigenvalue ${\kappa }_{1}$. As γ grows from 0 to ${\gamma }_{0}$, where ${\gamma }_{0}=\sqrt{q}$, the branch ${\kappa }_{1}$ descends, monotonically, from $\sqrt{{\gamma }_{c}}-\tfrac{1}{2}\gt 0$ to zero, where the eigenvalue immerses in the continuous spectrum.

When $q\gt \tfrac{1}{4}$, we have two branches. As γ changes from zero to ${\gamma }_{c}$, the branch ${\kappa }_{1}$ descends from $\sqrt{{\gamma }_{c}}-1/2$ to ${\kappa }_{c}$, where ${\kappa }_{c}=\tfrac{1}{2}(\sqrt{2{\gamma }_{c}}-1)\gt 0$. In the interval $({\gamma }_{0},{\gamma }_{c})$, there also is a monotonically growing branch, ${\kappa }_{2}(\gamma )$. As γ varies from ${\gamma }_{0}$ to ${\gamma }_{c}$, the eigenvalue ${\kappa }_{2}$ grows from zero and merges with ${\kappa }_{1}$.

Substituting (27) into (6) gives the flux across the gain–loss interface, associated with the two eigenfunctions:

$$\begin{eqnarray}&&{J}_{1}(\gamma )=\displaystyle \frac{2{\kappa }_{2}+1}{\sqrt{\pi }}\displaystyle \frac{{\rm{\Gamma }}\left({\kappa }_{1}+\tfrac{1}{2}\right)}{{\rm{\Gamma }}({\kappa }_{1})},\quad {J}_{2}(\gamma )=\displaystyle \frac{2{\kappa }_{1}+1}{\sqrt{\pi }}\displaystyle \frac{{\rm{\Gamma }}\left({\kappa }_{2}+\tfrac{1}{2}\right)}{{\rm{\Gamma }}({\kappa }_{2})}.\end{eqnarray} \tag{ 28 }$$

The dependencies (28) are shown in the bottom row of figure 4. In particular, when $q\lt \tfrac{1}{4}$, we have

$$\begin{eqnarray*}&&\displaystyle \frac{{{\rm{d}}{J}}_{1}}{{\rm{d}}\gamma }\to \displaystyle \frac{q}{q-1/4}\lt 0\ {\rm{as}}\ {\kappa }_{1}\to 0.\end{eqnarray*}$$

As in the double-delta potential (figure 2), the behaviour of the flux is anomalous near the point $\gamma ={\gamma }_{0}$ where the eigenvalue immerses in the continuous spectrum.

Zoom In Zoom Out Reset image size

The double-delta well and the Scarff potential exhibit a similar behaviour of the flux near the exceptional point; namely, one branch of eigenvalues has $J(\gamma )$ approaching $J({\gamma }_{c})$ at the infinite negative and the other one at the infinite positive slope. Our objective here is to show that this square-root behaviour is common for all ${ \mathcal P }{ \mathcal T }$-symmetric potentials with exceptional points. (In particular, this implies that an exceptional point cannot be a cusp of $J(\gamma )$.)

Consider an eigenvalue problem

$$\begin{eqnarray}&&H\psi =E\psi ,\end{eqnarray} \tag{ 29 }$$

for the ${ \mathcal P }{ \mathcal T }$-symmetric Schrödinger operator

$$\begin{eqnarray}&&H=-\displaystyle \frac{{{\rm{d}}}^{2}}{{{\rm{d}}{x}}^{2}}+U(x)+{\rm{i}}\gamma { \mathcal W }(x),\end{eqnarray} \tag{ 30 }$$

where the functions $U(x)$ and ${ \mathcal W }(x)$ are real, decay to zero as $| x| \to \infty$ and satisfy $U(-x)=U(x)$ and ${ \mathcal W }(-x)=-{ \mathcal W }(x)$.

Let the operator (30) have two real eigenvalues, ${E}^{(a)}$ and ${E}^{(b)},{E}^{(a,b)}\lt 0$, for all $\gamma \lt {\gamma }_{c}$ in some neighbourhood of ${\gamma }_{c}$, and denote the corresponding eigenfunctions ${\psi }^{(a)}(x)$ and ${\psi }^{(b)}(x)$. Assume that the two eigenfunctions coalesce as γ approaches ${\gamma }_{c}$:

$$\begin{eqnarray*}&&{\psi }^{(a)},{\psi }^{(b)}\to {\psi }_{0},\quad {E}^{(a)},{E}^{(b)}\to {E}_{0}\quad {\rm{as}}\ \gamma \to {\gamma }_{c}.\end{eqnarray*}$$

Therefore, the value $\gamma ={\gamma }_{c}$ is an exceptional point of the operator H, while E₀ is a repeated eigenvalue of algebraic multiplicity two:

$$\begin{eqnarray}&&{H}_{0}{\psi }_{0}={E}_{0}{\psi }_{0},\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&({H}_{0}-{E}_{0}){\phi }_{0}={\psi }_{0}.\end{eqnarray} \tag{ 32 }$$

Here

$$\begin{eqnarray*}&&{H}_{0}=-\displaystyle \frac{{{\rm{d}}}^{2}}{{{\rm{d}}{x}}^{2}}+U(x)+{\rm{i}}{\gamma }_{c}{ \mathcal W }(x);\end{eqnarray*}$$

${\psi }_{0}(x)$ is the eigenvector associated with the eigenvalue E₀ and ${\phi }_{0}(x)$ is the generalised eigenvector.

Without loss of generality we can choose ${\psi }_{0}$ to be ${ \mathcal P }{ \mathcal T }$-symmetric:

$$\begin{eqnarray*}&&{\psi }_{0}^{* }(-x)={\psi }_{0}(x).\end{eqnarray*}$$

(Indeed, if $\psi (x)$ is an eigenvector associated with a real eigenvalue E, then so is the sum $\psi (x)+{\psi }^{* }(-x)$. This sum defines a ${ \mathcal P }{ \mathcal T }$-symmetric function.) Furthermore, if ${\phi }_{0}(x)$ is a solution to equation (32) with a ${ \mathcal P }{ \mathcal T }$-symmetric right-hand side ${\psi }_{0}$, then so is ${\phi }_{0}(x)+{\phi }_{0}^{* }(-x)$. Therefore ${\phi }_{0}$ can also be chosen to be ${ \mathcal P }{ \mathcal T }$-symmetric:

$$\begin{eqnarray*}&&{\phi }_{0}^{* }(-x)={\phi }_{0}(x).\end{eqnarray*}$$

The ${ \mathcal P }{ \mathcal T }$-symmetric solution ${\phi }_{0}$ of equation (32) is defined up to the addition of an arbitrary real multiple of ${\psi }_{0}$. It is convenient to choose this multiple in such a way that

$$\begin{eqnarray}&&\int {\psi }_{0}^{* }{\phi }_{0}\;{\rm{d}}{x}=0.\end{eqnarray} \tag{ 33 }$$

We assume ${\psi }_{0}$ to be normalised to unity,

$$\begin{eqnarray}&&\int {\psi }_{0}^{* }{\psi }_{0}{\rm{d}}{x}=1,\end{eqnarray} \tag{ 34 }$$

and note the following identity which follows from equation (32):

$$\begin{eqnarray}&&\int {\psi }_{0}^{2}\;{\rm{d}}{x}=0.\end{eqnarray} \tag{ 35 }$$

Letting

$$\begin{eqnarray}&&{\epsilon }^{2}={\gamma }_{c}-\gamma \gt 0,\end{eqnarray} \tag{ 36 }$$

we expand the eigenvalue and eigenfunction in powers of $\epsilon$:

$$\begin{eqnarray}E & = & {E}_{0}+\epsilon {E}_{1}+{\epsilon }^{2}{E}_{2}\;+\;\ldots ,\\ \psi (x) & = & {\psi }_{0}(x)+\epsilon {\psi }_{1}(x)+{\epsilon }^{2}{\psi }_{2}(x)+\ldots .\end{eqnarray} \tag{ 37 }$$

Substituting these expansions in (29), we equate coefficients of like powers of $\epsilon$. At the order ${\epsilon }^{1}$, we have

$$\begin{eqnarray*}&&({H}_{0}-{E}_{0}){\psi }_{1}={E}_{1}{\psi }_{0}.\end{eqnarray*}$$

Using equation (32) we obtain a solution:

$$\begin{eqnarray}&&{\psi }_{1}(x)={E}_{1}{\phi }_{0}(x)+\eta {\psi }_{0}(x),\end{eqnarray} \tag{ 38 }$$

where η is an arbitrary constant coefficient. In order to ensure that ${\psi }_{1}$ is ${ \mathcal P }{ \mathcal T }$-symmetric, this coefficient has to be real.

The order ${\epsilon }^{2}$ gives

$$\begin{eqnarray*}&&({H}_{0}-{E}_{0}){\psi }_{2}-{\rm{i}}{ \mathcal W }{\psi }_{0}={E}_{2}{\psi }_{0}+{E}_{1}{\psi }_{1}.\end{eqnarray*}$$

The solvability condition for ${\psi }_{2}$ gives

$$\begin{eqnarray}&&{E}_{1}^{2}=\displaystyle \frac{1}{{\rm{i}}}\displaystyle \frac{\int { \mathcal W }{\psi }_{0}^{2}\;{\rm{d}}{x}}{\int {\psi }_{0}{\phi }_{0}\;{\rm{d}}{x}},\end{eqnarray} \tag{ 39 }$$

where we have made use of (35) and (38).

Using the symmetry of ${\psi }_{0}(x),{\phi }_{0}(x)$ and ${ \mathcal W }(x)$, one can readily check that the quotient (39) is real. The energies ${E}^{(a)}$ and ${E}^{(b)}$ are given by ${E}_{0}+\epsilon {E}_{1}$ and ${E}_{0}-\epsilon {E}_{1}$, respectively, where E₁ is one of the two opposite roots in (39). By assumption, ${E}^{(a)}$ and ${E}^{(b)}$ are real; hence the quotient (39) should be positive. (The negative quotient (39) would simply imply that the real eigenvalues ${E}^{(a,b)}$ pertain to the negative, not positive, parameter ${\epsilon }^{2}$ in (36).)

For the purposes of the flux analysis, it is essential that the eigenfunction $\psi (x)$ be normalised to unity:

$$\begin{eqnarray}&&\int {\psi }^{* }\psi \;{\rm{d}}{x}=1.\end{eqnarray} \tag{ 40 }$$

Substituting from (37) gives

$$\begin{eqnarray*}&&\int {\psi }^{* }\psi \;{\rm{d}}{x}=1+2\eta \epsilon +O({\epsilon }^{2}),\end{eqnarray*}$$

where we have used the identity (33). Thus if we choose $\eta =0$, we will ensure the normalisation condition (40) to order $\epsilon$.

Substituting the expansion of the eigenfunction (37) in (6) we obtain the dependence of the interfacial flux on γ in the vicinity of the exceptional point:

$$\begin{eqnarray*}&&{J}^{(a,b)}={J}_{0}\pm {({\gamma }_{c}-\gamma )}^{1/2}{J}_{1}+O({\gamma }_{c}-\gamma ),\end{eqnarray*}$$

where

$$\begin{eqnarray}{J}_{0} & = & {\left.{\rm{i}}\left(\displaystyle \frac{{\rm{d}}{\psi }_{0}^{* }}{{\rm{d}}{x}}{\psi }_{0}-{\psi }_{0}^{* }\displaystyle \frac{{\rm{d}}{\psi }_{0}}{{\rm{d}}{x}}\right)\right|}_{x=0},\\ {J}_{1} & = & {\rm{i}}{{E}}_{1}{\left[\left(\displaystyle \frac{{\rm{d}}{\phi }_{0}^{* }}{{\rm{d}}{x}}{\psi }_{0}+\displaystyle \frac{{\rm{d}}{\psi }_{0}^{* }}{{\rm{d}}{x}}{\phi }_{0}\right)-{\rm{c.c.}}\right]}_{x=0}.\end{eqnarray} \tag{ 41 }$$

(In (41), the c.c. stands for the complex conjugate of the preceding term.)

The curves ${J}^{(a)}(\gamma )$ and ${J}^{(b)}(\gamma )$ have opposite slopes in the vicinity of ${\gamma }_{c}$:

$$\begin{eqnarray*}&&\displaystyle \frac{{{\rm{d}}{J}}^{(a)}}{{\rm{d}}\gamma }=-\displaystyle \frac{{J}_{1}}{2\sqrt{{\gamma }_{c}-\gamma }}+O(1),\quad \displaystyle \frac{{{\rm{d}}{J}}^{(b)}}{{\rm{d}}\gamma }=\displaystyle \frac{{J}_{1}}{2\sqrt{{\gamma }_{c}-\gamma }}+O(1).\end{eqnarray*}$$

One of the two dependencies is anomalous, ${\rm{d}}{J}/{\rm{d}}\gamma \lt 0$. This is exemplified by figure 2 (second, third and forth panels in the bottom row) and figure 4 (bottom right panel).

The upshot of this analysis is that there is always a jammed branch near the exceptional point.

Does every ${ \mathcal P }{ \mathcal T }$-symmetric potential exhibit an interval of the gain and loss amplitude with the anomalous behaviour of the interfacial flux? In this section we produce an exactly solvable counter-example.

More specifically, we consider a ${ \mathcal P }{ \mathcal T }$-symmetric harmonic oscillator of the form

$$\begin{eqnarray}&&V(x)={x}^{2}-2{\rm{i}}\gamma x,\quad \gamma \gt 0.\end{eqnarray} \tag{ 42 }$$

The Schrödinger operator (4) with this potential supports an infinite sequence of real eigenvalues $-{\kappa }_{n}^{2}=2n+1+{\gamma }^{2}$, where $n=0,1,2,\ldots$ [11, 23]. The corresponding L²-normalized eigenvectors are given by [15]

$$\begin{eqnarray}&&{\psi }_{n}(x)=\displaystyle \frac{{\pi }^{-1/4}}{\sqrt{{2}^{n}n!}}\displaystyle \frac{{{\rm{e}}}^{-{\gamma }^{2}/2}}{\sqrt{{L}_{n}(-2{\gamma }^{2})}}{H}_{n}(x-{\rm{i}}\gamma ){{\rm{e}}}^{-{(x-{\rm{i}}\gamma )}^{2}/2}.\end{eqnarray} \tag{ 43 }$$

In (43), ${H}_{n}(x)$ and ${L}_{n}(x)$ are the n-th order Hermite and Laguerre polynomials, respectively.

Substituting (43) in (6) one can determine the value of the flux associated with each eigenvector. The first three values are given by

$$\begin{eqnarray}{J}_{0}(\gamma ) & = & \displaystyle \frac{2}{\sqrt{\pi }}\;\gamma ,\quad {J}_{1}(\gamma )=\displaystyle \frac{4}{\sqrt{\pi }}\;\displaystyle \frac{\gamma ({\gamma }^{2}+1)}{2{\gamma }^{2}+1},\\ {J}_{2}(\gamma ) & = & \displaystyle \frac{1}{\sqrt{\pi }}\;\displaystyle \frac{\gamma (4{\gamma }^{4}+12{\gamma }^{2}+5)}{2{\gamma }^{4}+4{\gamma }^{2}+1}.\end{eqnarray} \tag{ 44 }$$

As one can readily check, each of the three expressions (44) defines a monotonically growing function of γ. Hence no jamming behaviour is exhibited by the ${ \mathcal P }{ \mathcal T }$-symmetric harmonic oscillator.

This conclusion is consistent with the results of the previous section where we have identified exceptional points as one source of the jamming anomalies. In the system at hand, the ${ \mathcal P }{ \mathcal T }$-symmetry remains unbroken for an arbitrarily large γ and the spectrum does not feature any exceptional points.

Eigenvalue immersions in the continuous spectrum were also seen to be accompanied by jamming (figures 2, 3, and 4). The operator (4) with the parabolic potential (42) does not have any continuous spectrum; hence this source of anomalous behaviour was not available to the oscillator either.

The purpose of this paper was to describe an anomalous phenomenon occurring in ${ \mathcal P }{ \mathcal T }$-symmetric systems of optics and atomic physics. The counter-intuitive effect consists in the reduction of the power gain (or the particle influx) in the active part of the system accompanied by the reduction of the power loss (or particle leakage) in its dissipative part as a result of the increase of the gain–loss coefficient. The reduction of the power gain and loss in their respective parts of the system manifests itself in the drop of the flux between the two parts.

We have been referring to this phenomenon as jamming because of the analogy it bears to the behaviour of the traffic flow through a road network. Here, we explain the traffic analogy in some more detail.

8.1. Jamming in a parallel road network

8.2. Summary and concluding remarks

We have demonstrated that the jamming anomaly in the ${ \mathcal P }{ \mathcal T }$-symmetric linear Schrödinger equation may be anticipated in two broad classes of situations. First, the jamming occurs in the vicinity of the exceptional point, where two real eigenvalues coalesce and acquire imaginary parts. This general conclusion was exemplified by the double-delta well potential with a large distance between the wells ($L\gt 1$) and the deep single-well Scarff potential ($q\gt 1/4$). In contrast, the ${ \mathcal P }{ \mathcal T }$-symmetric harmonic oscillator—the potential that does not support any exceptional points—was shown to be anomaly-free.

Second, the flux associated with an eigenfunction $\psi (x)$ decreases as the corresponding eigenvalue approaches the edge of the continuous spectrum. This law admits a simple explanation. As $\kappa \to 0$, the effective width of the corresponding eigenfunction grows as $1/\kappa$. Since the L²-norm of the eigenfunction $\psi (x)$ remains equal to 1, this implies that $| \psi (0)|$, the maximum value of the single-hump eigenfunction tends to zero (in proportion to ${\kappa }^{1/2}$). As a result, the flux decreases: $J(\gamma )-J({\gamma }_{0})\sim \kappa$.

Thus, if it is the increase of γ that drives the eigenvalue to the continuum, the slope of the $J(\gamma )$ curve near the immersion point will be anomalous. The examples are the double-delta well potential with $L\lt 1$ and the Scarff potential with $q\lt 1/4$.

We have also demonstrated that the jamming anomaly can occur in nonlinear systems. Specifically, we computed a branch of nonlinear modes supported by the ${ \mathcal P }{ \mathcal T }$-symmetric double-delta potential and observed the nonmonotonic dependence of the flux on the gain–loss amplitude. It is worth noting that the nonlinear modes exhibiting the anomalous behaviour can be dynamically stable. The corresponding physical regimes are robust and should be experimentally detectable.

We conclude this section with two remarks. The first one is on the similarity and differences between the jamming anomaly and the macroscopic Zeno effect.

The macroscopic Zeno effect is the descendant of its celebrated quantum namesake [24]. It consists in the drop of the external current employed to compensate losses in the boson condensate— as a result of boosting the damping rate at some local sites [25, 26]. Experimentally, the effect is manifested in the reduction of the decay rate of atoms as the strength of the localised dissipation is increased [27]. Similar to the jamming anomaly, an increase of the dissipation coefficient produces a drop in losses here. However, there are notable differences as well. In particular, raising the damping coefficient in the ${ \mathcal P }{ \mathcal T }$-symmetric system requires the simultaneous and symmetric increase of its gain factor.

Since the light propagating in coupled waveguides obeys the same system of damped nonlinear Schrödinger equations as the boson condensate loaded in the double-well trap, the macroscopic Zeno effect has an optical analogue [28]. Its essential feature is the suppression of the light absorption in the waveguide coupler as the dissipation coefficient in one of its arms exceeds a critical value. A closely related phenomenon was observed in coupled silica micro-toroid resonators [29]. When the parameters of the two-resonator structure were chosen in the vicinity of its spectral exceptional point, the threshold of the Raman lasing was seen to be lowered by raising the dissipation coefficient.

The macroscopic Zeno effect, in condensates and in optics, as well as the enhancement of lasing in microresonators [29], are similar to the jamming anomaly in that the increase of the dissipation coefficient leads to a decrease of losses. One dissimilarity between the condensate traps [25, 26] and optical couplers [28, 29] on the one hand, and our ${ \mathcal P }{ \mathcal T }$-symmetric system on the other, is that the former are purely dissipative set-ups that rely upon the energy or particle influx 'from infinity' (i.e. from outside). In contrast, the ${ \mathcal P }{ \mathcal T }$-symmetric system balances its losses with an internal gain. Jamming is an internal property of the bicentric gain–loss configuration stemming from its short-range structure.

Another difference is that the Zeno effect and the lasing enhancement can be understood using just two modes, whereas the jamming anomaly is a collective phenomenon which requires both active and lossy subsystems to have an internal structure. There is no jamming in the two-mode system; see section 2.2.

As for the second remark, we note that the jamming anomaly provides a simple checking mechanism for the energy captured in the bicentric structure. Consider, for instance, the double-well potential (12) and assume $L\lt 1$. The quantities $| \psi (\pm L/2){| }^{2}$ give the intensity of light in the waveguides with loss and gain, respectively. Using equation (7), these are expressible as

$$\begin{eqnarray*}&&| \psi (\pm L/2){| }^{2}=\displaystyle \frac{J(\gamma )}{2\gamma }.\end{eqnarray*}$$

In the range of gain and loss where ${dJ}/d\gamma \leqslant 0$, that is, in the range ${\gamma }_{* }\leqslant \gamma \lt {\gamma }_{0}$, the intensity is bounded by its value at the point where $J(\gamma )$ is maximum:

$$\begin{eqnarray*}&&| \psi (\pm L/2){| }^{2}\leqslant \displaystyle \frac{J({\gamma }_{* })}{2{\gamma }_{* }}.\end{eqnarray*}$$

We thank Nora Alexeeva for computational assistance and Boris Malomed for useful discussions. This work was funded by the NRF of South Africa (grants UID 85751, 86991, and 87814), the National Institute for Theoretical Physics, and the FCT of Portugal through the grant PTDC/FIS-OPT/1918/2012. The financial support from the UCT Open Access Journal Publications Fund is also gratefully acknowledged.