Optical waveguide arrays: quantum effects and PT symmetry breaking

Over the last two decades, advances in fabrication have led to significant progress in creating patterned heterostructures that support either carriers, such as electrons or holes, with specific band structure or electromagnetic waves with a given mode structure and dispersion. In this article, we review the properties of light in coupled optical waveguides that support specific energy spectra, with or without the effects of disorder, that are well-described by a Hermitian tight-binding model. We show that with a judicious choice of the initial wave packet, this system displays the characteristics of a quantum particle, including transverse photonic transport and localization, and that of a classical particle. We extend the analysis to non-Hermitian, parity and time-reversal ($\mathcal{PT}$) symmetric Hamiltonians which physically represent waveguide arrays with spatially separated, balanced absorption or amplification. We show that coupled waveguides are an ideal candidate to simulate $\mathcal{PT}$-symmetric Hamiltonians and the transition from a purely real energy spectrum to a spectrum with complex conjugate eigenvalues that occurs in them.


Introduction
Historically, light and matter have been considered two quintessentially different entities. Since the advent of quantum theory, which elucidates the wave nature of material particles and the particle nature of electromagnetic waves, properties of quantum system of particles are described by a (possibly many-body) wave function whose time evolution is determined by the Schrödinger equation [1,2]. Such many-body condensed matter systems support collective excitations whose energy is linearly proportional to the momentum, and thus allow one to mimic lightlinearly dispersing massless excitations -in material systems [3]. However, due to the unique nature of electromagnetic waves, namely the lack of a rest-frame or, equivalently, zero rest mass, they were not considered useful for simulating the behavior of quantum particles with nonzero rest mass [4].
Over the past decade, the tremendous progress in fabrication and characterization of semiconductor heterostructures has made it possible to create arrays of evanescently coupled optical waveguides with numbers varying from a few to a few hundred [5,6]. The resulting "diffraction management" [7] makes evanescently coupled waveguides a paradigm for the realization of a quantum particle hopping on one or two dimensional lattices, and permits the observation of quantum and condensed matter phenomena in macroscopic samples using electromagnetic waves. One can engineer such a waveguide array to model any desired form of tight-binding, non-interacting Hamiltonian, because the local index of refraction and the width of the waveguide determine the on-site potential for the Hamiltonian while the tunneling amplitude from one site to its adjacent site can be changed by changing the separation between adjacent waveguides [8,9,10]. A variation in the index of refraction or the tunneling amplitude, both of which can be introduced easily, permit the modeling of a tight-binding Hamiltonian with site or bond disorders respectively. Due to this versatility, many quantum and condensed matter phenomena -Bloch oscillations [11], Dirac zitterbewegung [12], and increased intensity fluctuations [13,14] of light undergoing Anderson localization [15,16,17] -have been theoretically predicted to occur or experimentally observed in waveguide arrays. They have been used to investigate solitonic solutions that arise due to nonlinearities in the dielectric response [18,19]. Such arrays of coupled waveguides have also been used to simulate the quantum walks of a single photon [20,21], correlated photons [22], and Hanbury Brown and Twiss (HBT) correlations [23]. Most recently, they have been used to create a "topological insulator", an exotic state of matter in which the bulk is an insulator, but the two surfaces are conductors [24,25].
There are several advantages to using waveguides to investigate quantum behavior and statistics. First, the quantum effects are measurable over much longer distances than those in condensed matter systems with electrons or in cold-atom systems in electromagnetic traps. Second, instead of an indirect measurement through observables such as conductivity or other response functions [26], in optical waveguides, one can directly measure the timeevolution of a wave function via the time-and-space dependent probability distribution, since it is identical to the light intensity distribution. For lattice models realized via electronic or cold-atom systems, typically, eigenstates in a small fraction of the energy band near the Fermi energy are experimentally probed [27]; in contrast, the ability to create an initial wave packet localized to a single waveguide -by coupling light into a single waveguidemeans that quantum effects across the entire energy band of the tight-binding model can be investigated in optical waveguide arrays.
In the past fifteen years, there has been significant theoretical research on properties of non-Hermitian Hamiltonians that, sometimes, show purely real spectra [28,29,30,31]. In continuum models, such Hamiltonians usually consist of a Hermitian kinetic energy term and a complex potential that is invariant under the combined operation of parity and time-reversal (PT ), such as V (x) = x 2 (ix) or V (x) = n R (x) + in I (x) where n R (x) and n I (x) are even and odd functions of x, respectively. The region of parameter space where all energy eigenvalues of a PTsymmetric Hamiltonian are real is traditionally called the PT -symmetric region, and the emergence of complex conjugate eigenvalues that accompanies departure from this region is called PT -symmetry breaking. Since the effective potential in an optical waveguide array is given by the local (complex) index of refraction, properties of PT Hamiltonians have led to predictions of new optical phenomenon such as Bloch oscillations in complex crystals [32], an optical medium that can simultaneously act as an emitter and a perfect absorber of coherent waves [33,34], PTsymmetric Dirac equation [35], and induced quantum coherence between Bose-Einstein condensates [36]. The PTsymmetry breaking has recently been experimentally observed in two coupled waveguides [37,38], silicon photonic circuits [39], and optical networks [40]. Thus, coupled optical waveguide arrays are also an ideal candidate to simulate the quantum dynamics of a non-Hermitian, PTsymmetric Hamiltonian.
In this paper, we review properties of coupled optical waveguides. In the absence of any loss or gain in a waveguide, the effective Hamiltonian of such an array is Hermitian. In Sec. 2 we present the basics of such Hermitian, tight-binding models, and discuss quantum photonic transport (Sec. 2.1), continuum quasiclassical limit (Sec. 2.2), arrays with position-dependent nearest-neighbor tunneling (Sec. 2.3), and the effects of on-site and tunneling disorder (Sec. 2.4). Section 3 focuses on PT -symmetric tight-binding models where the non-Hermitian, PT -symmetric potential corresponds to loss in one waveguide and an equal gain in its mirror-symmetric counterpart waveguide. We introduce the terminology, present the PT -symmetric phase diagram for arrays with open boundary conditions (Sec. 3.1), discuss the salient features of non-unitary time evolution in such systems (Sec. 3.2), and compare the effects of Hermitian vs. non-Hermitian, PT -symmetric disorder on intensity correlations (Sec. 3.3). We conclude this review with a brief discussion of open questions in Sec. 4. Fig. 1. Schematic of an array of evanescently coupled optical waveguides. The height h and the width w of the waveguide determine the spatial profile of electromagnetic modes inside it, along with the effective potential β j in waveguide j, and the distance d between the centers of adjacent waveguides determines the effective tunneling C j+1,j between them. Due to its constant speed, the motion of light along the waveguide is equivalent to its time evolution whereas motion across different waveguides simulates a quantum particle on a tight-binding lattice.

