Control of power in parity-time symmetric lattices

We investigate wave transport properties of parity–time (PT) symmetric lattices that are periodically modulated along the direction of propagation. We demonstrate that in the regime of unbroken PT-symmetry, the system Floquet–Bloch modes may interfere constructively leading to either controlled oscillations or power absorption and unlimited amplification occurring exactly at the phase-transition point. The differential power response is affected by the overlap of the gain and loss system distribution with wave intensity pattern that is formed through Rabi oscillations engaging the coupled Floquet–Bloch modes.


Introduction
A special class of non-Hermitian Hamiltonians that are invariant under simultaneous parity and time reversal and possess a purely real eigenvalue spectrum was initially investigated in the context of quantum mechanics [1]. The PT-symmetry condition requires the real part of the complex potential in the Schrödinger equation to be symmetric, and its imaginary part representing the gain-loss distribution to be antisymmetric. These systems undergo a PT-symmetry breaking transition at some critical peak value of the imaginary part of the potential. Above this threshold value, the spectrum ceases to be real, implying instability for the quantum system. The idea of PT-symmetry has been applied to many areas such as cavity quantum electrodynamics [2], classical mechanics [3], and magnetohydrodynamics [4,5] yet perhaps the most intriguing features of PT-symmetry are being investigated in the context of optics [6][7][8][9][10][11][12][13][14][15][16][17]. Optical systems with a PT-invariant complex refractive index supporting stable propagation of waves were predicted theoretically [6,7] and demonstrated experimentally [8,9]. The features of the linear regime include double refraction, power oscillations, nonreciprocal diffraction [6][7][8][9], amplification of the Goos-Hänchen effect [10], and unidirectional invisibility [11][12][13]. On the other hand, implementation of unidirectional dynamics [14,15], nonreciprocal soliton scattering [16], and switching [17] require combined action of PT-symmetry and nonlinearity. Experimental studies of active PT-symmetric electric circuits [18] and theoretical investigations of PT-symmetric nonlinear magnetic metamaterials [19] were also reported recently.
Longitudinally modulated PT-symmetric structures, which have been attracting increasing attention, exhibit even richer behavior and potential for a broad range of applications. Recent investigations in this direction include theoretical studies of wave dynamics in complex crystals with PT-symmetry subjected to a sinusoidal ac force [20] that predicted dynamic localization of a wave packet due to the band collapse in the regime of unbroken PT-symmetry and nonreciprocal Bragg scattering in the vicinity of the breaking point. In [21], the authors introduce the concept of local PT-symmetric invariance exemplified by two coupled waveguides whose longitudinally modulated gain-loss profile remains antisymmetric at any propagation distance. Depending on the initial conditions, such systems can support symmetric or antisymmetric power evolution and unidirectional phase exchange between waveguides. The phase-transition in PT-symmetric tight-binding lattices with periodic longitudinal modulations of the gain-loss distribution was shown to depend on the ratio between frequency and amplitude of the modulations [22]. In the unbroken PT regime, these systems exhibited hyperballistic transport properties, while in the vicinity of phase-transition, they supported diffraction-free propagation. A stochastic PT-symmetric coupler with fluctuating parameters, such that the gains and the losses are exactly balanced on average, was investigated in [23]. It was shown that even in the regime of unbroken PT symmetry, the statistically averaged intensity of the wave in such systems grows. Nonlinear properties of longitudinally modulated PT-symmetric systems (dual-core nonlinear waveguides) include stabilization of the solitons by application of management, specifically periodic simultaneous switching of the sign of the gain, loss, and inter-core coupling [24]. Recently it was also shown that several classes of longitudinally modulated non-PT-symmetric waveguides can possess an all-real spectrum [25,26].
Among other effects, optical Rabi oscillations, i.e., resonant power transitions between different light modes originally proposed [27,28] and observed [29] in Hermitian waveguides and wave-guide arrays, were also predicted and periodically modulated in the direction of propagation PT-symmetric waveguides [30]. In the present work, we focus directly on these phenomena, specifically on propagation dynamics in complex PTsymmetric wave-guide arrays periodically modulated along the direction of propagation. We show that in the regime of unbroken PT-symmetry when all Floquet-Bloch (FB) modes remain stable, mode interference results in a wealth of behaviors including damping, amplification, or even unlimited amplification of incident beam power distributed linearly or quadratically in the propagation direction. The different power regimes depend on the resonant frequency overlap of the gain-loss distribution with the intensity pattern formed through the interference of Rabi-coupled FB modes. A noteworthy feature is the appearance of an intensity-locked region at the phase-transition point that is bounded by two beams formed through the double refraction of the unique incident beam.

