Insight of the Green's function as a defect state in a boundary value problem
IInsight of the Green’s function as a defect state in a boundary value problem
Jose D. H. Rivero and Li Ge
College of Staten Island, CUNY, Staten Island, New York 10314, USA andThe Graduate Center, CUNY, New York, NY 10016, USA (Dated: February 5, 2021)A new perspective of the Green’s function in a boundary value problem as the only eigenstate inan auxiliary formulation is introduced. In this treatment, the Green’s function can be perceived as adefect state in the presence of a δ -function potential, the height of which depends on the Green’sfunction itself. This approach is illustrated in one-dimensional and two-dimensional Helmholtzequation problems, with an emphasis on systems that are open and have a non-Hermitian potential.We then draw an analogy between the Green’s function obtained this way and a chiral edge statecircumventing a defect in a topological lattice, which shines light on the local minimum of the Green’sfunction at the source position. I. INTRODUCTION
Initially distributed to just 51 private members of asubscription library in 1828 [1], the essay of George Greenon the application of mathematical analysis to the theoriesof electricity and magnetism has inspired generations ofphysicists and mathematicians. By studying the three-dimensional Laplace equation governing the static electricpotential, Green recognized the unique utility of a functioninversely proportional to the distance between two chargedbodies, which now bears his name and preceded Dirac’sintroduction of the δ -function by more than a century.The importance of the Green’s function, together withthe Green’s theorem, was recognized by masterminds ofmodern physics such as Lord Kelvin and Julian Schwinger.Besides its effectiveness as a theoretical and numericaltool to solve linear differential equations [2, 3], the Green’sfunction provides a fundamental connection between thesource and the field [4]. When extended to address dy-namical equations and their spectral representations, theGreen’s function becomes the propagator that is funda-mental to field theories [5]. As such, it is a powerful andessential tool to study a variety of transport and scat-tering problems, ranging from condensed matter physics[6, 7], optics and photonics [8–10] to high energy physics[11].In optical physics, one area of interest in this article, theGreen’s function has been routinely employed to analyzerich optical phenomena including light-matter interactionin lasers and other light sources [12–14], as well as todevise applications for optical computing, informationprocessing, and telecommunications. What we offer hereis a novel perspective of the Green’s function, i.e., treatingit as the single eigenstate in an auxiliary boundary valueproblem. In addition to further enrichment and shapingof our physical intuition through the Green’s function, wefind exceptional parallels between the Green’s functionand defect states due to a local potential, including achiral edge state circumventing a defect on its path in atopological lattice.Below we introduce formally the Green’s function prob-lem and lay out the fundamentals of our approach. TheGreen’s function of an operator L in a variety of physics problems can be defined by[ z − L ] G ( r , r (cid:48) ; z ) = δ ( r − r (cid:48) ) (1)with a proper boundary condition. Here z is a parameter(e.g., the energy) and r , r (cid:48) are the coordinates of thefield and the source, respectively. A new perspective ofthe Green’s function can be obtained using the followingauxiliary eigenvalue problem[ z − L ] ψ m ( r , r (cid:48) ) = λ m ( r (cid:48) ) f ( r , r (cid:48) ) ψ m ( r , r (cid:48) ) . (2)Below we suppress the r (cid:48) -dependence of ψ m and λ m forconciseness. As we will show, this auxiliary problem hasa single eigenvalue λ and eigenstate ψ ( r ), when f ( r , r (cid:48) )is chosen to be δ ( r − r (cid:48) ). The Green’s function is thenuniquely determined by G ( r , r (cid:48) ) = ψ ( r ) λ ψ ( r (cid:48) ) , (3)and λ − gives the value of the Green’s function at thesource r = r (cid:48) . As a bonus, we obtain directly the localdensity of states (LDOS) that is proportional to the imagi-nary part of G ( r , r ; z ) [7], i.e., Im[ λ − ]. The reciprocity ofthe Green’s function, though implicit in Equation (3), canbe easily verified in the absence of an effective magneticfield, as shown in Appendix A].Before we apply this approach to various Hermitianand non-Hermitian problems, we first prove Equation (3)and discuss the general properties of ψ m ( r ). With thechoice of f ( r , r (cid:48) ) mentioned above, the right hand side ofEquation (2) for ψ ( r ) is simply λ δ ( r − r (cid:48) ) ψ ( r ) = λ δ ( r − r (cid:48) ) ψ ( r (cid:48) ) , (4)from which Equation (3) follows directly by comparingwith Equation (1), with the requirement ψ ( r (cid:48) ) = 0. Ifthere were another eigenstate ψ m =0 ( r ), then by repeatingthe same procedure we would find that the Green’s func-tion is proportional to ψ m =0 ( r ) as well, which contradictsthe uniqueness of the Green’s function with a properlyimposed boundary condition. When implemented numeri-cally, there do exist spurious eigenvectors ψ m =0 ( r ), which a r X i v : . [ c ond - m a t . m e s - h a ll ] F e b however can be easily discarded due to their ill-behaved λ m ’s as we will discuss in Section II A.The auxiliary eigenvalue approach equips us with aconceptually different way to treat the Green’s function,i.e., as a defect state [15] emerging due to the δ -functionpotential: [ L + λ δ ( r − r (cid:48) )] ψ ( x ) = zψ ( x ) . (5)As we will show, this point of view is particularly interest-ing and helpful in a topological system with chiral edgestates [16–20]. For example, if a point source is placedat the edge of a two-dimensional (2D) topological insula-tor, the auxiliary eigenvalue approach indicates that ananalogy exists between the Green’s function and a chiraledge state circumventing a defect at the same location.This interpretation provides an intuitive understanding ofthe local minimum of the Green’s function at the sourceposition, which we will illustrate using a 2D square latticewith a π/ PT ) symmetry [21, 22], focusingon the Green’s function at an exceptional point (EP). AnEP is a unique degeneracy found only in non-Hermitiansystems, where two or more eigenstates of the systemcoalesce. It has led to a plethora of intriguing phenomena,and in particular, it was shown recently that the Green’sfunction can be fully decoupled from the coalesced eigen-state in a photonic system, which is instead given by theJordan vector or the “missing dimension” of the Hilbertspace [23]. We show that our method based on Equa-tion (2) captures this extraordinary behavior nicely in aring cavity, besides describing correctly the Green’s func-tion in a PT -symmetric photonic molecule. We furthervalidate our method in quasi-1D waveguides, which arefrequently employed to study disordered mesoscopic andoptical systems [24–26]. In Section III, we first discusshow the defect state corresponding to the Green’s functionis conceptually different from previous studies of defectstates in a 1D continuous Hermitian system and a 1Dnon-Hermitian lattice. We then highlight the intriguingmanifestation of the linkage between the Green’s functionand a defect state in the aforementioned 2D topologicallattice. II. VALIDATIONA. 1D Hermitian case
We start with an exmaple where the analytical formof the Green’s function is available. Consider the scalarHelmholtz equation in 1D with a uniform refractive index n ∈ R : L = − n ∂ x , z = k . (6)Here k is the real-valued wave vector in free space. Be-low we take the speed of light in vacuum to be 1 anddo not distinguish k from the (circular) frequency. Weimpose the Dirichlet boundary conditions G ( x, x ) = 0, ψ m ( x ) = 0 at x = 0 , L , which renders the system Her-mitian. Consequently, it can be shown that the Green’sfunction at the source position (i.e., λ − ) is real: λ = kn (cid:20) nk ( x − L )] − nkx ) (cid:21) . (7)When solved using a finite difference scheme [27], thenumerical value of λ given by Equation (3) shows agood agreement with the analytical result given by Equa-tion (7). The corresponding Green’s function, which weshown in Figure 1(a), also agrees nicely with its analyticalexpression G ( x, x ) = ( sin[ nk ( x − L )] λ sin[ nk ( x − L )] , ( x > x ) sin( nkx ) λ sin( nkx ) . ( x ≤ x ) (8)The numerical implementation of the δ -function is usu-ally taken to be the limit of a sharp analytical distribution,such as a Gaussian with the standard deviation σ → x . If we choose σ > ∆ x , nkL/ π D O S / L (a)(b) FIG. 1. The Green’s function (a) and DOS (b) in a 1Ddielectric cavity with perfect mirrors. The green dashed linesare obtained using the analytical expression (8), while theblack solid lines are from our auxiliary eigenvalue approach(3). n = 3, k = 30 /L and x = 0 . L (marked by the arrow)are used in (a), and 2000 grid points are used for the finitedifference implementation of Equation (2). we find spurious eigenstates ψ m =0 ( x ) with almost identi-cal spatial dependence away from the source position (notshown), but they have structures that reside in the finiteextension of the approximated δ -function (e.g., fast oscilla-tions). This problem is remedied by letting σ < ∆ x , whichpractically leads to an approximation of the δ -functionthat has a single non-zero element at the position of thesource, the value of which is given by 1 / ∆ x . This choicewarrants that the integration of the δ -function is 1 overany range enclosing the source and 0 otherwise. Withthis choice, we find that all spurious eigenvalues of Equa-tion (2) approaches infinity (e.g., | λ m =0 L | > in thecase shown in Figure 1), and the sole, physical one λ iseasily obtained by setting the numerical routine to searchfor the eigenvalue with the smallest modulus.As mentioned in the introduction, our auxiliary eigen-value approach also produces the LDOS directly. For aHermitian system, a limiting procedure is needed to regu-late the singularity at real-valued resonant frequencies:LDOS( x ; k ) = lim s → + − kπ Im[ G ( x, x ; k + is )] , (9)i.e., a small positive imaginary part is added to the fre-quency here, constructing in this way the retarded Green’sfunction [See Appendix A]. Equation (9) also applies tonon-Hermitian cases, as long as the complex resonant fre-quencies are on or below the real axis. The integration ofLDOS( x ; k ) over the whole system then gives the densityof states (DOS) as a function of the frequency.To calculate DOS using our approach based on Equa-tion (2), we choose a small s = 0 . /L and calculate − kπ Im[ λ − ] numerically. The result agrees well with theanalytical result [see Figure 1(b)], where λ given byEquation (7) is used. The latter leads toDOS( k ) = X m δ (cid:16) k − mπnL (cid:17) (10)as expected once s →
0, where k m = mπ/nL ( m = 1 , . . . )are the real-valued resonant frequencies. B. 1D non-Hermitian cases
For our next validation, we study the Green’s functionat an EP in a photonic molecule [28]. Such a case presentsa serious challenge to the standard approach based on theeigenvalues of L , i.e., the bilinear expansion G ( r , r (cid:48) ; z ) = X m φ m ( r ) φ m ( r (cid:48) )( z − z m )( m, m ) , (11) L φ m ( r )= z m φ m ( r ) . (12)Depending on the symmetry of L , the partner function φ m ( r ) may or may not be the same as φ ∗ m ( r ), and ( m, n )is the resulting inner product of φ m ( r ) and φ n ( r ). Atan EP, the inner product ( m, m ) becomes zero due to the coalescence of two or more eigenstates of the system[22, 29, 30]. Although this divergence can be eliminatedusing a non-Hermitian perturbation theory [10, 23], itrequires a priori knowledge of the EPs and the generalizedeigenvectors, which adds to the complexity of the problem.Our auxiliary eigenvalue formulation, on the other hand,does not suffer from this drawback.Below we exemplify our method and compare it to theresult of the perturbation theory, without the need toinvoke the Jordan vector that completes the Hilbert spaceat the EP [23]. Our photonic molecule is composed of twohalf-wavelength cavities coupled by a distributed Braggreflector (DBR) placed in air [See Figure 2(a); inset]. If anunbalanced loss distribution is introduced in the two half-wavelength cavities, an emerging effective PT symmetrygoverns the system [31].Here we consider the outgoing boundary condition at areal-valued frequency, which corresponds to the laser andis different from quasi-bound states or resonances with acomplex frequency throughout the entire space [12]. Theresulting complex eigenvalues of L given in Equation (6)are found in the lower half of the complex frequencyplane, as a result of the non-Hermiticity caused by thecavity openness. These complex eigenvalues are knownas continuous flux (CF) states [12, 27], which have alsobeen used to study nuclear decays [32, 33]. We note thatthe EPs of CF states have not been studied before, unliketheir counterparts of quasi-bound states or resonances.In the vicinity of an EP of frequency k , the CF statesin our system can be expressed in terms of waves con-fined in the left- and right-cavity, i.e., ψ ( x ) ≈ a l ψ l ( x ) + a r ψ r ( x ), where the amplitudes a l,r are determined by theHelmholtz equation. Without the PT -symmetric pertur-bation, our photonic molecule is symmetric and hence a l = ± a r in the symmetric and anti-symmetric modeswith CF frequencies ˜ k S , ˜ k A . The introduction of a weak PT -symmetric perturbation in the dielectric function cou-ples the amplitudes a l and a r , determined by their spatialoverlap C with the non-Hermitian perturbation, whichrepresents the strength of gain and loss.The eigenfrequencies of the perturbed system are thenfound to be q ± = k (1 ± ∆ ς ). Here k is the CF frequencyof a single half-wavelength cavity sandwiched by twoDBRs and ∆ is a dimensionless detuning defined by (˜ k S − ˜ k A ) / k . We have also defined ς ≡ p − β , where β = C/ ∆. One EP is reached when ς = 0, resulting in q ± = q . The corresponding eigenstates ψ ± ( x ) are givenby ψ ± ( x ) = ψ l ( x ) + β ± ψ r ( x ) , (13)where β ± ≡ iβ ± ς . The inner product of the eigenstatesis defined as( i, j ) ≡ Z (cid:15) ( x ) ψ i ( x ) ψ j ( x ) dx, ( i, j = ± ) (14)where (cid:15) ( x ) is the dielectric function before the PT -symmetric perturbation is introduced. In other words, x/L | G ( x , x ' ) | / L x/L | G ( x , x ' ) | / L (a)(b) | Δ G | / L (c) -4 -2 Number of CF States rl FIG. 2. Green’s function of a PT photonic molecule at an EP.The source is placed at the center of the right cavity in (a)and in between the two cavities in (b). Inset in (a): Schematicof the photonic molecule and the imposed boundary condition,and labels denoting the refractive indices of each componentof the heterostructure. Green dashed lines correspond to theperturbative expression in Equation (16) and black solid linescorrespond to our auxiliary eigenvalue approach. The locationsof the sources are marked with arrows. (c) Difference betweenthe bilinear expansion with the perturbative correction andthe auxiliary method in (b) at the source position. the partner functions in Equation (11) are chosen as ψ ± ( x ) = (cid:15) ( x ) ψ ± ( x ), similar to the Hermitian case dis-cussed earlier. This definition of the inner product war-rants the biorthogonality (+ , − ) = ( − , +) = 0, and wefind (+ , +) = 2 β + ς →
0, ( − , − ) = − β − ς → , +) , ( − , − ) here at the EP cause a catastrophe inthe calculation of the Green’s function, because the twocorresponding terms diverge independent of the frequency: G ( x, x ; k ) ≈ ψ + ( x ) ψ + ( x )( q − q )(+ , +) + ψ − ( x ) ψ − ( x )( q − q − )( − , − ) . (15)However, it can be shown that the diverging behaviors inthe two terms cancel each other precisely [10, 23], leading to G ( x, x ; k ) ≈ (cid:15) ( x ) ψ l ( x ) ψ l ( x ) + ψ r ( x ) ψ r ( x ) k − k − i ∆ k (cid:15) ( x ) ψ EP ( x ) ψ EP ( x )( k − k ) (16)to the leading order of the small perturbation parameter ς . Here ψ EP = ψ ± | β =1 is the coalesced eigenstate at theEP. As expected [9, 34], a second-order pole appears inthe second term due to this coalescence. Details on theperturbative analysis can be found in Appendix B.To verify the robustness of our method based on Equa-tion (3) in the vicinity of the EP, we choose the het-erostructure of length L shown in Figure 2. It con-sists of DBRs of refractive indexes n = 2 + 0 . i and n = 3 + 0 . i , and each layer accommodates aquarter of their respective wavelengths at ka = 1 . a is the lattice constant. The two half-wavelengthcavities with loss are fine-tuned to achieve an EP at k EP a = 1 . − . i using n l = 2 + 0 . i , n r =1 . . i , and k = Re[ k EP ]. We note that the lossin the right cavity is stronger than that in the left cav-ity, creating an effective PT -symmetric system with itscomplex eigenvalue spectrum shifted along the imaginaryaxis [35].Figure 2 shows the Green’s function of this systemwhen the source is placed at different locations. First, weplace the source at the center of the right cavity, wherewe expect the approximation (16) using just ψ ± ( x ) [and ψ l,r ( x )] to hold. As Figure 2(a) shows, it is nearly identicalto the result of our auxiliary eigenvalue approach (3), andthe inclusion of more CF states away from the EP barelychanges the Green’s function (not shown). If instead,we place the source at the center of the heterostructure,Equation (16) alone is insufficient to capture the featuresof the Green’s function [see Figure 2(b)]. However, once alarge set of additional eigenfunctions of L are also included,a good agreement to our auxiliary eigenvalue approach isagain observed [see Figure 2(c)]. These results show thatEquation (3) provides a reliable and convenient methodto calculate the Green’s function, even in the presence ofEPs.A similar but more striking behavior of the Green’sfunction at an EP was recently reported in an effective PT -symmetric ring cavity with refractive index [23] n ( θ ) = ( n + iδn ) + δn (cos 2 lθ + i sin 2 lθ ) . (17)Here θ is the azimuthal angle, and below we will use thearc length x = Rθ ∈ [0 , L ] as the coordinate, where R isthe radius of the ring and L = 2 πR is the circumference. n and δn are the real and imaginary parts of the back-ground index, and the latter includes both absorption andradiation losses and is positive. l is a positive integer, andthe complex index grating, proportional to e ilθ , scattersthe clockwise (CW) wave of angular momentum − l tothe counterclockwise (CCW) wave of angular momentum l but not vice versa. Note that the chirality and the sign - π π P ha s e ang l e (a) - π /2 π /2 ψ EP G ( x,x ’ ) (b) FIG. 3. Chirality-reversal Green’s function in a PT -symmetricring cavity. (a) Phase of the Green’s function (solid line) andthe coalesced eigenstate (dotted line). Dashed lines point tothe phase of the Green’s function at the source (marked by thearrow). Shaded regions show the half periods with higher loss.(b) False color plots showing the constant amplitudes of thecoalesced eigenstate (top) and the Green’s function (bottom).A finite width is imposed on the ring for visual clarity. Here n = 3, δn = 0 . l = 1, k EP L = 2 . − . i , kL =2 . x = L/ of Im[ n ] here are defined with respect to the temporaldependence e − iωt .Consequently, an EP appears at k EP R = l/ ( n + iδn )with the coalesced CCW eigenstate ψ ( x ) = e ilθ [36]. TheGreen’s function, on the other hand, can be fully decou-pled from this mode even on resonance, if the source isplaced at θ = ( m + 1 / π/l where m is a non-negativeinteger [23], i.e., at one of the most lossy spots in thepassive ring cavity. The Green’s function is given by thecorresponding Jordan vector J ( x ) ∝ e − ilθ instead, i.e.,the “missing dimension” of the Hilbert space at the EPin the CW direction.This extraordinary behavior is captured nicely using ourauxiliary eigenvalue approach (see Figure 3). In addition,a perturbative approach shows that the Green’s functionat the source is given by [23] G ( x , x ; k ) ≈ k − k EP ) k EP L , (18)which is almost imaginary on resonance (i.e., k = Re[ k EP ])for a high-Q resonance with | Im[ k EP ] | (cid:28) Re[ k EP ]. Inthe case shown in Figure 3, this value is ( − . i +0 . R and nearly identical to that given by our auxiliaryeigenvalue approach, i.e., ( − . i + 0 . R . C. Quasi-1D waveguides
Quasi-1D waveguides are frequently used in the studyof disordered mesoscopic and optical systems [24–26], andthe Green’s function plays a crucial role to construct thescattering and transfer matrices [37, 38]. Here again weconsider the scalar Helmholtz equation in a waveguide with background refractive index n ( r ): L = − n ( r ) ( ∂ x + ∂ y ) , z = k . (19)A finite width L y in the transverse direction and theDirichlet boundary conditions at y = 0 , L y lead to a setof transverse modes (“channels”) f m ( y ) = sin[ nk ( y ) m y ],where nk ( y ) m = mπ/L y is the transverse wave numberand m is a positive integer. At a given frequency k , thelongitudinal wave number in the m th channel is given by nk ( x ) m = n { k − [ k ( y ) m ] } / , and this channel is propagating(evanescent) if k ( x ) m is real (imaginary).To validate our auxiliary eigenvalue approach in quasi-1D waveguides, we first consider a uniform waveguidewith the outgoing boundary condition. In this case ananalytical expression exists for the Green’s function, whichcan be written as the following infinite sum using thechannel functions: G ( r , r (cid:48) ; k ) = X m n sin[ nk ( y ) m y ] sin[ nk ( y ) m y ] ik ( x ) m L y e ± ink ( x ) m x . (20)Here the source point r (cid:48) = ( x , y ) is placed at x = 0, y = 0 , L y , and the “+ ( − )” sign in the exponent appliesto a positive (negative) x . Clearly, only the propagatingchannels of a finite number affect the farfield behavior ofthe Green’s function, while the logarithmic divergence ofthe Green’s function at the source, a generic property in2D (including quasi-1D), is reflected by the infinite numberof evanescent channels in the summation. Numerically,this divergence is truncated either by the inclusion ofonly a finite number of evanescent channels or the finiteresolution of the spatial discretization. We also note thatthe Green’s function is dimensionless in quasi-1D. -4 -2 0 2x ( m)012 | G ( x , x ' ) | | G ( y , y ' ) | y ( m)0.5 1 1.5 FIG. 4. Green’s function in a uniform waveguide with L y =2 µ m and n = 1 . x and y at the source position x = 0 , y = 3 / µ m.The free-space wavelength is 1550 nm. Figure 4 shows one example where our auxiliary eigen-value approach, implemented using the finite differencemethod [27, 39] (solid line), is compared with the an-alytical result given by Equation (20). With just thethree propagating channels available in this case (dot-ted line), Equation (20) describes the Green’s functionwell far from the source ( > µm ); with more channelsincluded (e.g., 100; dashed line), a good agreement be-tween Equation (20) and our auxiliary eigenvalue ap-proach is observed, where the grid spacings ∆ x = 1 / µ m,∆ y = 3 / µ m used in the finite difference method arecomparable to the shortest evanescent tail in the summa-tion (i.e., 1 /k ( x ) m =100 ≈ . µ m.)For structured or disordered quasi-1D waveguidesconnected to two semi-infinite regions (“leads”), theHelmholtz equation is no longer separable in x and y ,and an analytical expression for the Green’s function inthe form of Equation (20) does not exist. Nevertheless,one can still compare to the bilinear expansion (11) as weshow in Figs. 5(a) and (b) for a PT -symmetric waveguide,despite its slow convergence [see Figure 5(c)].Here we also briefly review a standard technique usedto calculate the Green’s function in quasi-1D waveguides,i.e., the recursive Green’s function method [40, 41]. Thecomparison of the Green’s functions obtained in our pro-posed approach and by the recursive Green’s functionmethod will be present in Section III.The heart of the recursive Green’s function method liesin the celebrated Dyson’s equation: G = G + G V G . (21) G is the Green’s matrix expressed in the spatial basishere, i.e., its element G i,j gives the value of the Green’sfunction at point i when the source is placed at point j .In a quasi-1D system with N segments, V represents thecouplings between the ( N − N th segments dueto L , and G is the value of G when V is taken to bezero. In other words, G contains the Green’s matricesof two separate systems, one for the first N − y . We note, however,these two methods serve different purposes. The auxil-iary eigenvalue approach gives the values of the Green’sfunction inside the entire waveguide for a given positionof the source, which can be both in the interior of thewaveguide or on its boundaries. The recursive Green’sfunction method, on the other hand, is particularly suit-able for solving the Green’s function connecting the twoboundary layers, with the source placed in either one ofthem [See Appendix C]. Extra steps are needed to solve the Green’s function inside the waveguide, which is moreintense numerically. III. NEW INSIGHTS
Having validated the proposed auxiliary eigenvalue ap-proach to calculate the Green’s function in various sys-tems, we now turn to the new insight this method provide,i.e., viewing the Green’s function as a defect state as man-ifested by Equation (5).We start by revisiting the 1D Hermitian system de-scribed by the scalar Helmholtz equation in Figure 1.The eigenstate ψ ( x ) in our auxiliary eigenvalue equa-tion (2) (and the Green’s function) for a given energy z = k then corresponds to the wave function given by theSchr¨odinger equation inside a closed box with a δ -functionpotential: (cid:20) − n ∂ x + λ δ ( x − x ) (cid:21) ψ ( x ) = zψ ( x ) , (22) -3 0 3x ( m)02 | G ( x , x ' ) | | G ( y , y ' ) | y ( m)213 500 1000(c) 0 | Δ G | Number of CF states -3 -1 FIG. 5. (a,b) Same as Figure 4 but for the Green’s functionin a PT -symmetric waveguide. Here L x = 9 µ m, L y = 3 µ m,and the source is at x = − / µ m, y = 1 . µ m. Insetin (a): Schematic of the waveguide. Im[ n ( x )] = ± .
05 for x ∈ [ − , µ m and [0 , µ m respectively, where Re[ n ( x )] = 2in the central half of the width. Re[ n ( x )] = 1 . x = 4 . µ m, y = 1 . µ m (lower). -60 -20 20 609001000 (L -1 ) z ( L - ) (a) 0 0.5 1x/L02040 V (b) ( L - ) FIG. 6. (a) Eigenvalues of the Schr¨odinger equation (22) asa function of λ for n = 3 and x = 0 . L . Only four ofthem are shown near z = 900 /L . (b) Schematics showing thepotential as a function of the position. Shaded areas indicatethe Dirichlet boundary condition and the δ -function potentialwith λ = 38 . /L . Two unnormalized wave functions | ψ ( x ) | are also shown. and the sign of λ determines whether the potential isattractive or repulsive, which can be found analyticallyusing Equation (7).This property of λ is a crucial difference between ourapproach and previous studies of defect states, where thedefect strength (i.e., λ ) is treated as a free parameterto find all possible z ’s. In a nutshell, a given defectstrength leads to an infinite set of wave functions withdifferent energies { z i } [42], while for a given z there is asingle wave function of Equation (22) [and more generally,Equation (5)], i.e., the Green’s function, the value of whichat the source gives λ − . Such a one-to-one mappingfrom z to λ as an eigenvalue problem has not beenrecognized previously to the best of our knowledge, whichhas prevented the association of a Green’s function as adefect state in the past.To illustrate this point, we show in Figure 6(a) severalenergy eigenvalues z i ’s of Equation (22) as a functionof λ . It is important to note that the “bandwidths”of the eigenstates do not overlap [shaded areas in inFigure 6(a)], even when the range of λ is extended from −∞ to ∞ . This feature is determined by the uniquenessof the Green’s function: for a given energy z and sourceposition x , the Green’s function is uniquely determinedby the boundary condition, and so is its inverse heightat x , which gives the only allowed value of λ for thispair of z and x . Each of the “bandgaps” in Figure 6(a)shrinks to a single point when λ is extended from −∞ to ∞ , and they correspond to the energy eigenvalues ofthe two subsystems separated by the δ -potential.A particularly interesting type of defect states has aspatially localized wave function. It takes place when theenergy of the defect state is inside a bandgap, causingthe defect state to be largely decoupled from the rest ofthe system. Here we show that the Green’s function ina 1D non-Hermitian lattice can have the same property,which further stresses the new perspective our method pro-vides. This 1D lattice features a non-Hermitian flatband[15], enabled by imposing non-Hermitian particle-hole (NHPH) symmetry [43–46] via staggered gain and lossalong a lattice of optical resonators. In the tight-bindingHamiltonian description, it is given by Hφ n = [ z + ( − n − iγ ] φ n + g ( φ n +1 + φ n − ) , (23)where n labels the lattice sites from 1 to N , g ∈ R isthe nearest neighbor coupling, and the alternate positiveand negative imaginary on-site detuning γ models gainand loss. Below we shift the energy reference such thatthe on-site energy z of a resonator mode is zero. NHPHsymmetry then warrants a symmetric energy spectrumabout the origin of the complex plane, i.e., z i = − z ∗ j . Ata critical value of the gain and loss strength γ = 2 g , thereal part of the energy bands collapse to 0, with all theeigenvalues { z i } completely imaginary.Now let us consider the Green’s function of such asystem, i.e., ( z − H ) G = δ nn . (24)In the auxiliary eigenvalue problem, it can be rewrittenas ( H + λ δ nn ) ψ n = zψ n , G = ψ n λ ψ n . (25)In other words, the Green’s function corresponds to asingle defect state, where the on-site energy at the positionof the source n is shifted by λ .We emphasize that this shift λ is complex in generaland hence different from the Hermitian case shown inFigure 6 and other typical considerations of a defectstate, including that in Reference [15]. One exampleis shown in Figure 7, where we place the source at theleft edge ( n = 1) of a lattice with 60 lattice sites. λ stays in the lower complex plane when z increases from − g to g , and the stronger defect strength (and lowerLDOS) at z = ± g gives the Green’s function a clearlylocalized profile, as reflected by the half participationratio (HPR) defined by P n ( | ψ n | ) / (2 P n | ψ n | ) inthe units of the lattice constant [see Figure 7(b)].We note that the HPR of an exponentially localizedintensity | ψ n | = e − n/ξ ( n = 1 , , . . . , N ) is approx-imately equal to its localization length ξ [47] when ξ, N (cid:29)
1. The localization of the Green’s function isweakened as the defect strength λ becomes smallernear z = 0, with its HPR now greater than 15 lattice sites.Finally, we reveal the most intriguing connection be-tween the Green’s function and a defect state, i.e., ina 2D topological lattice with a chiral edge state. Thissystem breaks Lorentz reciprocity [48, 49], which can beachieved in an electronic system by imposing a magneticfield. An analogy can be introduced to photonic systemsby imposing an artificial gauge field, achieved experimen-tally by asymmetric couplings between two neighboringlattice sites on a tight-binding lattice [16–19]: while thecouplings in both the x and y directions still have thesame amplitude, their phases are now different. Here weconsider a square lattice (Figure 8) with uniform verticalcoupling and horizontal asymmetric couplings of the sameamplitude g , realizing a Landau gauge with a π/ − g is also introduced to decouple itfrom the rest of the system.Due to its sublattice symmetry [50], the energy spec-trum of the main region is symmetric about z = 0, set atthe value of the identical on-site potential. It has an edgeband with CCW chiral edge states in z/g ∈ [1 . , . z/g ∈ [ − . , − . x = a, y = 10 a , where a = 1 is the lattice constant.As Figure 8(a) shows, the Green’s function with z atthe middle of the CCW chiral edge band displays thesame characteristic as the eigenstates of L , i.e., localizednear the edge of the topological region. Our auxiliaryeigenvalue approach gives identical result as the recur-sive Green’s function method [solid and dashed lines inFigs. 8(b) and (c)], and the Green’s function displaystwo noticeable features when compared with the bilin-ear expansion: while only a small number ( ∼
20) of edgestates are needed to capture the Green’s function on theopposite (right) edge of the system (dotted line), theycannot describe the behavior of the Green’s function nearthe source, where it displays a local minimum.One way to understand this behavior is offered byour auxiliary eigenvalue approach, where the Green’sfunction is effectively a defect state perturbed by the δ -function potential [see Equation (5)]. In this point ofview, the local minimum of the Green’s function is theconsequence of the topological protection of the chiraledge states: it circumvents this defect as can be clearlyseen in Figure 8(a). This change, of course, creates alocal perturbation so strong that can only be capturedby a large number of the unperturbed eigenstates of L in the bilinear expansion, including even the chiral edgestates in the opposite direction that are much different inenergy (not shown). The far-field, on the other hand, isnot affected by this local perturbation, as we have seenin Figure 8(b). We note that the Green’s function does -1 0 1 z/g H P R -1 0 1 Re[ /g] -1.4-1.2-1 I m [ / g ] increasing z (a) (b) FIG. 7. (a) Complex defect strength λ when z increases from − g to g . The lattice has 60 sites and γ = 2 g . (b) HPR asa function of z , showing a strong localization (HPR ∼
1) at | z | = t to a weaker localization (HPR >
15) at z = 0. not display a local minimum at all energies: if z is almoston resonance with an unperturbed eigenstate in the edgeband, the Green’s function is given essentially by thissingle chiral edge mode, which goes through the defectcreated by the δ -function potential with little scattering. IV. DISCUSSION AND CONCLUSIONS
In summary, we have introduced in this work a new per-spective of the time-independent Green’s function, whichis obtained as a single eigenstate of an auxiliary eigen-value formulation that embodies a defect state createdby the δ -function potential. The height of the δ -functionpotential is determined by the inverse value of the Green’sfunction at the source position, which is given directlyin the form of a generalized eigenvalue problem givenby Equation (2). It is the only well behaved and finite Φ = π /2 region(a) x/a y / a y/a | G ( y , y ' ) | (b) y/a | G ( y , y ' ) | (c) maxmin Φ = - π /2 region FIG. 8. Analogy between the Green’s function and a chiraledge state circumventing a defect in a topological squarelattice. (a) The Green’s function calculated with z = 1 . g and the source at x = a, y = 10 a . Solid, dashed, and dottedlines show the results of our auxiliary eigenvalue approach,the recursive Green’s function, and the bilinear expansion,respectively. (b,c) The values of the Green’s function along y at x = 40 a and a , respectively. eigenvalue, easy to identify numerically. Therefore, ourapproach differs both conceptually and computationallyfrom previous investigations of eigenstates of a δ -functionpotential, which were not related to the Green’s function.The uniqueness of the eigenstate that gives the Green’sfunction in our approach should be distinguished from thesingle bound state in an attractive 1D δ -function potential[51]. In our case the δ -function potential can be repulsiveor attractive in a Hermitian system [see Figure 6(a)], andit becomes complex in general in non-Hermitian systems.Furthermore, the eigenstate in our approach does notdepend on the original potential included in the oper-ator L , while the aforementioned bound state assumesa vacuum background. Finally, our eigenstate exists inhigher dimensions as well, with clearly different boundaryconditions and spatial profiles from a bound state [seeFigure 8(a), for example].We have also verified the Green’s function obtainedin our method by comparing both to analytical resultswhen available and to two other frequently used numeri-cal methods, i.e., the recursive Green’s function methodand the bilinear expansion in the basis of the system’seigenstates. At an EP, where a perturbative treatment ofthe bilinear expansion becomes necessary, our method isstill robust as seen from the examples in two PT -systems.Our method also gives identical result to the recursiveGreen’s function method, implemented by finite differenceon the same tight-binding lattice.Our defect-state approach can also be applied to morenumerically demanding cases, e.g., in the study of diffu-sive transport and wave localization [24, 25]. We haveinvestigated disordered quasi-1D waveguides over 100wavelengths long and with over 60 transverse channels(not shown). In these cases, the memory storage for therecursive Green’s function may become an issue, whenthe values of the Green’s function in the interior of thewaveguides are also computed for light-matter interactionor laser emissions. This is because the Green’s matrixobtained from the recursive procedure mingles the valuesof the Green’s function generated by sources placed acrossthe entire system. The defect-state approach excels inthis regard, because each source is treated independently.We should also mention that this advantage becomes lessnoticeable or even a disadvantage if the Green’s functionneeds to be evaluated with the source at many differentlocations.More importantly, our approach offers a previously un-explored physical insight that both the recursive Green’sfunction method and the bilinear expansion lack, i.e., thelinkage between the Green’s function and a defect state.We have exemplified an intriguing manifestation of thislinkage using a topological chiral edge state, where thelocal minimum of the Green’s function is analogous to achiral edge state circumventing a boundary defect. There-fore, even though our discussions have focused on theHelmholtz equation for scalar optical waves, we expectthis new perspective and the physical insight it offers tohave important applications in other related fields as well, including condensed matter systems, acoustics, electroniccircuits and so on. Acknowledgements
We thank Douglas Stone, Yidong Chong and AzrielGenack for helpful discussions. This project is supportedby the NSF under Grant No. PHY-1847240.
Appendix A: Reciprocity in the auxiliary eigenvalueapproach
We derive some relations between the Green’s functions G ( r , r (cid:48) ) and G ( r (cid:48) , r ), based on the symmetries of L . Westart from a matrix description in which the Green’sfunction operator G ( z ) is defined from G ( r , r (cid:48) ; z ) = h r | G ( z ) | r (cid:48) i . We use the same symbol L for the matrix rep-resentation of the system. The Green’s function satisfiesthe following equation:[ z − L ] G ( z ) = , G ( z )[ z − L ] = , (A1)where is the identity matrix. Here, we identify theGreen’s function as being both the right and left matrixinverses of [ z − L ]. Considering that in non-Hermitianphysics, there is a distinction between left and right sets ofeigenvectors [52], one may attempt to distinguish betweenthe left and right Green’s functions. However, it can beshown that they are identical, since [ z − L ] is a matrixof full rank. It is possible to obtain relations between G ( z ) and [ G ( z )] T , or [ G ( z )] † , when performing matrixtransposition or hermitian conjugation on Equation (A1).We start by analyzing the case L T = L : matrixtransposition of Equation (A1) leads to the relation G ( z ) T = G ( z ), or more explicitly G ( r (cid:48) , r ; z ) = G ( r , r (cid:48) ; z ),which is the usual Lorentz reciprocity condition. Thisrelation holds when L is a symmetric matrix, describinga system with symmetric couplings and hence withoutan effective magnetic field. It is valid for non-Hermitiansystems with this property as well, such as the case witha non-Hermitian flatband discussed in Figure 7 of themain text.In a general Hermitian system, L † = L , and uniquenessof the matrix inverse demands that [ G ( z )] T = [ G ( z ∗ )] ∗ , or G ( r (cid:48) , r ; z ) = [ G ( r , r (cid:48) ; z ∗ )] ∗ . It is customary to state thisrelation in terms of Green’s functions G ( ± ) with outgoing and incoming boundary conditions, respectively. Thisnomenclature is specially useful if z = k + iη , with η > k : G ( ± ) ( r , r (cid:48) ; k ) ≡ lim η → + G ( r , r (cid:48) ; k ± iη ) . (A2)Therefore, we can write the reciprocity condition as G ( ± ) ( r (cid:48) , r ; k ) = [ G ( ∓ ) ( r , r (cid:48) ; k )] ∗ . The solutions G ( ± ) are oftentimes referred to as the retarded and advanced Green’s functions, inspired from the association to sourcesgenerating outgoing or incoming waves. Systems that arenon-Hermitian and do not have symmetric couplings lead0to more complicated reciprocity relations, involving theGreen’s functions of two different systems. This is ex-pected, since the scattering properties derived from theGreen’s functions of systems related by a time-reversaltransformation are nontrivial and yet to be understood.Reference [53] provides interesting scattering relationsbetween two “time-reversal partners” in 1D and quasi 1Dsystems.Our Green’s function approach for time-reversal or non-Hermitian symmetric systems lead to a relation betweenthe single eigenstates: ψ ( r , r (cid:48) ) λ ( r (cid:48) ) ψ ( r (cid:48) , r (cid:48) ) = ψ ( r (cid:48) , r ) λ ( r ) ψ ( r , r ) , (A3)or ψ ( r (cid:48) , r ) = ψ ( r , r (cid:48) ) (cid:20) λ ( r ) ψ ( r , r ) λ ( r (cid:48) ) ψ ( r (cid:48) , r (cid:48) ) (cid:21) , (A4)whereas for the general non-Hermitian problem, we find¯ ψ ( r (cid:48) , r ) = [ ψ ( r , r (cid:48) )] ∗ (cid:20) ¯ λ ( r ) ¯ ψ ( r , r ) λ ( r (cid:48) ) ψ ( r (cid:48) , r (cid:48) ) (cid:21) , (A5)where we used the notation ¯ ψ ( r , r (cid:48) , z ) ≡ ψ ( r , r (cid:48) , z ∗ ).The same result is obtained if we identify ¯ ψ with theadvanced solution and ψ with the retarded solution, orvice-versa. Appendix B: Perturbation theory at an EP
Below we review a generic procedure to approach effec-tive PT -symmetric systems with the outgoing boundarycondition based on the coupled-mode theory. In par-ticular, we describe a photonic molecule that consistson identical components L and R that support localizedstates ψ l ( x ) and ψ r ( x ) in each component. The systemcan be, for instance, two identical whispering gallery res-onators that are brought to proximity, such that theirmodes overlap and produce a pair of bonding and anti-bonding states. Alternatively, one can think of a pair ofhalf-wavelength cavities coupled by a DBR, which is thesystem discussed in Section 2.2 of the main text. Notethat due to the outgoing boundary condition, such sys-tems are not truly PT -symmetric, because time reversalchanges the outgoing boundary condition to the incomingboundary condition. However, if the states are highlylocalized as in our discussion below, the minute flux in orout of the system can be omitted in the discussion of thecoupled mode theory.The left component satisfies the Helmholtz equation: (cid:2) ∂ x + (cid:15) l ( x ) k (cid:3) ψ l ( x ) = 0 , ( x ∈ L ) (B1) (cid:2) ∂ x + k (cid:3) ψ l ( x ) = 0 . ( x ∈ elsewhere) (B2)Here k is the real-valued external frequency at which theGreen’s function is going to be evaluated, and k ∼ k is one complex-valued CF frequency. The right component isdefined similarly, and the localized states are normalizedby R x ∈ L , R (cid:15) l,r ψ l,r dx = 1, respectively. The compositesystem is defined by (cid:15) ( x ) = (cid:15) l ( x ) + (cid:15) r ( x ) and satisfies (cid:2) ∂ x + (cid:15) ( x )˜ k (cid:3) ˜ ψ ( x ) = 0 , ( x ∈ L ∪ R ) (B3) (cid:2) ∂ x + k (cid:3) ˜ ψ ( x ) = 0 , ( x ∈ elsewhere) (B4)where ˜ k is the CF frequency of a “supermode”˜ ψ ( x ) = a l ψ l ( x ) + a r ψ r ( x ) (B5)to be determined. We note that ˜ ψ ( x ) automatically satis-fies the outgoing boundary condition at the two ends ofthe composite system due to the same property of both ψ l,r ( x ).We make the following assumption about the system:the spatial overlapping of the localized states in one cav-ity is very weak, i.e., R x ∈ L , R (cid:15) l,r ψ l ψ r dx ≡ g ( | g | (cid:28) | R x ∈ R (cid:15) r ψ l dx | = | R x ∈ L (cid:15) l ψ r dx | is even smaller, which isthen neglected in our analysis below.To derive the coupled-mode theory, we first insert Equa-tion (B5) into Eqs. (B3), (B4) and simplify the resultusing Eqs. (B1), (B2) for ψ l ( x ) and their counterparts for ψ r ( x ). Then, we multiply the resulting equation by either ψ l ( x ) or ψ r ( x ) and perform an integration in the left andright component of the system respectively, obtaining inthis way two independent equations: (cid:20) ˜ k − k g ˜ k − tk g ˜ k − tk ˜ k − k (cid:21) (cid:20) a l a r (cid:21) = 0 . (B6)It is a system of equations for the unknown amplitudes a l,r and frequency ˜ k , where t is defined by R x ∈ L ψ r dx = R x ∈ R ψ l dx . Since the matrix operator is symmetric andhas identical diagonal elements, a zero-determinant condi-tion implies that a r = ± a l , giving rise to the symmetric and anti-symmetric modes ˜ ψ S ( x ) and ˜ ψ A ( x ) respectively,with frequencies: ˜ k S,A = k ± tk ± g . (B7)With the aforementioned approximations | g | (cid:28) k ∼ k , these frequencies reduce to˜ k S ≈ (1 + t − g ) k , ˜ k A ≈ (1 − t + g ) k , (B8)and their splitting equal to ˜ k S − ˜ k A = 2( t − g ) k ≡ k .Now, we introduce a PT -symmetric perturbation i(cid:15) ( x ), which is an odd function of x . The effective PT -symmetric system satisfies[ ∂ x + ( (cid:15) ( x ) + i(cid:15) ( x )) q ] ψ ( x ) = 0 , (B9)where q is the perturbed CF frequency to be determinedtogether with the amplitudes A and S in the expansionof the perturbed eigenstate ψ ( x ) = A ˜ ψ A ( x ) + S ˜ ψ S ( x ).1Following an analogous procedure to the one describedabove, we obtain again a system of equations that can berepresented by the following matrix equation for A , S : (cid:20) q − ˜ k S iCq iCq q − ˜ k A (cid:21) (cid:20) SA (cid:21) = 0 , (B10)where C ≡ R (cid:15) ( x ) ˜ ψ A ( x ) ˜ ψ S ( x ) dx is integrated over thewhole system. We note that this equation and the follow-ing perturbative results apply when (cid:15) ( x ) is complex, i.e.,the system does not need to have balanced gain and lossto be effectively PT symmetric [35].To rewrite it in a more familiar form, we note that sincethe PT -symmetric perturbation is assumed to be weak,we approximate Cq by Ck that leads to˜ H (cid:20) SA (cid:21) ≡ (cid:20) k S − iCk − iCk k A (cid:21) (cid:20) SA (cid:21) = q (cid:20) SA (cid:21) . (B11)The perturbed frequency q is then found to be: q ± ≈ k (cid:16) ± p ∆ − C (cid:17) , (B12) where we have used the forms of ˜ k S,A from Equation (B8).The EP then emerges when the C = ∆ , where q ± ≡ k EP ≈ k . In the main text we have introduced ς ≡ p − β and β ≡ C/ ∆.This system of equations can be transformed to thebasis of left- and right-components, in terms of the am-plitudes a l , a r : H (cid:20) a l a r (cid:21) ≡ k (cid:20) − iC ∆∆ 1 + iC (cid:21) (cid:20) a l a r (cid:21) = q (cid:20) a l a r (cid:21) . (B13)It has the most familiar form of a PT -symmetric Hamil-tonian. This is the model that we have used to describethe photonic molecule in Section 2.2.The perturbative result of the Green’s function givenby Equation (16) in the main text is derived from thefollowing expression: G ( x, x ; k ) = (cid:15) ( x )( k − k ) − ∆ ς { ( k − k )[ ψ l ( x ) ψ l ( x ) + ψ r ( x ) ψ r ( x )] − ∆ k [ iβ [ ψ l ( x ) ψ l ( x ) − ψ r ( x ) ψ r ( x )] − [ ψ l ( x ) ψ r ( x ) + ψ r ( x ) ψ l ( x )] ] } (B14)by dropping terms proportional to ς , including that in β = √ − ς ≈
1. We note that Equation (B14) includesterms to all orders of ς . Here it is unnecessary to invokethe Jordan vector J ( x ) that completes the Hilbert spaceat the EP [23], but to compare with the result given inReference [10], we note that it can be chosen as J ( x ) = iψ l ( x ) + ψ r ( x )2 √ ∆ k e − iπ/ , (B15)which satisfies [ H − k ] J ( x ) = ψ EP ( x ) (B16)with the coalesced eigenstate ψ EP ( x ) = 2 √ ∆ k e − iπ/ [ ψ l ( x ) + iψ r ( x )] (B17)at the EP. Here is the identity matrix, and ψ EP ( x )is normalized differently from the main text to satisfy( ψ EP , J ) = 1 [see the definition of the inner product givenby Equation (14) in the main text]. We further note thatthe PT -symmetric Hamiltonian given by Equation (B13)above is symmetric, and hence it is unnecessary to dis-tinguish left and right eigenstates (and Jordan vectors)because they are identical. Using J ( x ) and ψ EP ( x ) speci-fied above, we find that the Green’s function derived inReference [10] at the EP agrees with our result given by Equation (16) in the main text, once the additional factor (cid:15) ( x ) is accounted for that originates from rewriting theHelmholtz equation as Equation (6) in the main text. Appendix C: The recursive Green’s function
Below we briefly review the recursive Green’s functionmethod and our implementation of it using finite differ-ence. Assuming that the Green’s function defined byEquation (1) has been obtained at every grid point for aquasi-1D waveguide with N − N y rows, wecan write down a Green’s matrix G ( N − that satisfies h z ( N − − L ( N − i G ( N − = ( N − . (C1)Its individual element G ( N − i,j gives the value of theGreen’s function at point i when the source is placedat point j . We reserve the notation G ( N − i,j without theunderscore below the subscripts for the N y × N y block ofthe Green’s matrix evaluated in the i th column with thesource in the j th column, and we apply the same notationto other matrices as well. The ( N − N y grid points arelabeled from the first to the last point in the first column,continued in the same way onto the second and remainingcolumns. ( N − is the identity matrix with ( N − N y L ( N − is the finite differenceimplementation of the operator L in Equation (1), againwith ( N − N y rows and columns.To calculate the Green’s matrix with an additionalcolumn added from the right, we separate L ( N ) into L + V and construct an ancillary matrix G , where L ≡ L ( N − L ( N ) N,N ! , G ≡ (cid:18) G ( N − G N,N (cid:19) . V only contains the couplings between the ( N − N th columns in L ( N ) , and h z (1) − L ( N ) N,N i G N,N = (1) (C2)defines the Green’s matrix G N,N of the isolated N thcolumn. It is straightforward to show that h z ( N ) − L i G = ( N ) (C3)and we derive the well known Dyson’s equation G ( N ) = G + G ( N ) V G , (C4)where G ( N ) is defined by h z ( N ) − L ( N ) i G ( N ) = ( N ) (C5)similar to Equation (C1). Equation (C4) can be verifieddirectly by multiplying [ z ( N ) − L ( N ) ] from the left. Below we drop the superscript of G ( N ) . Equation (C4)differs from its more familiar form in mesoscopic physics[41], i.e., ˜ G = ˜ G + ˜ G V ˜ G , (C6)because there the Green’s matrix is defined by ˜ G [ z ( N ) − L ( N ) ] = ( N ) instead. We have already proven the equiva-lence between the left and right matrix inverses of [ z −L ]in Appendix A, thus ˜ G and G are identical.Using Equation (C4), the recursive relations betweenthe blocks of the Green’s matrices G i,j connecting thetwo boundary layers (i.e., the 1st and N th columns) canbe derived: G NN = A − , (C7) G N = G ,N − V N − ,N A − , (C8) G N = G N,N V N,N − G N − , , (C9) G = G , + G ,N V N,N − G N − , . (C10)A single matrix inversion for A = z (1) − L ( N ) N,N − V N,N − G N − ,N − V N − ,N (C11)is needed when adding one layer from the right, andits last term is often referred to as the self-energy [40].Note that besides the operator L ( N ) , all the pieces ofinformation needed for this recursive scheme are the samefour matrix for the system with ( N −
1) layers. As a result,the recursive Green’s function method is highly efficientto treat transport problems, which does not require theknowledge of the Green’s function inside the system. If weneed the full Green’s matrix to investigate, for example,light-matter interaction in the bulk of the system, moresteps are required in the recursive Green’s function, whichscale as N and are more intense numerically. [1] G. Green, Mathematical papers of the late George Green (Macmillan, London, 1871).[2] I. Stakgold and M. Holst,
Green’s functions and boundaryvalue problems , 3rd ed., Pure and Applied Mathematics(John Wiley & Sons, Hoboken, 2011).[3] Y. A. Melnikov and V. N. Borodin,
Green’s Functions ,Developments in Mathematics No. 48 (Springer, Cham,2017).[4] P. M. Morse and H. Feshbach,
Methods of TheoreticalPhysics , Vol. 1 (McGraw-Hill, New York, 1953).[5] J. Schwinger,
Particles, Sources, and Fields (Perseus,Reading, Mass, 1998).[6] G. Rickayzen,
Green’s functions and Condensed Matter (Academic Press, New York, 1980).[7] E. Economou,
Green’s functions in Quantum Physics , 3rded. (Springer, 2007).[8] M. Davy, Z. Shi, J. Wang, X. Cheng, and A. Z. Genack,Transmission Eigenchannels and the Densities of Statesof Random Media, Phys. Rev. Lett. , 033901 (2015).[9] Z. Lin, A. Pick, M. Lonˇcar, and A. W. Rodriguez, En- hanced Spontaneous Emission at Third-Order Dirac Ex-ceptional Points in Inverse-Designed Photonic Crystals,Phys. Rev. Lett. , 107402 (2016).[10] A. Pick, B. Zhen, O. D. Miller, C. W. Hsu, F. Hernandez,A. W. Rodriguez, M. Soljaˇci´c, and S. G. Johnson, Generaltheory of spontaneous emission near exceptional points,Opt. Express , 12325 (2017).[11] R. G. Newton, Scattering Theory of Waves and Particles ,2nd ed. (Springer, Berlin, 1982).[12] H. E. T¨ureci, A. D. Stone, and B. Collier, Self-consistentmultimode lasing theory for complex or random lasingmedia, Phys. Rev. A , 043822 (2006).[13] H. E. T¨ureci, L. Ge, S. Rotter, and A. D. Stone, Stronginteractions in multimode random lasers, Science ,643 (2008).[14] L. Ge, Y. D. Chong, and A. D. Stone, Steady-state abinitio laser theory: Generalizations and analytic results,Phys. Rev. A , 063824 (2010).[15] B. Qi, L. Zhang, and L. Ge, Defect States Emerging froma Non-Hermitian Flatband of Photonic Zero Modes, Phys. Rev. Lett. , 093901 (2018).[16] M. Hafezi, E. A. Demler, M. D. Lukin, and J. M. Taylor,Robust optical delay lines with topological protection,Nat. Phys. , 907 (2011).[17] M. Hafezi, S. Mittal, J. Fan, A. Migdall, and J. M. Taylor,Imaging topological edge states in silicon photonics, Nat.Photonics , 1001 (2013).[18] M. A. Bandres, S. Wittek, G. Harari, M. Parto, J. Ren,M. Segev, D. N. Christodoulides, and M. Khajavikhan,Topological insulator laser: Experiments, Science ,eaar4005 (2018).[19] H. Zhao, X. Qiao, T. Wu, B. Midya, S. Longhi, andL. Feng, Non-Hermitian topological light steering, Science , 1163 (2019).[20] S. Barik, A. Karasahin, C. Flower, T. Cai, H. Miyake,W. DeGottardi, M. Hafezi, and E. Waks, A topologicalquantum optics interface, Science , 666 (2018).[21] C. M. Bender and S. Boettcher, Real spectra in non-Hermitian Hamiltonians having PT symmetry, Phys. Rev.Lett. , 5243 (1998).[22] L. Feng, R. El-Ganainy, and L. Ge, Non-Hermitian pho-tonics based on parity–time symmetry, Nat. Photonics , 752 (2017).[23] H. Z. Chen, T. Liu, H. Y. Luan, R. J. Liu, X. Y. Wang,X. F. Zhu, Y. B. Li, Z. M. Gu, S. J. Liang, H. Gao, L. Lu,L. Ge, S. Zhang, J. Zhu, and R. M. Ma, Revealing themissing dimension at an exceptional point, Nat. Phys. ,571 (2020).[24] C. W. J. Beenakker, Random-matrix theory of quantumtransport, Rev. Mod. Phys. , 731 (1997).[25] Y. Imry, Introduction to Mesoscopic Physics , 2nd ed. (Ox-ford University Press, Oxford, 2008).[26] P. A. Mello and N. Kumar,
Quantum Transport in Meso-scopic Systems (Oxford University Press, 2004).[27] L. Ge, Constructing the scattering matrix for opticalmicrocavities as a nonlocal boundary value problem, Pho-tonics Res. , B20 (2017).[28] M. Liertzer, L. Ge, A. Cerjan, A. D. Stone, H. E. T¨ureci,and S. Rotter, Pump-induced exceptional points in lasers,Phys. Rev. Lett. , 173901 (2012).[29] W. D. Heiss, Exceptional points of non-Hermitian opera-tors, J. Phys. A , 2455 (2004).[30] M. V. Berry, Physics of nonhermitian degeneracies, Czech.J. Phys. , 1039 (2004).[31] L. Ge and R. El-Ganainy, Nonlinear modal interactionsin parity-time (PT) symmetric lasers, Sci. Rep. , 24889(2016).[32] P. L. Kapur and R. Peierls, The dispersion formula fornuclear reactions, P. Roy Soc. A-Math. Phy. , 277(1937).[33] M. L. Goldberger and K. M. Watson, Lifetime and decayof unstable particles in S-matrix theory, Phys. Rev. ,B1472 (1964).[34] W. D. Heiss, Green’s Functions at Exceptional Points,Int. J. Theor. Phys. , 3954 (2015).[35] A. Guo, G. J. Salamo, D. Duchesne, R. Morandotti,M. Volatier-Ravat, V. Aimez, G. A. Siviloglou, and D. N. Christodoulides, Observation of PT-symmetry breakingin complex optical potentials, Phys. Rev. Lett. , 93902(2009).[36] P. Miao, Z. Zhang, J. Sun, W. Walasik, S. Longhi, N. M.Litchinitser, and L. Feng, Orbital angular momentummicrolaser, Science , 464 (2016).[37] D. S. Fisher and P. A. Lee, Relation between conductivityand transmission matrix, Phys. Rev. B , 6851 (1981).[38] A. D. Stone and A. Szafer, What is measured when youmeasure a resistance? The landauer formula revisited,IBM J. Res. Dev. , 384 (1988).[39] L. Ge, K. G. Makris, D. N. Christodoulides, and L. Feng,Scattering in PT- and RT-symmetric multimode waveg-uides: Generalized conservation laws and spontaneoussymmetry breaking beyond one dimension, Phys. Rev. A , 062135 (2015).[40] D. J. Thouless and S. Kirkpatrick, Conductivity of thedisordered linear chain, J. Phys. C , 235 (1981).[41] F. Sols, M. Macucci, U. Ravaioli, and K. Hess, Theoryfor a quantum modulated transistor, J. Appl. Phys. ,3892 (1989).[42] S. H. Patil, Harmonic oscillator with a $ \ updelta$-function potential, Eur. J. Phys. , 899 (2006).[43] S. Malzard, C. Poli, and H. Schomerus, Topologically Pro-tected Defect States in Open Photonic Systems with Non-Hermitian Charge-Conjugation and Parity-Time Symme-try, Phys. Rev. Lett. , 200402 (2015).[44] L. Ge, Symmetry-protected zero-mode laser with a tun-able spatial profile, Phys. Rev. A , 023812 (2017).[45] J. D. H. Rivero and L. Ge, Chiral symmetry in non-hermitian systems: Product rule and clifford algebra,Phys. Rev. B , 014111 (2021).[46] K. Kawabata, S. Higashikawa, Z. Gong, Y. Ashida, andM. Ueda, Topological unification of time-reversal andparticle-hole symmetries in non-Hermitian physics, NatureCommun. , 297 (2019).[47] F. Baboux, L. Ge, T. Jacqmin, M. Biondi, E. Ga-lopin, A. Lemaitre, L. Le Gratiet, I. Sagnes, S. Schmidt,H. E. Tureci, A. Amo, and J. Bloch, Bosonic Condensa-tion and Disorder-Induced Localization in a Flat Band,Phys. Rev. Lett. , 066402 (2016).[48] L. D. Landau and E. M. Lifshitz, Electrodynamics ofContinuous Media (Pergamon Press, Oxford, 1960).[49] L. Ge and L. Feng, Optical-reciprocity-induced symmetryin photonic heterostructures and its manifestation in scat-tering PT-symmetry breaking, Phys. Rev. A , 043836(2016).[50] M. Z. Hasan and C. L. Kane, Colloquium: Topologicalinsulators, Rev. Mod. Phys. , 3045 (2010).[51] K. Gottfried and T.-M. Yan, Quantum Mechanics: Fun-damentals , 2nd ed. (Springer, New York, 2003).[52] J. D. H. Rivero and L. Ge, Pseudochirality: A Manifes-tation of Noether’s Theorem in Non-Hermitian Systems,Phys. Rev. Lett. , 083902 (2020).[53] J. D. H. Rivero and L. Ge, Time-reversal-invariant scalingof light propagation in one-dimensional non-hermitiansystems, Phys. Rev. A100