Hermitian Tight-binding Models
The Hamiltonian for a one-dimensional array with N identical, single-mode waveguides is given by (1) where = h/(2π) is the scaled Planck's constant, a † j (a j ) is the Bosonic creation (annihilation) operator for the single mode in waveguide j, β j is the effective potential on site j or equivalently, the propagation constant for waveguide j, and C j+1,j denotes the tunneling amplitude from site j to adjacent site j + 1. Based upon its geometry, the array can have open boundary conditions (C N = 0) or periodic boundary conditions, a † N +1 = a † 1 . It is straightforward to generalize this Hamiltonian to two-dimensional arrays. The on-site potential β j and the tunneling amplitude C j+1,j are determined by profile of the electric field u(r) in a single waveguide as where k 0 is the wavenumber for the incident light, c is the speed of light in vacuum, k j characterizes the wave vector for the single eigenmode in waveguide j, n j+1 and n b are refractive indices for waveguide j + 1 and the barrier between adjacent waveguides respectively, and denotes integral over the two-dimensional cross section of waveguide j. Thus, the potential β j is linearly proportional to the local index of refraction n j , whereas the tunneling amplitude is proportional to the overlap between the electric-field envelope functions in waveguides j and j + 1. Note that it is possible to create a non-Hermitian tunneling profile -C j+1.j = C j,j+1 -by varying the index of refraction and maintaining the waveguide geometry; we will, however, only consider waveguide arrays where the tunneling is Hermitian, C j+1,j = C j,j+1 = C j . The electromagnetic waves in dielectric media do not interact with each other when the light intensity is small and the effects of non-linear susceptibility χ 3 can be ignored [41]; therefore, there are no quartic "interaction terms" in the Hamiltonian. Thus, the Hamiltonian that describes the time-evolution of an electromagnetic pulse (with many photons) in an array of waveguides is equivalent to that of a single quantum particle hopping on a lattice with on-site potentials β j and tunneling amplitudes C j . This absence of interaction allows us to use the coupled waveguide array as an exquisite probe of competition among dispersion, disorder, quantum statistics, and boundary conditions.
When the on-site potential is constant β j = 0 and tunneling amplitudes are constant, C j = C, the Hamiltonian is translationally invariant. Therefore, it can be diagonalized using eigenfunctions ψ kn (j) characterized by eigenmomentum k n . The energy spectrum of the one-dimensional lattice is given by E(k n ) = −∆ B cos(k n )/2 where ∆ B = 4 C is the bandwidth and the dimensionless eigenmomenta are k n = nπ/(N + 1) (n = 1, . . . , N ) for open boundary conditions and k n = ±2nπ/N with n = 0, . . . , N/2 for periodic boundary conditions [42]. It follows then that for an array with N → ∞ sites and lattice spacing a, the permitted dimensionful wave vectors k form a continuum, bounded by −π/a < k ≤ π/a, known as the first Brillouin zone [43].
We emphasize that although electrons in condensed matter materials and light in optical waveguide arrays can both be described by Eq.(1), the relevant lattice-site numbers and energy scales in the two cases are vastly different. For electronic materials, the number of atoms or lattice sites is N 10 9 whereas for light, the number of coupled waveguides is N 100. For electrons, the tunneling amplitude C ∼ 1 eV or equivalently, C ∼ 240 THz and C/(2πc) ∼ 8000 cm −1 , whereas the on-site potential β ∼ E F ∆ B where E F is the Fermi energy; these parameters cannot be varied significantly (by orders of magnitude) since Coulomb interactions are the primary determinant for these parameters. For light, the tunneling amplitude, determined by the distance between adjacent waveguides, is C/(2πc) ∼ 3-50 cm −1 . Thus, the typical bandwidth of the waveguide array, ∆ B ∼ a few meV, is smaller than its electronic counterpart by orders of magnitude. In addition, the on-site potential in waveguide arrays can be comparable with the bandwidth, β ∼ ∆ B ∼ 10-100 cm −1 . This tremendous flexibility, present even in a small array with a constant tunneling, hints at the rich possibilities for designing waveguide arrays with dramatically different properties. In the following subsections, we will illustrate this point with a few examples.