Model
Paraxial beam propagation in the wave-guide arrays periodically modulated along the direction of propagation can be described through the Schrödinger equation: where y is the complex amplitude of the beam, x and z are transverse and longitudinal coordinates, respectively, the amplitude of longitudinal modulations e is taken to be small, i.e., e  1,b 0 is frequency of longitudinal modulations, and the complex PT-symmetric potential is given by where D is the lattice period. The real part of U that describes the refractive index profile is symmetric in the transverse direction while the imaginary part representing the gain-loss distribution is antisymmetric, i.e., fulfills the PT-symmetry condition.
The PT-symmetry of the Hamiltonian is necessary but not a sufficient condition for transmission spectrum to be real. PT-symmetry is broken when imaginary amplitude of the refractive index W becomes higher than its real amplitude V [6,7]. In the broken PT-symmetry regime, propagation constants merge together forming pairs of complex conjugate values and thus making the wave propagation unstable. In what follows, we consider two regimes of the wave propagation. In the first (under-critical) regime < W V , whereas in the second (phasetransition) regime = W V . Both considered regimes belong to the unbroken (stable) phase when separate FB modes exhibit stable propagation.

Coupled mode theory
The total field distribution in the wave-guide arrays can be approximated by a superposition of normalized FB modes, i.e., , e x pi d , The latter expression was obtained utilizing the bi-orthogonality condition [6,7] * ò j Equations (4) and (5) were derived with an additional condition imposed on eigenvectors 6 : For the regime of unbroken PT-symmetry (W/V=0.8), the first three bands of the transmission spectra are depicted in figure 1(a). As in the case of real periodic potential [28], longitudinal modulation can support Rabi conversions between two different FB modes provided the resonant condition is satisfied; i.e., the frequency of longitudinal modulations b 0 must be matched to the difference of propagation constants of the corresponding FB modes. The transitions between first and second band modes for the right and left incidences (positive and negative Bloch wave numbers, respectively) are schematically shown by green arrows in figure 1(a). Following the standard procedure [28], one can substitute expansion (2) into equation (1), apply bi-orthogonality condition (5) (bi-orthogonality should be used because skewed FB modes of PT-symmetric lattices are not orthogonal), and, after neglecting off-resonant terms, obtain coupled mode equations for the evolution of mode population coefficients where the coupling coefficients given by the overlap integral remain purely real and positive if condition (6) is satisfied.

Propagation dynamics in the under-critical regime
The coupling coefficients for 1→2 transitions in the under-critical (V=1 and W=0.8) regime are shown in figure 1 The striking difference with the case of Hermitian (real) Hamiltonians is asymmetry of the coupling coefficients with respect to k, i.e.,  (7), whereas the dark-colored thin lines were obtained from equation (4) with field distribution y ( ) which is identical to the frequency of Rabi oscillations in real periodic potential [28]. However, the asymmetry of coupling coefficients results in unequal amplitudes of mode population oscillations which ratio is given by 21 12 The physical reason for different amplitudes is that the exchange of energy between two FB modes is augmented by their amplification (damping) due to the positive (negative) overlap of the intensity pattern with the gain-loss distribution. The system can experience overall amplification (damping) if the integral of intensityweighted gain-loss distribution (distribution of imaginary part of refractive index) is positive (negative). For example, in the case of the left tilted incidence > ( ) M M 21 12 in the first half-period, the second band mode accepts power from the first band mode and is also amplified by the positive overlap, whereas in the second halfperiod it loses its power to the first band mode yet it is also dampened by negative overlap. The first band mode in 6 In [6, 7] expression for modes population and bi-orthogonality condition include the coefficients d kn =±1 Our corresponding expressions (equations (4, 5)) do not include these coefficients due to additional condition given by equation (6).
the first half-period loses its power to the second band mode yet it is also amplified by the positive overlap; whereas in the second half-period it accepts power from the second band mode but is also damped by the negative overlap. The evolution of the total power is given by where highly oscillating terms resulting from the non-orthogonality of FB modes were dropped, amplitude of power oscillations is given by 21 21 12 and ò ò 2 / / / / represents the ratio of power contents of different modes. In contrast with the case of real potential, this ratio is not necessarily equal to one because FB modes were normalized using the bi-orthogonality condition (5). The power evolution given by equation (9)  In this case, both modes and total power are amplified during the first half-period and dampened in the second half-period. For the right tilted incidence h < 0 21 so dynamics is opposite: both modes and total power are dampened during the first half-period and amplified in the second half-period. The blue curves in figures 1(e) and (f) represent power evolution obtained by direct numerical solution of equation (1). Numerically obtained power exhibits fast oscillations around an analytical solution represented by equation (9) due to the nonorthogonality of FB modes [6,7]. These fast oscillations diminish in the vicinity of points where population of one of the modes becomes small. The dynamics of the 1→2 transitions that are described previously can be also observed in figures 1(g) and (h) depicting numerically found intensity evolution patterns for the cases of the left and right tilted incidences, respectively. In the case of normal incidence = M 0 21 and = M 0 12 therefore the system cannot support 1→2 transitions.
The coupling coefficients for 1→3 transitions are shown in figure 2(a). The dynamics in the case of 1→3 transitions are very similar for the tilted incidences and therefore are not shown here. The only difference with respect to the 1→2 transitions is that q > M M 1 31 31 13 for the right tilted incidence and q < M M 1 31 31 13 for the left tilted incidence. Therefore, the power is initially increasing for the right tilted incidence and vice versa. However for the normal incidence ( = ) k 0 the dynamics of 1→3 transitions is different. The transition between first and third band modes for the normal incidence is schematically shown by the red arrow in . Evolution of power for the 1→3 transition governed by equation (9) with index 2 replaced by 3 is shown in figure 2(c) by the red curve. Power oscillations result from unequal power content of different modes, specifically q ¹ 1. 31 The blue curve in figure 2(c) represents power evolution obtained by direct numerical solution of equation (1). The fast oscillations of the numerical solution result from non-orthogonality of FB modes. The dynamics of the 1→3 transitions that are described above can be also observed in figure 2(d) depicting numerically found intensity evolution pattern for the case of normal incidence.

Propagation dynamics in the phase-transition regime
The most intriguing feature involves the propagation in the phase-transition regime when W/V=1. The first three bands of the transmission spectra in this regime (V=W=1) are depicted in figure 3(a) and the coupling and third (green curves) band mode populations for the normal incidence; light-colored thick lines denote analytical results, dark-colored thin lines denote numerical results; (c) evolution of power normalized to input power for the normal incidence; red curve denotes analytical results, blue curve denotes numerical results; (d) numerically found intensity evolution pattern for the normal incidence. coefficients for 1→2 transitions are shown in figure 3(b). One can observe that It follows then from equation (7) that for the left tilted incidence = A const,  (c) evolution of the first (blue curves) and second (green curves) band mode populations; light-colored thick lines-analytical results, dark-colored thin lines-numerical results; (d) evolution of power normalized to initial power; red curves-analytical results, blue curves-numerical results; (e)-(h) input: Gaussian beam ψ(x, z=0)=exp(−x 2 /σ 2 )exp(ik 0 x), k 0 D/2π=−0.15, σ=60; (e) evolution of the first (blue curves) and second (green curves) band mode populations; light-colored thick lines-analytical results, dark-colored thin lines-numerical results; (f) numerically found intensity evolution pattern; (g) numerical (blue curve) and analytical (red curve) output (z=400) intensity profile; inset shows zoom area; (h) evolution of power normalized to initial power; red curves-analytical results, blue curves-numerical results.
In the derivation of equation (12) we assumed that the resonance condition for first and second band modes respectively. Dark-colored thin lines were obtained by substituting numerically found field distribution y ( ) x z , into equation (4). The second band mode is amplified at a constant rate due to the constant population of the first band mode and positive overlap of the intensity pattern with gain-loss distribution. The total power, averaged over the fast oscillations, grows quadratically according to where V q = M 4. 21