Phase controlled photonic transport
In free space, the change in the momentum of a particle under constant force, or equivalently, a potential that varies linearly with position, is proportional to the time and thus can increase continuously. In sharp contrast, when a particle on a lattice is acted upon by a constant force, its momentum change is bounded by the size of the Brillouin zone. Physically, the particle can transfer its momentum to the underlying lattice as long as the transferred momentum is equal to one of the reciprocal lattice vectors, and therefore, the momentum of the particle is only defined within the bounds of the first Brillouin zone. This surprising result, which occurs only due to the presence of the lattice, implies that the velocity of the particle oscillates about zero in the presence of a constant force, and is called Bloch oscillations. For electronic materials in constant electric field, the time required for the requisite change of momentum is given by t B = ∆p/(qE) ∼ /(qEd) where q is the electronic charge, E is the applied, constant electric field, and d ∼ fewÅ is the lattice constant. For typical fields E ∼ 10 3 V/m, this time is orders of magnitude longer than the typical time t s between electron-lattice scatterings, t B ∼ 10 −8 sec t s ∼ 10 −14 sec [43,44]. Therefore, although long predicted in electronic systems, Bloch oscillations have not been and are unlikely to be observed in them. In addition, due to the large number of lattice sites, the effects of boundary on Bloch oscillations cannot be explored in electronic materials. Since there is no interaction of light with the dielectric and therefore no scattering that can randomize the transverse momentum of a wave packet in a lattice of waveguides, they provide an ideal platform to study Bloch oscillations and other energy-band related quantum phenomena in finite lattices where boundary effects can be prominent [11].
To this end, we consider the waveguide array with a linear ramp in the on-site potential given by β j = β 0 +δβj with δβ/β 0 1. Since β 0 only shifts the zero of the energy spectrum, we will ignore it in the subsequent treatment. This system is created by using variable-width waveguides with variable spacing between them to ensure constant tunneling and a linear gradient with δβ/β 0 ∼ 10 −4 [11]. The equation of motion for the electric-field creation operator is given by i ∂a † j /∂t = [H, a † j ] and reduces to where one of the tunneling terms is absent when the site index j corresponds to the first or the last waveguide in an array with N waveguides. In the limit N → ∞, this equation can be exactly solved by using Fourier transform [14] and we get the following expression for the time-evolution The right-hand column shows the corresponding differences ∆I(p, t) between the exact solution for a finite array and the analytical result for an infinite array. Note that, on average, ∆I(p, t) increases with time, but becomes appreciable only after the ballistically expanding wave packet has reached the boundaries. operator Note that as the potential gradient vanishes, δβ → 0, we recover the propagator for a uniform lattice with bandwidth ∆ B = 4 C. The time-evolution operator allows us to obtain the time and site-dependent intensity for an arbitrary normalized initial state |ψ(0) = m α m a † m (0)|0 , where the sum of weights is unity, m |α m | 2 = 1. If the initial input is confined to a single waveguide, α m = δ m,m0 , the intensity profile becomes This analytical result for the site and time-dependent intensity has the following features: It is symmetrical about the initial wave packet location; it is periodic in time with a period given by T = 2π/δβ; its maximum spread occurs at time t = T /2 and is determined by the ratio of the nearest-neighbor tunneling to the potential gradient C/δβ. We emphasize that this result is valid only for an infinite array where the effects of boundaries can be ignored. On the other hand, since the tight-binding Hamiltonian Eq.(1) for a finite array corresponds to a finite, tri-diagonal, Hermitian matrix, one can obtain the timeevolved wave function p|ψ(t) exactly by straightforward numerical evaluation of the time-evolution operator G(t).
The left-hand panels in Fig. 2 show the numerically obtained intensity I(p, t) for an N = 21 waveguide array with initial input in the central waveguide m 0 = (N + 1)/2 = 12; we use t 0 = /∆ B as the unit of time. For δβ = C (bottom panel), the period of Bloch oscillations is given by T /t 0 = 8πC/δβ = 8π and the maximum spread of intensity is small compared to the size of the array. When δβ = 0.5C (center panel), the period is doubled, T /t 0 = 16π and so is the vertical maximum spread of intensity. For δβ = 0.25C (top panel), the estimated wave packet spread is greater than the size of the array, and the open boundaries destroy Bloch oscillations although the intensity profile continues to remain symmetric about the center of the array. The right-hand panels in Fig. 2 show the difference between numerically obtained intensity profile and the analytical result that is valid only for an infinite array, ∆I(p, t) = I(p, t) − I a (p, t). When δβ = 0.25C (top panel), the wave packet reaches the boundaries and thus the difference between the exact solution and the analytical result is the greatest, although we point out that this difference becomes appreciable only after the ballistically expanding wave packet has reached the array boundaries. When δβ = 0.5C (center panel), the maximum intensity difference is approximately 1% of the total intensity, although it increases with subsequent reflections from the boundaries of the finite array. When δβ = 0.25C (bottom panel), the intensity difference ∆I is essentially zero. Thus, although the analytical result for the site and time dependent intensity is ideally applicable only for an infinite array, it accurately describes the dynamics of a finite array as long as the maximum spread of the wave packet does not detect the array boundaries.
The symmetrical intensity distribution in Bloch oscillations seen in Fig. 2 is because all momenta within the Brillouin zone have equal weight in an input state that is localized to a single site. Next, we consider an initial state that is localized to two adjacent waveguides, with a phase difference φ between the two, α m = cos θδ m0 +sin θe iφ δ m1 . The analytical result for the site-and time-dependent intensity is given by where τ (t) = (4C/δβ) sin(δβt/2) and the last term in the intensity arises as a result of the interference between the two inputs. To quantify this interference, we consider the time-dependent average and standard deviation of the position, which, for an infinite array, can be simplified to Note that when the input is only confined to the central, zeroth waveguide, sin θ = 0, we recover j mean (t) = 0 and j 2 std = τ 2 /2, and when the light is completely confined to the first waveguide, sin θ = 1, we obtain the expected results. At small times, since the function τ (t) ≈ 2Ct, Eqs.(10)- (11) imply that the mean position and its standard deviation both change linearly with time except when φ = {0, π}; in those two cases, they change quadratically with time. At large times, the mean position and standard deviation both oscillate due to the periodic nature of the function τ (t) = τ (t + 4π/δβ). These results are only valid for an infinite array and, as we have seen earlier, they remain applicable to a finite array only if the maximum spread of the wave packet is smaller than the size of the array. Figure 3 shows the effects of the relative phase φ and initial wave packet on the intensity profile I(p, t) (lefthand panels) and the mean position j mean (t) (right-hand panels) for an N = 21 waveguide array with δβ = 0.5C, period T /t 0 = 16π, and equally distributed weight on the adjacent sites, θ = π/4. In the left-hand column, the top panels shows the asymmetrical intensity profile that results from an initially symmetric state |ψ(0) = (|m 0 + |m 0 + 1 )/ √ 2 with m 0 = 11. The center panel shows corresponding intensity profile for |ψ(0) = (|m 0 + i|m 0 + 1 )/ √ 2, with m 0 = 11, where the asymmetry in the intensity profile switches direction with time. Both of these numerically obtained results are virtually identical with those obtained from Eq.(9) that is valid for an infinite array. The bottom panel shows that the same wave function, with m 0 = 16, gives rise to an aperiodic intensity profile due to the presence of the boundary. The righthand panels in Fig. 3 show corresponding mean position of the wave packet. When φ = 0 (top panel), the mean position is confined to the region of lower index of refraction and changes quadratically with time. When φ = π/2 (center panel) we see that j mean (t) oscillates about the initial mean position, and changes linearly with time at small times. The bottom panel shows that when the initial position is close to the boundary, the periodic behavior is destroyed due to the added interference with partial waves that are reflected from one edge of the array. Thus, the direction of the lateral photonic transport can be tuned by the relative phase difference φ between inputs at adjacent waveguides.