2
The power evolution given by equation (13) is depicted by the red curve in figure 3(d). The blue curve in figure 3(d) represents numerically found power evolution. The fast oscillations of numerical solution result from the beating of non-orthogonal FB modes.
In the following we will focus on the evolution of the Gaussian incident beam which can be represented as superposition of modes with different Bloch wave numbers. The propagation dynamics of the left tilted incident beam y , is shown in figures 3(e)-(h). The evolution of the first (blue curves) and second (green curves) band mode populations at = k k 0 are depicted in figure 3(e). The light-colored thick lines represent analytical results given by = A const 1 and equation (11) for first and second band modes respectively. Dark-colored thin lines were obtained by substituting numerically found field distribution ψ(x, z) into equation (4).
The population of the first band mode remains constant upon propagation and the population of the second band mode grows almost linearly (a small deviation of linear growth at the beginning is due to the ¹ ( ) A k , 0 0 2 0 term of equation (11)). The numerically obtained evolution of intensity pattern is shown in figure 3(f). The video file describing the evolution of the intensity when Gaussian beam is excited at the input is presented in the online supplementary material. We observe the phenomenon of double refraction which manifests itself as splitting of the single incident beam into two diverging beams propagating with the group velocities of the first and second band modes, i.e., b [6,7]. In addition to the double refraction, the intensity pattern exhibits an intensity-locked region between these diverging beams, which develops due to the positive overlap of the intensity pattern with the gain-loss distribution. Intensity oscillations in this region conform to the periodicity of the lattice while the amplitude of oscillations remains constant between the two diverging beams. The amplitude of the intensity oscillation in this region remains constant upon propagation whereas its width grows linearly as the distance between the diverging beams increases. This intensity-locked region corresponds to the first term on the RHS of equation (12). Its intensity (absolute value squared) at the output plane (z=400) is shown in figure 3(g) by the red curve along with numerically found intensity (blue curve).
The evolution of the total power is given by 2 / / and the highly oscillating terms were dropped. The power evolution given by equation (14) is depicted by the red curve in figure 3(h). The blue curve in figure 3(h) represents the numerically found total power of the beam. As previously, the fast oscillations result from the beating of non-orthogonal FB modes while the linear power growth is due to the linearly increasing width of the intensity-locked region.
In the case of right tilted incidence the linear amplification can be implemented through the 1→3 transition since Evolution of the third band mode population, power, and field for the 1→3 transition is governed by equations (10)-(13) with index 2 replaced by 3. In the case of normal incidence (k=0) second and third bands coalesce forming an exceptional point. Corresponding FB states become identical and self-orthogonal at this point [31]. This implies that in order for these states to remain normalized to unite, their amplitudes (components of eigenvectors ) u l 2 should become extremely large as  k 0. Our coupled mode analysis fails at k=0 because self-orthogonality precludes normalization of FB modes. However direct simulation of beam propagation shows that the amplitude of the intensity-locked region remains finite and continuously changes as 

Control of power
Importantly, in the regime of unbroken PT-symmetry (W<V) one can manipulate both the sign and magnitude of power oscillations by changing the angle of incidence and corresponding adjustment of longitudinal frequency modulations to satisfy the phase matching condition. This can be observed in figure 4(a), which depicts amplitudes of power oscillations h 21 and h 31 versus Bloch momentum k for W/V=0.8. However at the phase-transition regime W=V only the rate of always positive power amplification can be manipulated through angular dependence of growth coefficients V 21 or V . 31 Figure 4(b) shows V 21 and V 31 versus k for W=V.

Propagation dynamics in two-dimensional lattices
The phenomena described so far are not particular to one-dimensional lattices. We have also investigated propagation in longitudinally modulated two-dimensional PT-symmetric lattices with potential , 0 x y x y y Figure 5(a) depicts the intensity pattern at the output plane (z=480). The intensity-locked regions in the twodimensional case correspond to two 'ridges' that are directed along the x and y axes. Intensity oscillations in these 'ridges' conform to the periodicity of the lattice while the amplitude of oscillations remains constant between the beams that diverge in the x and y direction from the main beam propagating along the z axis. The flat 'plateau' region between the two 'ridges' develops due to the secondary generation of intensity-locked regions in perpendicular directions: the 'x-ridge' radiates in the y direction and the 'y-ridge' radiates in the x direction. The development of this square 'plateau' region is responsible for quadratic amplification of power shown in figure 5(b).

Conclusions
While power amplification has been predicted in the gain-loss media [32,33], in the present work we have investigated power control phenomena in the complex PT-symmetric wave-guide arrays that are periodically modulated along the direction of propagation. We demonstrated periodic power oscillations as well as linear or quadratic amplification of the power that occurs in the regime of unbroken PT-symmetry at certain resonant frequencies of longitudinal modulation. This differential power response stems from the overlap of the gain-loss distribution with the intensity pattern formed through interference of Rabi-coupled FB modes. Linear amplification that occurs at the phase-transition regime is due to the development of the intensity-locked region between the two diverging beams induced through double refraction of the unique incident beam. However, when a single Bloch mode is excited at the input its power oscillates in the under-critical regime, whereas in the phase-transition regime it grows quadratically. In the phase-transition regime, only the rate of power amplification can be controlled through the angular dependence of the coupling coefficients. On the other hand, both the sign and magnitude of power oscillations that occur in the under-critical regime can be manipulated through the incidence angle of the beam. In both cases, the frequency of longitudinal modulations should be adjusted accordingly to satisfy the phase matching condition. We used coupled mode theory to analyze quantitatively the beam dynamics and verified the results numerically. We found similar effects in the twodimensional longitudinally modulated PT-symmetric lattices. The phase-transition point in the twodimensional case is characterized by the development of perpendicular intensity-locked 'ridges' with a flat 'plateau' region between them and the quadratic amplification of the power. The efficient power control technique introduced here should pave the way for new discoveries in the flourishing field of PT-symmetric optics and may have significant applications in integrated PT photonic devices as well as in gain-loss metamaterials and metasurfaces.