Continuum limit: non-relativistic particle
In the last subsection, we considered the time evolution of a wave packet that is initially localized to one or two sites. Due to this extreme localization in real space, such a wave packet has components with all momenta (or equivalently, energies) across the entire bandwidth of the onedimensional lattice. Due to the presence of these dimensionless momenta −π < k ≤ π, the time evolution of the wave packet is dominated by quantum interference. On the other hand, by an appropriate choice of initial state that has energy components only near the bottom or the top of the cosine-band E(k) = −2 C cos(k), one can mimic the behavior of a non-relativistic particle on a line segment.
To formalize this mapping from a lattice to the continuum, let us consider lattice with sites N → ∞ and site-to-site distance d → 0 such that N d → L [42]. We will choose a continuum co-ordinate system such that site m = 1 maps to x = −L/2 whereas site m = N maps to x = +L/2. In this limit, the nearest-neighbor tunneling term in Eq.(4) translates into a spatial second-derivative with effective mass m * given by Therefore, time evolution of an initial state |ψ e with components only near the bottom of the band, k ∼ 0, in the presence of a linearly varying potential V (x) = 2 δβx/L for |x| ≤ L/2 should correspond to the time-evolution of a classical particle of mass m * = + /(2Cd 2 ) in the presence of a constant force F 0 = 2δβ/L along the −x direction. Borrowing the terminology from condensed matter physics, we call such a wave packet with positive effective mass electron-type or "e-type". Equivalently, an initial state |ψ h with components near the top of the band, k ∼ ±π, corresponds to a classical particle with mass m * = − /(2Cd 2 ) and will be called hole-type or "h-type". We remind the reader that choosing purely real components α m for the initial wave packet ensures that the initial velocity of the classical particle is zero. Based upon this analysis, it follows that the average position of the wave packet x(t) will satisfy where the negative sign is for an "e-type" wave packet, the positive sign is for an "h-type" wave packet, and x(0) is the initial location of the wave packet. Figure 4 shows the numerically obtained results for time evolution of a wave packet in an array with N = 201 waveguides and δβ = 5C. The top left panel shows the site and timedependent intensity I(p, t) of an "e-type" wave packet with initial Gaussian profile of size M = N/10 1 at the center of the array. We see that, in a sharp contrast with earlier results, the wave packet largely maintains its shape and moves toward the region with lower potential or, equivalently, smaller waveguide index, in a parabolic manner. The top right panel shows corresponding results for an identical "h-type" wave packet; it, too, maintains the shape, but moves towards larger waveguide index in a parabolic manner. We emphasize that in both cases, the external linear potential is identical; the opposite motions of the "e-type" and "h-type" wave packets arise due to their equal but opposite effective masses, and subsequent accelerations. These observations are quantified in the bottom panel where we plot the mean position of the wave packet, j mean (t) as a function of normalized time for the "e-type" (dashed red) and "h-type" (solid blue) wave packets. It is clear that they follow Eq. (13) where the magnitude of dimensionless acceleration is given by , and matches the acceleration obtained from a quadratic fit to the data shown in the bottom panel. We emphasize that as the wave packet gets closer to the edge, the contribution from reflected partial waves increases and destroys its mapping onto a classical non-relativistic particle.
These results show that a waveguide array with constant nearest-neighbor tunneling can be used to investigate properties of a quantum particle or a non-relativistic classical particle in an external potential. It also has the special property that the bandwidth of its corresponding Hamiltonian, ∆ B = 4 C, does not depend upon the number N 1 of waveguides in that array; this Nindependence ensures the existence of the thermodynamic limit for such a lattice. However, as we discussed in the introduction, waveguide arrays offer the possibility of a sitedependent, nearest-neighbor tunneling C k,k+1 = C k+1,k = C(k). In the following subsection, we present the properties of arrays with such position-dependent tunneling profiles. The left-hand panel shows the intensity I e (p, t) for the "e-type" wave packet, which simulates a particle with positive effective mass, whereas the right-hand panel shows the corresponding result I h (p, t) for an "h-type" wave packet, which simulates a particle with negative effective mass. In contrast with the earlier results, here the wave packets (mostly) maintain their shape as they move towards lower or higher potential, respectively, in a parabolic manner. The bottom panel shows that the mean positions of the two wave packets, j mean (t), obtained from the time-dependent intensity distributions, follow the trajectory of a non-relativistic particle with constant acceleration and zero initial velocity. quantum particle is shifted towards the end of the array with site index N . On the other hand, if C(k) is a rapidly oscillating function of site index, C(2k) C(2k +1), then the N -site array is best understood in terms of N/2 weakly coupled dimers with tunneling profile C(2k+1) where each dimer represents two adjacent waveguides with a strong tunneling C(2k) between them [45]. In general, the tunneling profiles in both of these models break the paritysymmetry about the center of the array, C(k) = C(N −k), and thus prefer one end of the array over the other.

Arrays with site-dependent tunneling profiles
To maintain the equivalence between two ends of a finite, N -site array, we restrict ourselves to Hermitian tunneling profiles that obey C(k) = C(N −k). In the simplest case, this constraint implies that the tunneling profile has either a single maximum or a single minimum at the center of the array. Therefore, we consider single-parameter tunneling functions When α > 0, the tunneling rate at the center of the array is (N/4) α/2 times larger than the tunneling near its edges, whereas when α < 0, the converse is true; when α = 0, we recover the constant-tunneling case. Since the tunneling amplitude C(k) can be varied by a factor of hundred in a single material [5,6,46,47], establishing such tunneling profile constraints the size of the array to (N/4) |α|/2 ∼ 100 or, equivalently, N ≤ 10 4 for |α| = 1, N ≤ 200 for |α| = 2, and N ∼ 20 for |α| = 3. These numbers show that it is feasible to fabricate waveguide arrays with a reasonable number of waveguides for tunneling profiles up to |α| ≤ 3.
The Hamiltonian for such an N -site array is given by We remind the reader that when α = 0, due to the loss of translational invariance, the eigenstates of the Hamiltonian are not labeled by momentum and, in general, it is not possible to obtain analytical solutions for the eigenvalues and eigenfunctions. The sole, notable exception is the case with α = 1, where analytical solutions for the eigenvalues and eigenfunctions are possible [48,49,50]. One can, however, show that energy eigenvalues of H α for any α occur in pairs ±E n and that the corresponding eigenfunctions are related by a simple transformation [51]. Figure 5 shows the typical properties of Hamiltonian H α , for an array with N = 500 and |α| ≤ 2 obtained numerically. The left-hand four-panel figure shows the energy eigenvalues normalized by their respective maximum for α = {0, 1, 2, −1} (clockwise). For α = 0, we get the wellknown cosine-band. When α = 1, we obtain a spectrum with equidistant energy eigenvalues, maximum eigenenergy E max = (N − 1) C, and level spacing ∆E = 2 C; for α = 2, the spectrum is linear near the band edges, with a flatter region in between. For α = −1, the spectrum consists of a few localized states near the band edges (shown by the blue oval) along with a bulk of extended states [50]. The four panels on the right-hand side show the unnormalized density of eigenstates D( ), which provides a measure of number of eigenstates available in a small interval δ around energy for α = {0, 2, −2, −1} (clockwise). For α = 0, we recover the well-known result for a one-dimensional lattice with van-Hove singularities, signaled by a diverging D( ), at the band edges [43,52]. For α = 1, due to the equidistant energy levels, the density of states is a constant. When α = 2, the density of states has a maximum near zero energy, consistent with the small slope of the corresponding energy spectrum near = 0. When α = −1, the D( ) has two distinct features. The first is a two-peaked structure that represents the density of bulk, extended states; the second is the presence of discrete, localized states near the band edges (shown by the blue oval). When α = −2, these features are preserved, but there are a number of localized states at different energies; note that the logarithmic vertical scale in this panel shows the distributed weight of such states. These results show that arrays with α-dependent tunneling have widely tunable spectra.
We define the energy bandwidth as ∆ α (N ) = E max − E min = 2E max . When α = 0, the bandwidth is independent of the array size for N 1, ∆ α=0 (N ) → ∆ B = 4 C, whereas for α = 0, the bandwidth depends upon the size  (15) with N = 500 waveguides. The α = 1 spectrum is exactly linear, whereas for α = 2, it is linear near the edges. When α = −1, the tunneling at the edge of the array is higher than that at its center, and the spectrum has discrete, localized states with energies near the band edges (the blue oval). On the right, when α = 0, D( ) is maximum at = ±1 whereas for α = 2, it is maximum at = 0. The quasilinear behavior of the α = 2 spectrum near the band edges is reflected in the flat D( ) near = ±1. For α < 0, the presence of discrete, localized states at the bottom and the top of the energy band is reflected in the finite, but vanishingly small, density of states away from the center of the band.
of the array and is essentially determined by the maximum tunneling element in the array. Thus, ∆ α (N ) ∼ N α for α > 0 and ∼ N −|α|/2 for α < 0. In the following, we use inverse-bandwidth as the characteristic unit of time for an array with a given tunneling profile α and number of waveguides N , τ α (N ) = /E max = 2 /∆ α (N ). Thus, as α > 0 increases, the characteristic time τ α and the characteristic length l α = cτ α /n both decrease, where c/n is the (constant) speed of light along the waveguide with index of refraction n. Thus, in a sample with a given physical length, long-time dynamics are easily observed as α increases, whereas short-time dynamics become accessible for α < 0 [53]. Now we consider the time evolution of a wave packet in such an array. For an arbitrary initial state |ψ(0) , the time-evolved state is obtained by |ψ(t) = G α (t)|ψ(0) where the time-evolution operator G α (t) = exp [−iH α t/ ] is obtained numerically. Since we have discussed the timedependent intensity profiles of wave packets that are localized to a single or two sites in Sec. 2.1, here we choose a broad initial state that is equally distributed across all waveguides, |ψ(0) = 1/ √ N . Fig. 6 shows the intensity I(p, t) = | p|G α (t)|ψ(0) | 2 in an N = 50 array with α = {0, 1, 2, −1}; the horizontal axis denotes time normalized by the α-and N -dependent time-scale τ α (N ). Note that, due to the symmetries of the Hamiltonian and the initial state, the intensity satisfies I(p, t) = I(N + 1 − p, t) and that the average intensity per site is I a = 0.02 = 1/N . When the tunneling is constant, the effects of interference and reflection at the boundaries lead to a suppression of the intensity at the edges, and a modest enhancement, by a factor of five, near the cen- For α = 0, we get larger intensity in the central region due to edge-reflection and interference. α = 1 shows periodic behavior due to the equally-spaced eigenvalues of the underlying Hamiltonian. When α = 2, the quasilinear energy spectrum and the edge-reflections contribute to the quasi-periodic larger intensity in the central region. The bottom panel shows that for α = −1, eigenstates localized at the two edges lead to a larger intensity at the edge instead of in the central region.
ter of the array (α = 0 top panel). When α = 1 (second panel) the constant spacing between the energy levels implies that the intensity profile is periodic in time, I(p, t) = I(p, t + π/C). In contrast to the constant tunneling case, we also observe that the maximum intensity at the center of the array is enhanced by a factor of 20. For α = 2 (third panel) due to the quasilinear nature of the energy spectrum, we see approximate reconstruction of the intensity profile, and the maximum intensity at the center is again significantly enhanced from its initial value. In all the three cases, since the tunneling at the center is maximum, we see that the intensity profile I(p, t), in general, is largest at the center of the array and reduces symmetrically on the two sides. The bottom panel in Fig. 6 shows the intensity evolution for an array with α = −1, which has localized eigenstates at the two ends of the array. In a sharp contrast with the earlier results, we see that I(p, t) now shows symmetrical maxima near the two edges of the array, with a broad minimum near the central region. These results show that identical initial states give rise to strikingly different intensity profiles in tunable waveguide arrays with a position-dependent tunneling profiles.

Disorder induced localization
In the past three subsections, we have focused on the properties of waveguide arrays with constant or position dependent tunneling profiles and constant or linearly varying on-site potentials; we implicitly assumed that it was possible to fabricate a waveguide array with the exactly specified Hamiltonian. This is, of course, an approximation. In real samples, disorder is always present through variations in the tunneling amplitudes C j,j+1 or on-site potentials β j in the tight-binding Hamiltonian, Eq.(1). The effect of such disorder on the transport properties of lattices was first investigated in the context of electronic systems [54,55], and then extended to classical waves [15,16,17]. In one dimension, all eigenstates of a disordered Hamiltonian are exponentially localized in the limit of an infinite system size, N → ∞ irrespective of the strength of the disorder v d . This non-analytical result -exponential localization at infinitesimal disorder -is due to the subtleties associated with the order of limits N → ∞ and v d → 0 [56,57,58].
In a finite array of N ∼ 10 2 coupled waveguides, localization refers not to an exponential localization of all eigenstatesà la electronic systems, but rather to the development of a "steady-state" intensity profile I(p) that contrasts the ballistic expansion and edge-reflection present in a clean system. The time required for the emergence of the steady-state profile is inversely proportional to the strength of the disorder. In another sharp contrast, the typical strength of disorder in (weakly conducting) electronic materials is v d E F whereas in waveguide arrays, the disorder strength can be comparable to the tunneling, v d ∼ C [13,14].
In this subsection, we present the effects of disorder on the time-evolution of a uniform initial state. We consider two distinct disorders. The diagonal disorder ran-domly modulates the on-site potential β i → β i + v i where v i is a random variable with zero mean and variance v ds . The off-diagonal disorder randomly modulates the tunneling C i → C i + v i where v i is a zero-mean random variable with variance v dt . We use uniformly distributed random variables to ensure that the modulated tunneling rates remain strictly positive, although the results are independent of the type of distribution used as long as any such distribution has zero mean and identical variance [14,59]. The resultant intensity distribution is averaged over multiple M ∼ 10 4 realizations to ensure that the final results are independent of the number of disorder realizations and the probability distribution of the site or tunneling disorder. Figure 7 shows the intensity profile I(p, t) for an array with N = 50, uniform initial state, and α = {0, 2, −1} where · · · denotes disorder average. We remind the reader that the average intensity per site is I a = 1/50. The left-hand column has results for on-site disorder v ds and the right-hand column has results for the tunneling disorder of equal strength, v ds = v dt = 0.1∆ α (N ). The top line, α = 0, shows that for both disorders the initial interference pattern is replaced at later times by a steady state intensity that is suppressed at the edges. The center line, α = 2, shows the same qualitative behavior, but also shows slight difference between the the two intensity profiles, particularly at small times. The bottom line, α = −1, shows steadystate profiles that have maxima near the two edges. In all cases, the differences between the left-hand and righthand panels for a given tunneling profile C α (j) decrease with increasing time, measured in units of τ α (N ).   In all cases, the interference pattern at small times is replaced by quasi steady-state intensity at large times. For α = 0 (top line) and α = 2 (center line), I(p) has a maximum near the central region, whereas for α = −1 (bottom line) the intensity has multiple maxima near the two edges of the array. This emergence of steady state profiles shows that "extended" initial states also undergo disorder-induced "localization" as it is defined here.
Lastly, we compare the cross-section of the intensity profiles at t/τ α = 500 for the same array with on-site disorder (solid symbols) and tunneling disorder (open symbols) of equal strength, v ds = v dt = 0.1∆ α . When α = −1 (circles), the intensity profile shows a minimum at the center and multiple, symmetric maxima at the two edges, whereas for α = 2 (squares), the intensity is maximum at the center and monotonically decays away from it. Note that the multiple maxima near the two edges show up as striations in the intensity profiles for α = −1 in Fig. 7. The (gray) dashed line shows the average intensity I a = 1/N = 0.02 per site. We point out that the intensity profiles for on-site and tunneling disorders coincide with each other at sufficiently long times, although the time required for such a match depends upon the tunneling profile α and the initial state. For example, Fig. 8 shows virtually identical intensity profiles for α = 2, whereas for α = −1, the intensity suppression due to the tunneling disorder (black open circles) is larger than that by the on-site disorder (blue solid circles). It is also worth emphasizing that the disorder-averaged intensity profile recovers the underlying parity-symmetry shared by the clean Hamiltonian and the initial state, I(p, t) = I(N + 1 − p, t) .
We end this section with another phenomenon due to the parity-symmetric tunneling profile in a finite array of waveguides. Fig. 9 shows the time-and site-dependent intensity evolution in an array with N = 100 waveguides, a small on-site disorder v ds /∆ α = 0.05, and tunneling profiles with α ≥ 0. The initial wave packet is localized at a single site m 0 = 15. The top panel (α = 0) shows that I(p, t) changes from interference-dominated behavior at short times to disorder-dominated steady-state behavior at longer times; the steady-state intensity is maximum at site m 0 and decays exponentially with distance from m 0 [13]. The center panel (α = 1) shows that, at short times, the wave packet partially reconstructs at the parity-symmetric sitem 0 = N + 1 − m 0 . The steady-state intensity profile in this case has two peaks, at m 0 and m 0 , and their relative weights are tuned by the disorder strength and the distance between the two peaks. The bottom panel (α = 2) shows a qualitatively similar result. Thus, a position-dependent, parity-symmetric tunneling in a finite array of waveguides leads to effective localization at two waveguide locations, even if the initial wave packet is introduced in a single waveguide [60].

Non-Hermitian, PT -symmetric Models
In the last section, we only considered Hermitian Hamiltonians, Eq.(1), which modeled waveguides that have no loss or amplification of the input signal. The ubiquitous losses that are present in real waveguides are phenomenologically taken into account by adding a negative imaginary part to the real eigenvalues of the Hermitian Hamiltonian, E n → E n − iΓ n [61,62]. This imaginary part Γ n > 0 leads to an exponential decay of the total intensity and therefore represents dissipation, absorption, or friction [3]. Nominally, if we assign a positive imaginary part to the energies of a Hermitian Hamiltonian, E n → E n + iΓ n with Γ n > 0, the total intensity of an initially normalized wave packet  Fig. 9. Intensity I(p, t) in an array with N = 100, a weak disorder v ds /∆ α = 0.05, and initial wave packet |m 0 with m 0 = 15. Top panel shows that for constant tunneling, the steady-state intensity profile is maximum at m 0 , with exponential decay on the two sides. The center (α = 1) and bottom (α = 2) panels show that, at short times, the wave packet partially reconstructs at the mirror-symmetric sitem 0 = N + 1 − m 0 = 86. Thus, in sharp contrast with the traditional localization, α = 0 case, the steady-state intensity profiles for α ≥ 1 have a two peaks, one at the initial wave packet location and the other at its parity-symmetric counterpart. Note that in all cases, the average intensity is I a = 0.01 = 1/N and thus, the localization enhancement is only by a factor of two.
will increase, and will therefore represent gain or amplification. Such a phenomenological model breaks down at long times, when the power required to maintain the exponential intensity increase cannot be supplied by the "reservoir". In this section, we will focus on non-Hermitian Hamiltonians that represent balanced, spatially separated loss and gain. In a waveguide-array realization of such a Hamiltonian, one of the waveguides is lossy, its parity-symmetric counterpart has gain, and the rest of the waveguides are neutral [37,38]. To get a feel for properties of such a system and to define the terminology, let us start with the simplest example with N = 2 waveguides. The tunneling Hamiltonian for this system is given by H t = − C(a † 1 a 2 + a † 2 a 1 ). The non-Hermitian, PT -symmetric potential, which represents gain in the first waveguide and loss in the second, is given by V = i γ(a † 1 a 1 − a † 2 a 2 ). In a matrix notation, the total Hamiltonian becomes Although H = H t + V is not Hermitian, it is invariant under the combined parity (P : 1 ↔ 2) and time-reversal (T : i → −i) operations [2]. It is straightforward to obtain the eigenvalues λ ± and (right) eigenvectors |± R of the Hamiltonian (16). We remind the reader that since the matrix H is not Hermitian, its left-eigenvectors and right-eigenvectors are not Hermitian conjugates of each other [63,64]. For a small non-Hermiticity, γ ≤ C, the eigenvalues of H are purely real, and given by λ ± = ± = ± C 2 − γ 2 . The corresponding right-eigenvectors are given by where sin θ = γ/C ≤ 1. Thus, R +|− R = (−i)e iθ sin θ = 0. Since the matrix H is symmetric, H = H T , the lefteigenvectors are obtained by taking the transpose of the right-eigenvectors. |± R are simultaneous eigenvectors of the combined PT operation as well, and each of them has equal weight on the gain and the loss site. When γ = 0 the inner product is zero, whereas for γ → C, the two eigenvalues become degenerate and the two eigenvectors become parallel to each other. For γ ≥ C, the eigenvalues are purely imaginary complex conjugates, λ ± = ±i Γ = ±i γ 2 − C 2 . The corresponding right-eigenvectors are now given by where cosh φ = γ/C ≥ 1. Thus, the inner product of the two eigenvectors is equal to 1/ cosh φ ≤ 1. Note that now the eigenvectors are not simultaneous eigenvectors of the PT -operation; the |− R eigenvector has higher weight on the gain site and the |+ R eigenvector has higher weight on the loss site. The region of parameter space where all eigenvalues are real and the eigenvectors are simultaneous eigenvectors of the PT operation, γ/C ≤ 1, is traditionally called the PT -symmetric region, and γ P T = C is called the threshold loss-and-gain strength. For γ/γ P T > 1, complex conjugate eigenvalues emerge and the PT -symmetry of the Hamiltonian H is not shared by its eigenvectors with complex eigenvalues. Therefore, the emergence of complex eigenvalues is called PT -symmetry breaking. In the following subsections, we present the properties of Nwaveguide arrays with Hermitian, position-dependent tunneling profiles C α (j) and a single pair of non-Hermitian, PT -symmetric loss and gain potentials.

PT symmetric phase diagram
We begin with the Hamiltonian for an N -site array with open boundary conditions, where H α is the Hermitian tunneling Hamiltonian, Eq. shows that γ P T (µ)/∆ α = 1 is maximum at µ = 0.5, and remains relatively constant over a wide range of µ as µ → 0 for α ≥ 1. In contrast, for constant tunneling, the threshold loss-and-gain strength decays rapidly with decreasing µ, but increases again as µ → 0. The right-hand panel shows corresponding, qualitatively similar, results for an odd array with N = 101. For an odd array, the smallest separation between loss and gain waveguides is D = 2, instead of D = 1 in an even array, and therefore, the maximum threshold value for an odd array near µ = 0.5 is γ P T = 0.5∆ α .
P : a k → ak. Thus, it follows that the Hermitian part of the Hamiltonian is PT -symmetric, C α (k) = C α (N − k), and so is the non-Hermitian potential term. Thus, to obtain the PT -symmetric phase diagram, we need to obtain the eigenvalues of the Hamiltonian H P T and then locate the threshold loss and gain strength γ P T (µ) as a function of the relative location µ = m/N of the gain waveguide. It is possible to obtain this threshold analytically only in the case of constant tunneling, α = 0 [65,66]; however, for an arbitrary α, a numerical approach is most fruitful. By numerically tracking the emergence of complex eigenvalues of the tridiagonal matrix H P T α , we obtain the typical phase diagram, shown in Fig. 10. Note that µ = 1/N corresponds to largest distance between the loss and gain waveguides, whereas µ ∼ 0.5 corresponds to the shortest separation between them. Due to the constraint of paritysymmetric locations, in an even N -array this separation is unity, and for an odd N -array, the loss and gain locations have to be separated by a single waveguide between them.
The left-hand panel in Fig. 10 shows the threshold strength measured in units of the lattice bandwidth as a function of relative location of the gain waveguide for an N = 100 array with α = 0 (blue circles), α = 1 (red squares), and α = 2 (beige diamonds); all eigenvalues of H P T α are real for values of γ below the curve for that α. Note that we use quarter-bandwidth, ∆ α = ∆ α /(4 ), as the relevant scale in the phase diagram. For α ≥ 1, the threshold strength is maximum γ P T /∆ α = 1 is at µ = 0.5, when the loss and gain waveguides are nearest neighbors. It reduces to γ P T /∆ α ∼ 0.3 and remains approximately constant for 0.15 ≤ µ ≤ 0.45, and is monotonically suppressed with the separation D = 1+N (1−2µ) between the loss and gain waveguides. Note that the behavior γ P T (µ) for an array with constant tunneling amplitude, α = 0 is dramatically different. Starting from the maximum value of γ P T /C = 1 for closest loss and gain, the threshold strength first drops rapidly with increasing d, but is again enhanced as the loss and gain sites approach the two edges of the array. Thus, for moderate separations µ ∼ 0.25 and number of waveguides N ∼ 100, the PT -symmetric phase in an array with non-constant tunneling amplitudes is substantially stronger than in an array with constant tunneling amplitude. The right-hand panel shows the PT -phase diagram for an array with an odd number of waveguides, N = 101. We see that the robust nature of the PT -symmetric phase for α ≥ 1 is maintained, although the threshold for smallest separation µ = (N − 1)/2N is reduced to γ P T = 0.5∆ α [66,67].
We emphasize that although the qualitative form of the PT -phase diagram is the same for different N , as N increases, the threshold strength γ P T (µ)/∆ α decreases for all separations except when the loss and gain are the closest (µ ∼ 0.5) or the farthest (µ = 1/N ). Thus, rigorously, γ P T /∆ α (N ) → 0 as N → ∞; however, this is of no concern for experiments where the number of sites in an array -whether the "site" be an optical waveguide [37,38,40], an RLC circuit [68], or a pendulum [69] -is typically N 100. Since the PT -symmetry breaking occurs when two adjacent eigenvalues E n , E n+1 become degenerate and then complex, and since the eigenvalues of H α occur in pairs ±E n , it follows that, for a generic position µ of the gain waveguide, N − 4 eigenvalues of the Hamiltonian H P T α remain real while four eigenvalues become complex conjugate pairs. The remarkable exception to this rule is the case of nearest-neighbor loss and gain waveguides in an even N array. In this case, since the array can be effectively divided into two systems, one with the loss and the other with the gain, all N eigenvalues of H P T α become complex simultaneously [67,70]. Thus, the implications of PT -symmetry breaking are determined by both the threshold loss-and-gain strength γ P T and the location and number of eigenvalues that become complex at the threshold.

Time evolution across the PT threshold
In the previous section with Hermitian Hamiltonians, we presented intensity profiles I(p, t) for various initially normalized states, ψ(0)|ψ(0) = 1. Since the time evolution operator in these cases is unitary, G † (t) = exp +iH † t/ = exp [+iHt/ ] = G −1 (t), the total intensity of the timeevolved wave packet remains unity, N p=1 I(p, t) = 1. For a non-Hermitian Hamiltonian, since H † P T = H P T , the corresponding time evolution operator is not unitary. Therefore, the norm of an initially normalized state is not preserved and the total intensity is a function of time, I(t) = N p=1 I(p, t) = 1. Note that G(t) = exp [−iH P T t/ ] is not a unitary operator irrespective of whether the system is in the PT -symmetric phase or has complex conjugate eigenvalues.
To get a better feel for this non-unitary time evolution operator, let us calculate it for the two-site Hamiltonian, Eq. (16). From the completeness property of its left and right eigenvectors, it follows that where the left eigenvectors L ±| are obtained by transposing the right eigenvectors |± R . In the PT -symmetric phase, γ/C ≤ 1, Eq.(17) implies that G ≤ (t) = cos τ + γ sin τ +i C sin τ +i C sin τ cos τ − γ sin τ = G T ≤ (t) (21) where τ = t/ is the dimensionless time. We leave it to the reader to verify that G ≤ (t) is not unitary, but its eigenvalues have unit modulus and are given by e ±iτ . Therefore the non-unitary time evolution operator satisfies det G ≤ (t) = 1. In the PT -symmetry broken phase, γ/C ≥ 1, a corresponding calculation using Eq. (18) gives where τ = Γ t/ . The reader can verify that G ≥ (t) is not unitary, its eigenvalues are e ±τ , and thus, det G ≥ (t) = 1. We note that the matrix elements of G ≤ (t) are bounded, those of G ≥ (t) diverge with increasing time, and that the time evolution operator is continuous across the PTsymmetry threshold. At the threshold γ = C, since the Hamiltonian is singular, H 2 = 0, the exponential expansion for the time-evolution operator truncates at the linear order and gives Since the time evolved state is given by |ψ(t) = G(t)|ψ(0) , the change in net intensity is proportional to unitary deficit, G † (t)G(t) − 1. Equations (21)- (23) show that, for a PTsymmetric Hamiltonian (16), the net intensity I(t) in the PT -symmetric phase remains bounded, increases exponentially with time in the PT -symmetry broken phase, and exactly at the threshold, varies quadratically with time at long times [71]. Figure 11 shows the evolution of net intensity I(t) in an N = 40 waveguide array with constant tunneling, α = 0, the loss-and-gain waveguides farthest apart (m = 1) or closest together (m = N/2 = 20) as a function of γ/γ P T . These numerically obtained results are for an initial state localized at the first waveguide, |ψ(0) = |1 . We remind the reader that the crucial difference between the m = 1 case and the m = N/2 case is that only four eigenvalues, at the center of the cosine-band become complex for m = 1, whereas all eigenvalues simultaneously become complex when m = N/2 [66,67,70]. The top-left and bottom-left panels show that in the PT -symmetric phase, γ/γ P T < 1, the net intensity I(t) oscillates but remains bounded, and its time-average increases monotonically with its proximity to the PT -symmetric phase boundary. In addition, they show that the average and fluctuations in the m = N/2 case are smaller than those in the m = 1 case. The top-right panel shows I(t) at the threshold, γ/γ P T = 1, for the two cases; note the logarithmic scale on both axes. At small times, the orderof-magnitude difference between intensities for the two gain-waveguide locations is consistent results in the lefthand panels. At longer times, we see that the net intensity scales quadratically with time, although the prefactor of this quadratic dependence is greater for the m = 1 case.
The bottom-right panel shows I(t) in the PT -symmetry broken phase, γ/γ P T = 1.01; note the logarithmic scale on the vertical axis. These results show that, as expected, the net intensity diverges exponentially, but with a larger exponent for loss and gain waveguides at the two ends of the array, m = 1. We emphasize that this qualitative trend is valid for arbitrary location and shape of the initial wave packet. The results in Fig. 11 show that the simple 2 × 2 non-Hermitian Hamiltonian, Eq.(16), captures the timedependence of the net intensity in a large tight-binding array, although it does not capture the full gamut of PTsymmetry breaking signatures [67].

Intensity correlations with Hermitian or PT -symmetric disorders
We have seen in Sec. 2.4 that (Hermitian) disorder leads to "localization" of an arbitrary initial state, that is characterized by a steady-state, disorder-averaged intensity profile I(p) . The steady-state intensity profile is solely determined by the initial state and the strength of the disorder potential, but is independent of whether the disorder is in the on-site potentials or tunneling amplitudes. Therefore, the site-dependent steady-state intensity measurements can only determine the strength of the disorder, but not the type of the disorder. These two disorders affect the particle-hole symmetric spectrum of the clean show that I(t) is remains bounded in the PT -symmetric region, γ/γ P T < 1. The top-right panel shows that I(t) ∝ t 2 at the threshold, γ/γ P T = 1. The bottom-right panel shows that in the PT -symmetry broken region, γ/γ P T > 1, the net intensity diverges exponentially with time. These results, obtained for |ψ(0) = |1 , have the same qualitative behavior for an arbitrary initial state.
lattice in qualitatively different manners: the on-site, diagonal disorder destroys this symmetry whereas the tunneling, off-diagonal disorder preserves it. Therefore, although intensity measurements are insensitive to it, it is known that intensity correlation function is able to distinguish between the on-site and tunneling disorders [72]. In contrast to the Hermitian potential, a non-Hermitian, PT -symmetric potential, in the PT -symmetric phase, preserves particle-hole symmetry of the resulting, purely real spectrum [51]. Therefore, in this section, we compare the steady-state intensity correlations from a Hermitian tunneling disorder and a non-Hermitian, PT -symmetric disorder, both with zero mean and equal variance. The PTsymmetric disorder potential is given by where the random, loss (or gain) potentials |γ m | ≤ γ P T (µ = m/N ) ensure that the system is in the PT -symmetric phase. The normalized correlation matrix is defined as [72] Γ jk (t) = I(j, t)I(k, t) I(j, t) I(k, t) t 1 (25) where I(j, t) is the intensity profile determined by the initial state |ψ(0) and the disorder potential. I(j, t) is the disorder-averaged intensity that becomes independent of time at long times (Sec. 2.4). The net intensity p I(p, t) is conserved at unity for a Hermitian disorder, but not for the PT -symmetric disorder. The intensity correlation function is defined as and represents the sum of weights along a diagonal that is shifted by r from the main diagonal of the steadystate correlation matrix, Eq. (25). Figure 12 shows the normalized, steady-state correlation matrix Γ ij and the intensity correlation function g(r) for an N = 20 array with constant tunneling, C α (j) = C, and initial state |ψ(0) = (|9 + |10 )/ √ 2. The top line shows the results for a PT -symmetric, on-site disorder, whereas the bottom line has results for a Hermitian, tunneling disorder; both disorders have zero mean, equal variance v d /∆ B = 0.02, and the results are averaged over M ∼ 10 4 disorder realizations. Panels (a) and (c) show that the full correlation matrix Γ jk is sensitive to the source of disorder. However, panels (b) and (d) show that the intensity correlation function g(r) = g(−r) cannot distinguish between the two. Thus, symmetry properties of the disorder-induced spectrum are reflected in the disorder-averaged intensity correlation function, and not the on-site or off-diagonal nature of disorder [60].
These results also suggest that although intensity distribution, or intensity correlation function is insensitive to the disorder distribution function, higher order intensity correlations may encode signatures of different disorder distributions that have zero mean and identical variance [14,59].  Fig. 12. Normalized correlation matrix Γ jk and intensity correlation function g(r) for an N = 20 array with constant tunneling, PT -symmetric on-site disorder (top line) and Hermitian, tunneling disorder (bottom line) with zero mean and equal variance v d = 0.02∆ B . The steady-state Γ jk , panels (a) and (c), are different for the two sources of disorder, whereas the steady-state intensity correlation function g(r), panels (b) and (d), is insensitive to them. Their similarity shows that the particle-hole symmetry of the disordered spectrum is instrumental to the correlation function properties.

Conclusion
In this article, we have presented the properties of coupled waveguide arrays. We have argued that they provide a versatile and robust realization of a tight-binding model, ideally suited for investigating many quantum, quasi-classical, and bandwidth effects that are not easily accessible in "naturally occurring" lattices in electronic materials. We have shown that finite arrays with small number of waveguides exhibit a rich variety of effects, such as localization in the parity-symmetric waveguide, that are absent in a lattice with sites N → ∞.
Due to the ease of introducing absorption or amplification, coupled optical waveguides are also well-suited to model open systems with spatially separated, balanced loss and gain. Such systems are formally described by non-Hermitian, PT -symmetric Hamiltonian. Since the spectrum of such Hamiltonian changes from purely real to complex, and since the time-evolution under such Hamiltonian is always non-unitary, we have discussed a few salient properties of PT -symmetric lattice models.
In this review, we have ignored nonlinear effects that arise at high intensities in a waveguide, and that are expected to play a large role in the PT -symmetry broken region where the net intensity increases exponentially with time. We have not considered the effects of shape-preserving solitonic solution that exist in the nonlinear regime on time evolutions discussed here. In addition, we have not discussed the effects of PT -symmetric, non-Hermitian disorder, including the fate of Anderson localization, in the PT -symmetry broken region. The investigation of these outstanding questions will further deepen our knowledge of this exciting research area.