The resonant nonlinear scattering theory with bound states in the radiation continuum and the second harmonic generation
aa r X i v : . [ m a t h - ph ] M a r The resonant nonlinear scattering theory with bound states in the radiationcontinuum and the second harmonic generation
R´emy F. Ndangali and Sergei V. ShabanovDepartment of Mathematics, University of Florida, Gainesville, FL 32611, USA
A nonlinear electromagnetic scattering problem is studied in the presence of bound states in theradiation continuum. It is shown that the solution is not analytic in the nonlinear susceptibility andthe conventional perturbation theory fails. A non-perturbative approach is proposed and appliedto the system of two parallel periodic arrays of dielectric cylinders with a second order nonlinearsusceptibility. This scattering system is known to have bound states in the radiation continuum.In particular, it is demonstrated that, for a wide range of values of the nonlinear susceptibility, theconversion rate of the incident fundamental harmonic into the second one can be as high as 40%when the distance between the arrays is as low as a half of the incident radiation wavelength. Theeffect is solely attributed to the presence of bound states in the radiation continuum.
I. INTRODUCTION
A conventional approach to nonlinear electromagnetic scattering problems is based on the power series expansionin a nonlinear susceptibility χ c . For example, for the 2nd order susceptibility, the physical parameter that determinesnonlinear effects is χ c E r ≪ E r is the electric field at the scattering structure. The smallness of χ c E r justifiesthe use of perturbation theory and the solution is analytic in χ c . The situation is different if the scattering structurehas resonances.Planar periodic structures (e.g., gratings) are known to exhibit sharp scattering resonances when illuminated byelectromagnetic waves (for a review see, e.g., [1, 2]). Furthermore, it is known (see, e.g., [2]) that in such structuresa local electromagnetic field E r is amplified if the structure has narrow resonances: E r ∼ E i / √ Γ, where E i is theamplitude of the incident wave and Γ is the width of the resonance. Consequently, optical nonlinear effects areamplified if the system has a sufficiently narrow resonance. For example, an amplification of the second harmonicsgeneration by a single periodic array of dielectric cylinders [3] and by other single-array periodic systems [4] have beenreported, but no significant flux conversion rate, comparable to that in conventional methods of second harmonicgeneration, has been found. Owing to the smallness of χ c and a finite / √ Γ, the condition χ c E r ≪ h , then it can be shownthat for each resonance associated with the single structure, the combined structure has two close resonances whosewidth depends continuously on h so that the width of one of the resonances vanishes, i.e., Γ( h ) → h → h b forsome discrete set of distances h b , for sufficiently large h [2, 5, 6]. This means that the system has bound states in theradiation continuum . Their existence was first predicted in quantum mechanics by von Neumann and Wigner [7] in1929 and later they were discovered in some atomic systems [8] (see also [9–11] for more theoretical studies). Theiranalog in Maxwell’s theory has only attracted attention recently [5, 12–14]. In particular, for a system of two parallelarrays of periodically positioned subwavelength dielectric cylinders (depicted in Panel (a) of Fig. 1), the existence ofbound states in the radiation continuum has been first established in numerical studies of the system [5]. A completeclassification of bound states as well as their analytic form for this system is given in [6] for TM polarization. It is alsoshown [6] that bound states exist in the spectral range in which more than one diffraction channel are open. Fromthe physical point of view, bound states in the radiation continuum are localized solutions of Maxwell’s equations likewaveguide modes, but in contrast to the latter their spectrum lies in the spectrum of scattering radiation (diffraction)modes.The perturbation theory parameter χ c E r ∼ χ c E i / p Γ( h ), as defined above, can no longer be considered small ifbound states in the radiation continuum are present. This qualitative assessment should be taken with a precaution.In the present study, a rigorous analysis of the nonlinear scattering problem by means of the formalism of Siegertstates (appropriately extended to periodic structures [15]) shows that no divergence of a local field occurs as h → h b .However, the conventional perturbative approach fails because the solution is not analytic in χ c . The situation canbe compared with a simple mechanical analog. Consider a scattering problem for a particle on a line in a hard corerepulsive potential V ( x ) = g/x , g >
0. No matter how small g is, the particle never crosses the origin x = 0 and afull reflection occurs, but it does so when g = 0 (a full transmission). So, the scattering amplitude is not analytic in g . Other, more sophisticated, examples of quantum systems with such properties are studied in [16].The purpose of the present study is twofold. First, the nonlinear scattering problem is studied in the presence ofbound states in the radiation continuum . A non-perturbative approach is developed to solve the problem. Second,as an application of the developed formalism, the problem of the second harmonic generation is analyzed with anexample of the system depicted in Fig. 1 (Panel (a)). A. An overview of the results
In Section II, the nonlinear resonant scattering problem with bound states in the radiation continuum is trans-formed into a system of integral equations. A non-perturbative method is proposed to solve these equations in theapproximation that takes into account two nonlinear effects: a second harmonic generation in the leading order of χ c ,and the fundamental harmonic generation by mixing the second and fundamental harmonics in the leading order in χ c . This is the second order effect in χ c known in the theory as the optical rectification . The latter is shown to benecessary to ensure the energy flux conservation.The formalism is illustrated with an example of two parallel periodic arrays of dielectric cylinders shown in Panel(a) of Fig. 1. The analysis is based on the subwavelength approximation (Section III) when the incident wave lengthis larger than the radius R of the cylinders. If k is the magnitude of the wave vector, then the theory has three smallparameters: δ ( k ) = ( kR ) ( ε c − ≪ , χ c ≪ , | ∆ h | = | h − h b | ≪ R < / δ ( k ) is the scatteringphase for a single cylinder, ε c is the linear dielectric susceptibility, and the amplitude of the incident wave is set toone, E i = 1. With this choice of units, all three parameters are dimensionless. The scattering amplitudes of thefundamental and second harmonics are explicitly found in Section IV.In Section V, it is shown that the ratio of the flux of the second harmonic along the normal direction and theincident flux is σ = Cχ c E r where E r = E r ( χ c , ∆ h, δ ) is the fundamental field on the cylinders, and C = C ( δ , ∆ h ) is some function. The function E r is non-analytic in the vicinity of zero values of its arguments. The non-perturbative approach of Section II is usedto prove that the generated flux of the second harmonic attains its maximal value when the small parameters satisfythe condition (∆ h ) δ ( k b ) = ϑ χ c (1)where ϑ ≪ k b is the magnitude of the wave vector of the bound state that occurs at h = h b . Under this condition, σ becomes analytic in the scattering phase δ so that in the leading order, σ , max ≈ πk − b δ ( k b )An interesting feature to note is the independence of the conversion efficiency on the nonlinear susceptibility χ c (inthe leading order in the scattering phase δ ). In other words, given a nonlinear susceptibility χ c , by a fine tuning ofthe distance between the arrays one can always reach the maximal value which is only determined by the scatteringphase at the wave length of a bound state. The lowest value of k b for the system considered occurs just below thefirst diffraction threshold (the wavelength is slightly larger than the structure period) [6], i.e., k b ≈ π . Taking, forexample, R = 0 .
15 and ε c = 2 (so that δ (2 π ) ≈ . σ , max ≈ .
44, that is, about 44% ofthe incident flux is converted into the second harmonic flux, which is comparable with the conversion rate achievedin slabs (crystals) of optical nonlinear materials [17].From the physical point of view, the scattering structure plays the role of a resonator with the quality factor inverselyproportional to Γ. The field in the resonator is not uniform and has periodic peaks of the amplitude E r ∼ E i / √ Γdue to a constructive interference of the scattered fundamental harmonic. The second harmonic is produced by theinduced dipole radiation of point scatters located at these peaks. The induced dipole strength is proportional to χ c E r . The dipoles are excited by the incident wave and, due to their periodic arrangement, they radiate in phaseproducing a plane wave in the asymptotic region (just like a phased array antenna). If the system has a resonancewhose width Γ can be continuously driven to zero by changing a physical parameter of the system, i.e., the system hasa bound state in the radiation continuum, then the strength of the induced dipoles radiating the second harmonicscan be magnified as desired, but the resonator cannot be excited by the incident radiation if Γ = 0 (a bound stateis decoupled from the radiation continuum). So, the optimal width Γ at which the second harmonic amplitude ismaximal occurs for some Γ = 0, which explains the existence of conditions like (1). Since the second harmonic isgenerated by point scatterers, the phase matching condition, needed for optically nonlinear crystals, is not required.The energy flux of the incident radiation is automatically redistributed and focused on the scatterers owing to theconstructive interference. Thanks to these physical features, an active length at which the conversion rate is maximalis close to 2 h b whose smallest value for the system studied is roughly a half of the wave length of the incident light[5, 6] (i.e. for an infrared incident radiation it is about a few hundreds nanometers). II. THE NONLINEAR RESONANT SCATTERING THEORY
Suppose that a scattering system has a translational symmetry along a particular direction and has non-dispersivelinear and second-order nonlinear dielectric susceptibilities, ε and χ , respectively. When the electric field is parallelto the translational symmetry axis (TM Polarization), Maxwell’s equations are reduced to the scalar nonlinear waveequation 1 c ∂ t (cid:16) εE + χ π E (cid:17) = ∆ E (2)Let the coordinate system be set so that the functions ε − ≥ χ ≥ z − directionand the system has the translational symmetry along the y − direction. In this case, ε , χ and E are functions of z and x . In the asymptotic regions | z | → ∞ , Eq. (2) becomes a linear wave equation. So, the scattering problemcan be considered for a plane wave of the frequency ω that propagates from the asymptotic region z → −∞ to theregion z → ∞ . Furthermore, it is assumed that the functions ε and χ are piecewise constant, i.e., ε = ε c = const and χ = χ c = const in regions occupied by the scattering system. A conventional treatment of the problem is based onthe assumption that the solution E is analytic in χ c and, therefore, can be represented as a power series expansion, E = 2Re (cid:8) E e − iωt + χ c E e − iωt + χ c (cid:0) E , e − iωt + E , e − iωt (cid:1) + . . . (cid:9) (3)where E is the amplitude of the fundamental harmonics in the zero order of χ c , E is the amplitude of the secondharmonics in the first order of χ c , and so on. This assumption is not true if the system has bound states in theradiation continuum. Indeed, a general solution has the form E = E L + E NL , where E L is the solution when χ c = 0and E NL is the correction due to nonlinear effects. Let χ be written as χ = χ c η , where η is the indicator functionof the region occupied by the scattering system, i.e., its value is 1 in that region and 0 elsewhere. Then, if b G isthe Green’s function of the operator εc ∂ t − ∆ with appropriate (scattering) boundary conditions, the function E NL satisfies the integral equation E NL = − χ c πc b G (cid:2) η∂ t ( E L + E NL ) (cid:3) The power series expansion (3) can be obtained by the method of successive approximations for this integral equation,provided the series is proved to converge. According to scattering theory [18, 19], the Fourier transform of b G ismeromorphic in k = ω c . As is clarified shortly (see discussion of Eq.(11)), its real poles correspond to boundstates in the radiation continuum. Hence, in the presence of a real pole k = k b , the kernel of b G is not summableand, therefore, the successive approximations produce a diverging series. This implies a non-analytic behavior of thesolution in χ c . Thus, when a bound state in the radiation continuum is present, the conventional perturbative approachbecomes inapplicable. Here, a non-perturbative approach is developed to obtain the solution to the scattering problemthat is valid in any small neighborhood of a real pole of the Fourier transform of b G .Suppose that the incident radiation is a plane wave E in ( r , t ) = 2 cos( k · r − ωt ) , k = k x e + k z e , ck = ω, where e i , i = 1 , ,
3, denote unit vectors along the x , y , and z coordinate axes, respectively. A general solution to Eq.(2) should then be of the form, E ( r , t ) = ∞ X l = −∞ E l ( r ) e − ilωt (4)where E ≡
0, and for all l , E − l = E l is the complex conjugate of E l (as E is real). Therefore it is sufficient todetermine only E l , l ≥
1. Next, it is assumed that the scattering structure is periodic in the x − direction (e.g., agrating). The units of length are chosen so that the period is one. Then the amplitudes E l satisfy Bloch’s periodicitycondition E l ( r + e ) = e ilk x E l ( r ) (5)This condition follows from the requirement that the solution E satisfies the same periodicity condition as the incidentwave E in : E in ( r + e , t ) = E in (cid:16) r , t − k x ω (cid:17) By Eq.(2), the amplitudes of the different harmonics satisfy the equations,∆ E l + l k εE l = − νl k ( ε − X p E p E l − p , ν = χ c π ( ε c − ν is often used in lieu of χ c . Since ν ∼ χ c , it is a small parameter in the system.The scattering theory requires that for l = ±
1, the partial waves E l e − ilωt be outgoing in the spatial infinity( | z | → ∞ ). The fundamental waves E ± e ∓ iωt are a superposition of an incident plane wave e ± i ( k · r − ωt ) and ascattered wave which is outgoing at the spatial infinity. In all, the above boundary conditions lead to a system ofLippmann-Schwinger integral equations for the amplitudes E l : ( E = b H( k )[ E + ν P p E p E − p ] + e i k · r E l = b H(( lk ) )[ E l + ν P p E p E l − p ] , l ≥ E − l = E l for l ≤ −
1, where b H( q ) is the integral operator defined by the relation b H( q )[ ψ ]( r ) = q π Z ( ε ( r ) − G q ( r | r ) ψ ( r ) d r (7)in which G q ( r | r ) is the Green’s function of the Poisson operator, ( q + ∆) G q ( r | r ) = − πδ ( r − r ), with the outgoingwave boundary conditions. For two spatial dimensions, as in the case considered here r = ( x, z ) and r = ( x , z ), theGreen’s function is known [20] to be G q ( r | r ) = iπH ( q | r − r | ) where H is the zero order Hankel function of thefirst kind.When ν = 0, the amplitudes of all higher harmonics ( l ≥
2) vanish. Therefore it is natural to assume that | E | ≫ | E | ≫ | E | ≫ · · · for a small ν . Note that this does not generally imply that the solution, as a function of ν , is analytic at ν = 0. Under this assumption, the solution to the system (6) can be approximated by keeping onlythe leading terms in each of the series involved. In particular, the first equation in (6) is reduced to E ≈ e i k · r + b H( k )[ E ] + 2 ν b H( k ) (cid:2) E E (cid:3) (8)while the second equation becomes E ≈ b H((2 k ) )[ E ] + ν b H((2 k ) ) (cid:2) E (cid:3) (9)It then follows that a first order approximation to the solution of the nonlinear wave equation (2) may be found bysolving the system formed by the equations (8) and (9). To facilitate the subsequent analysis, the system is rewrittenas ( [1 − b H( k )][ E ] = e i k · r + 2 ν b H( k ) (cid:2) E E (cid:3) [1 − b H((2 k ) )][ E ] = ν b H((2 k ) ) (cid:2) E (cid:3) (10)Solving the first of Eqs.(10) involves inverting the operator 1 − b H( k ), and therefore necessitates a study of thepoles of the resolvent [1 − b H( k )] − . Such poles are eigenvalues in the generalized eigenvalue problem, b H( k )[ E ] = E (11)for fixed k x . The corresponding eigenfunctions E = E s are referred to as Siegert states . In contrast to Siegert statesin quantum scattering theory [18], electromagnetic Siegert states satisfy the generalized eigenvalue problem (11) inwhich the operator is a nonlinear function of the spectral parameter k . Their general properties are studied in [15].Eigenvalues have the form k = k r − i Γ. If k r > k x , then, according to scattering theory, such a pole is a resonancepole. In the case of the linear wave equation ( χ c = ν = 0), the scattered flux peaks at k = k r indicating the resonanceposition, whereas the imaginary part of the pole Γ defines the corresponding resonance width (or a spectral width ofthe scattered flux peak; a small Γ corresponds to a narrow peak). If Γ = 0, the corresponding Siegert state is a boundstate . This is a localized (square integrable) solution of Eq.(11). Suppose that the scattering system has a physicalparameter h such that a pole k = k r ( h ) − i Γ( h ) of the resolvent [1 − b H ( k )] − depends continuously on h and thatthere is a particular value h = h b at which the pole becomes real, i.e., Γ( h b ) = 0. If k b = k r ( h b ) > k x , then thecorresponding bound state lies in the radiation continuum. Note that Eq.(11) may have solutions for real k < k x .These are bound states below the radiation continuum. Such states are not relevant for the present study and,henceforth, bound states are understood as bound states in the radiation continuum. As noted in the introduction,two periodic planar scattering structures separated by a distance 2 h have bound states in the radiation continuum.Suppose that for fixed k x , the set k b consists of isolated points (which is generally true) and the points (2 k b ) donot belong to it (which is fulfilled in a concrete example studied in the next section [6]). Consequently, the operator(1 − b H((2 k ) ) − is regular in a neighborhood of k b and, for k close to but not equal to k b , the operator 1 − b H( k ) isinvertible. It then follows that E satisfies the nonlinear integral equation, E = (cid:16) − b H( k ) (cid:17) − (cid:20) e i k · r + 2 ν b H( k ) (cid:20) E (cid:16) − b H((2 k ) ) (cid:17) − hb H((2 k ) ) (cid:2) E (cid:3)i(cid:21)(cid:21) (12)where, in accord with the notation introduced in (7), the function on which an operator acts is placed in the squarebrackets following the operator. The operator ν (1 − b H( k )) − that determines the “nonlinear” part of Eq. (12) isnot bounded as k → k b , no matter how small ν ∼ χ c is. This precludes the use of a power series representationof the solution in ν . To find the solution of Eq.(12) when ν is small, its property under parity transformations isestablished first.Suppose that the scattering system is such that the operator b H and the parity operator b P defined by b P[ E ]( x, z ) = E ( x, − z ) commute. This implies that the Siegert states have a specific parity: b P[ E s ] = pE s where p = ±
1. Considerthen the ratio µ ( x, z ) = E ( x, − z ) E ( x, + z ) = b P[ E ] E (13)It will be proved that µ ( x, z ) → p in the limit ( h, k ) → ( h b , k b ) along a certain curve. Indeed, it follows from themeromorphic expansion of [1 − b H( k )] − that near a pole k r ( h ) − i Γ( h ), E = iC ( h ) k − k r ( h ) + i Γ( h ) E s + O (1) (14)where C ( h ) is some constant depending on h , and E s is an appropriately normalized Siegert state [15]. Consider thenthe curve of resonances C : k = k r ( h ) in the ( h, k )-plane. Along C , E ( x, z ) = C ( h )Γ( h ) E s ( x, z ) + O (1) (15)Now, as h → h b , the width Γ( h ) goes to 0, and the Siegert state E s becomes a bound state E b in the radiationcontinuum. Equation (15) shows that if C ( h ) does not go to zero faster than Γ( h ) as h → h b , i.e., the pole still givesthe leading contribution to E in this limit, then µ ( x, z ) → E b ( x, − z ) E b ( x, + z ) = b P[ E b ] E b = p = ± E b is even or odd in z . For the linear wave equation ( ν = 0), the constant C ( h ) is shown to be proportional to p Γ( h ) [15]. Therefore, for a small ν , the assumption that C ( h ) does not go tozero faster than Γ( h ) as h → h b is justified.Based on the limit (16), the following (non-perturbative) approach is adopted to solve Eq.(12) near a bound state.First, the curve of resonances C in the ( h, k )-plane is found. Using the relation (13) in the right side of (12), the field E is expressed via the ratio µ with the pair ( h, k ) being on the curve C . Next, the principal part of the amplitude E relative to ∆ h = h − h b is evaluated near a critical point ( h b , k b ) on C by taking µ to its limit value (16). Thisapproach reveals a non-analytic dependence of the amplitude E on the small parameters ∆ h and χ c of the systemand allows to obtain E and E when a bound state is present in the radiation continuum. As the technicalities ofthe proposed non-perturbative approach depend heavily on peculiarities of the scattering system, the procedure isillustrated with a specific example. III. A PERIODIC DOUBLE ARRAY OF SUBWAVELENGTH CYLINDERS
The system considered is sketched in Fig. 1(a). It consists of an infinite double array of parallel, periodicallypositioned cylinders. The cylinders are made of a nonlinear dielectric material with a linear dielectric constant ε c > χ c ≪
1. The coordinate system is set so that the cylinders are parallel to the y-axis,the structure is periodic along the x-axis, and the z-axis is normal to the structure. The unit of length is taken to bethe array period, and the distance between the two arrays relative to the period is 2 h .The solution of the integral equation (12) is obtained for k near k b in the limit of subwavelength dielectric cylinders.The approximation is defined by a small parameter δ ( q ) = ( qR ) ε c − ≪ q on a single cylinder of radius R . For sufficientlysmall R , this approximation is justified. The integral kernel of b H( q ) is defined by (7) and has support on the region FIG. 1. Panel (a): Double array of dielectric cylinders. The unit of length is the array period. The axis of each cylinder isparallel to the y -axis, and is at a distance h from the x -axis.Panel (b): The scattering process for the normal incident radiation ( k x = 0). The scattered fundamental harmonic is symbolizedby a single headed arrow while the (generated) second harmonic radiation is symbolized by a double headed arrow. The incidentradiation wave length is such that only one diffraction channel is open for the fundamental harmonic while three diffractionchannels are open for the second harmonic. The flux measured through the faces L ± / : x = ± / ω r = ck r ) of scattering resonances as functions of thedistance between the arrays, k = k r ( h ). The dots on the curves indicate positions of bound states in the radiation continuum(i.e., the values of h at which a resonance turns into a bound state). The solid curve connects bound states symmetric relative tothe reflection z → − z . The dashed line connects the skew symmetric bound states. The curves are realized for R = 0 . , ε c = 2,and k x = 0 (normal incidence). occupied by cylinders. The condition (17) implies that the wavelength is much larger than the radius R , and thereforefield variations within each cylinder may be neglected, so that ψ ( x, z ) ≈ ψ ( n, ± h ) where ( n, ± h ) are the positions ofthe axes of the cylinders ( n is an integer). The integration in b H( q )[ ψ ] yields then an infinite sum over positions ofthe cylinders. By Bloch’s condition, ψ ( n, ± h ) = e inq x ψ (0 , ± h ), so that the function b H( q )[ ψ ]( x, z ) is fully determinedby the two values ψ (0 , ± h ). In particular, b H( q )[ ψ ](0 , ± h ) ≈ αψ (0 , ± h ) + βψ (0 , ∓ h ) (18)where the coefficients α and β are shown to be [6] α ( q, q x ) = 2 πiδ ( q ) ∞ X m = −∞ (cid:18) q z,m − πi ( | m | + 1) (cid:19) + iπ ln(2 πR ) ! β ( q, q x , h ) = 2 πiδ ( q ) ∞ X m = −∞ e ihq z,m q z,m where q z,m = p q − ( q x + 2 πm ) with the convention that if q < ( q x + 2 πm ) , then q z,m = i p ( q x + 2 πm ) − q .To obtain the energy flux scattered by the structure, the action of the operator b H( q ) on ψ must be determined inthe asymptotic region | z | → ∞ . It is found that for | z | > h + R , b H( q )[ ψ ]( x, z ) ≈ πiδ ( q ) ∞ X m = −∞ (cid:16) ψ (0 , h ) e i | z − h | q z,m + ψ (0 , − h ) e i | z + h | q z,m (cid:17) e ix ( q x +2 πm ) q z,m (20) IV. AMPLITUDES OF THE FUNDAMENTAL AND SECOND HARMONICS
Now that the action of the operator b H( q ) has been established in (18) and (20), the amplitudes E and E of thefundamental and second harmonics can be determined by solving the system (10). As noted earlier, this will be donealong a curve C in the h, k -plane defined by k = k r ( h ) where k r ( h ) is the real part of a pole of [1 − b H( k )] − , orequivalently, when the incident radiation has the resonant wave number k = k r ( h ). To find the curve, the eigenvalueproblem (11) is solved in in the approximation (18):[1 − H ] (cid:18) E b + E b − (cid:19) = (cid:18) (cid:19) , H = (cid:18) α ββ α (cid:19) (21)where E b ± = E b (0 , ± h ) and the functions α = α ( k, k x ) and β = β ( k, k x ) have been defined in the previous section.In particular, bound states occur at the points ( h b , k b ) at which the determinant det(1 − H ) vanishes,det (cid:18) − α − β − β − α (cid:19) = (1 − α − β )(1 − α + β ) = 0It follows from Eq.(21) that the bound states for which 1 − α − β = 0 are even in z because E b + = E b − in thiscase. Similarly, the bound states for which 1 − α + β = 0 are odd in z . More generally, the poles of the resolvent[1 − b H( k )] − are complex zeros of det(1 − H ). They are found by the conventional scattering theory formalism.Specifically, the resonant wave numbers k = k r ( h ) are obtained by solving the equation Re { − α ± β } = 0 for thespectral parameter k . According to the convention adopted in the representation (14), the corresponding resonancewidths are defined by Γ( h ) = − Im { − α ± β } ∂ k Re { − α ± β } (cid:12)(cid:12)(cid:12)(cid:12) k = k r ( h ) where ∂ k denotes the derivative with respect to k . This definition of the width corresponds to the linearizationof Re { − α ± β } near k = k r ( h ) as a function of k in the pole factor [1 − α ± β ] − . The curves of resonances k = k r ( h ) > k x come in pairs. There is a curve connecting the symmetric bound states in the h, k -plane, and anothercurve that connects the odd ones.In what follows, only the curve connecting symmetric bound states will be considered. The other curve can betreated similarly. Panel (c) of Fig. 1 shows that the first symmetric bound state occurs when the distance 2 h isabout half the array period, while the skew-symmetric bound states emerge only at larger distances. This feature isexplained in detail in [6]. So, the solution obtained near the first symmetric bound state corresponds to the smallestpossible transverse dimension of the system (roughly a half of the wave length of the incident radiation). Thus, fromnow on the curve of resonances C refers to the curve in the h, k -plane defined by the equation Re { − α − β } = 0. Tosimplify the technicalities, it will be further assumed that only one diffraction channel is open for the fundamentalharmonics, i.e., k x < k < π − k x .Let k = k r ( h ) be the solution of Re { − α − β } = 0. By making use of the explicit form of the functions α and β for one open diffraction channel, one infers that along the curve k = k r ( h ),1 − α − β = i Im { − α − β } = − i πδ ( k b ) k z ϕ , ϕ = cos( hk z ) , k z = q k b − k x Bound states in the radiation continuum occur when the distance h satisfies the equation ϕ = 0, i.e., h p k r ( h ) − k x =( n − / π with n being a positive integer. Its solutions h = h b ( n ) define the corresponding values of the wave numbersof the bound states, k b ( n ) = k r ( h b ( n )). So, the sequence of pairs { ( h b ( n ) , k b ( n )) } ∞ n =1 indicates positions of the boundstates on the curve C . In the limit h → h b ( n ) along C , the function ϕ has the asymptotic behavior, ϕ = ( − n k z,b ∆ h + o (∆ h ) , k z,b = k z | h = h b ( n ) , ∆ h = h − h b The objective is to determine the dependence of the amplitudes E and E on the parameters ∆ h and χ c which areboth small.To this end, let E ± = E (0 , ± h ) be the values of the field E on the axes of the cylinders at (0 , ± h ), and E ± = E (0 , ± h ) be the values of the field E on the same cylinders. In the subwavelength approximation, thesevalues determine the scattered field because the latter is produced by the radiation of point dipoles induced by theincident wave on the scatterers and the strength of the dipoles is proportional to E ± for the fundamental harmonicand E ± for the second harmonic. Applying the rule (18) to evaluate the action of the operator b H( q ) in the system(10), the first equation of the latter becomes,[1 − H ] (cid:18) E E − (cid:19) = 2 ν H (cid:18) E E E − E − (cid:19) + (cid:18) e + ihk z e − ihk z (cid:19) (22a)Similarly, the second of Eqs.(10) yields,[1 − H ] (cid:18) E E − (cid:19) = ν H (cid:18) E E − (cid:19) , H = H (2 k, k x ) (22b)As stated above, the resolvent [1 − b H((2 k ) )] − is regular in a neighborhood of k b so that Eq. (22b) can be solvedfor E ± , which defines the latter as functions of E ± . The substitution of this solution into Eq.(22a) gives a systemof two nonlinear equations for the fields E ± . Adding these equations and replacing the field E − by its expression E − = µ (0 , h ) E in terms of the field ratio of Eq.(13) yields the following implicit relation between the field E and its amplitude | E | : E = − ϕ ν ζ | E | + ϕ ξ (23)where ϕ = cos( hk z ) and ν = χ c π ( ε c − are small and, in terms of the field ratio µ ≡ µ (0 , h ), the values of ζ and ξ read, ξ = i πδ ( k ) k z (1 + µ ) , ζ = (cid:18) i πδ ( k ) k z ϕ (cid:19) (cid:0) a + bµ + µ (cid:0) b + aµ (cid:1)(cid:1) (24)with a and b being defined by the relation, [1 − H ] − H = (cid:18) a bb a (cid:19) (25)In particular, ζ and ξ are continuous functions of µ and ϕ . In Appendix A it is shown that if ζ b and ξ b are the respectivelimits of ζ and ξ as h → h b along the curve of resonances C , then these limits are nonzero. It follows then that Eq.(23)for E is singular in both ν and ϕ when these parameters are small, i.e., in the limit ( ν, ϕ ) → (0 , X + 2 ϕ ν Re { ζξ } X + ϕ ν | ζξ | X − ϕ ν | ζ | = 0 , X = | E | (26)The solution to this cubic equation is obtained in Appendix C. It is proved there that Eq. (26) admits a uniquereal solution so that there is no ambiguity on the choice of E . In the vicinity of a point ( h b , k b ) along the resonancecurve C , the field E is found to behave as, | E | = | ∆ h | / χ c τ (∆ h, χ c ) (27)Recall that ∆ h = h − h b . An explicit form of the function τ (∆ h, χ c ) is given in Appendix C (see Eq. (C3)). It involvescombinations of the square and cube roots of functions in ∆ h and χ c and has the property that τ (∆ h, χ c ) → h, χ c ) → (0 ,
0) (in the sense of the two-dimensional limit). In the limit h → h b , the field ratio µ approaches 1 fora symmetric bound state as argued earlier. Therefore it follows from Eqs.(22b) that E ± ∼ νE because the matrix(25) exists at h = h b . Since ν ∼ χ c , relation (27) leads to the conclusion that E ± E = O ( | ∆ h | / )Thus, the approximation | E | ≫ | E | used to truncate the system (10) remains valid for h close to the critical value h b despite the non-analyticity of the amplitudes at (∆ h, χ c ) = (0 , V. FLUX ANALYSIS: THE CONVERSION EFFICIENCY
For the nonlinear system considered, even though Poynting’s theorem takes a slightly different form as comparedto linear Maxwell’s equations, the flux conservation for the time averaged Poynting vector holds. The scatteredenergy flux carried across a closed surface by each of the different harmonics adds up to the incident flux across thatsurface. The flux conservation theorem is stated in Appendix B. Consider a closed surface that consists of four faces, L ± = { ( x, z ) | − ≤ x ≤ , z → ±∞} and L ± / = { ( x, z ) | x = ± / } as depicted in Fig. 1(b). As argued in AppendixB, the scattered flux of each l th -harmonics across the union of the faces L ± / vanishes because of the Bloch condition(and so does the incident flux for any k x ). Therefore only the flux conservation across the union of the faces L ± hasto be analyzed. If σ l designates the ratio of the scattered flux carried by the l th -harmonics across the faces L ± to the0incident flux across the same faces, then P l ≥ σ l = 1. Thus, for l ≥ σ l defines the conversion ratio of fundamentalharmonics into the l th -harmonics.In the perturbation theory used here, only the ratios σ and σ may be evaluated. By laborious calculations it canbe shown that σ + σ ≤ σ as a function ofthe parameter h at a given value of the nonlinear susceptibility χ c .The ratio σ is defined in terms of the scattering amplitudes of the second harmonic, i.e., by the amplitude of E in the asymptotic region | z | → ∞ : E ( r ) → X m op,sh R shm e i r · k − m,sh , z → −∞ X m op,sh T shm e i r · k + m,sh , z → + ∞ (28)where k ± m,sh = (2 k x + 2 πm ) e ± k shz,m e is the wave vector of the second harmonic in the m th open diffraction channel.Recall that the m th channel is open provided (2 k ) > (2 k x + 2 πm ) and in this case k shz,m = p (2 k ) − (2 k x + 2 πm ) ,while if the channel is closed, then k shz,m = i p (2 k x + 2 πm ) − (2 k ) . In the asymptotic region | z | → ∞ , the field inclosed channels decays exponentially and, hence, the energy flux can only be carried in open channels to the spatialinfinity. The summation in Eqs.(28) is taken only over those values of m for which the corresponding diffractionchannel is open for the second harmonic, which is indicated by the superscript “ op, sh ” in the summation index m op,sh . Note that there is more than one open diffraction channel for the second harmonic even though only onediffraction channel is open for the fundamental one. For instance, if the x − component of the wave vector k , i.e., k x , is less than π , there are 3 open diffraction channels for the second harmonic, the channels m = 0 , m = −
1, and m = 1. These three directions of the wave vector of the second harmonic propagating in each of the asymptotic regions z → ±∞ are depicted in Fig. 1(b) by double-arrow rays. Thus, in terms of the scattering amplitudes introduced inEqs. (28), the ratio of the second harmonic flux across L ± to the incident flux is σ = 12 k z X m op,sh k shz,m (cid:0) | R shm | + | T shm | (cid:1) The scattering amplitudes R shm and T shm are inferred from Eq.(9) in which the rule (20) is applied to calculate theaction of the operator b H((2 k ) ) in the far-field regions | z | → ∞ : R shm = 2 πiδ (2 k ) k shz,m h ( E + νE ) e ihk shz,m + ( E − + νE − ) e − ihk shz,m i T shm = 2 πiδ (2 k ) k shz,m h ( E + νE ) e − ihk shz,m + ( E − + νE − ) e ihk shz,m i Since for ν = 0, the amplitudes E ± remain finite as h → h b along C , and since µ → σ in a vicinity of a bound state along the curve C is obtained by setting E − = E in Eq. (22b), solving thelatter for E ± , and substituting the solution into the expression for σ . The result reads σ = C b ν | E | (29)where C b is a constant obtained by taking all nonsingular factors in the expression of σ to their limit as h → h b ,which gives C b = " (16 πδ ( k )) k z | a + b | X m op,sh cos ( hk shz,m ) k shz,m ( h,k )=( h b ,k b ) for a and b defined in Eq.(25). Using the identity | E | = | E | | E | , and substituting Eq.(23) into one of thefactors | E | , the conversion ratio σ is expressed as a function of a single real variable, σ ( u ) = C ′ b u | u + ζ b ξ b | , u = (cid:18) ν | E | ϕ (cid:19) (30)where C ′ b = C b | ζ b | is a constant, and ξ and ζ in Eq.(24) have been taken at their limits as h → h b to obtain theprincipal part of σ . The function u σ ( u ) on [0 , ∞ ) is found to attain its absolute maximum at u = | ξ b ζ b | . This1condition determines the distance 2 h between the arrays at which the conversion rate is maximal for given parameters R, ε c and χ c of the system. Indeed, since u ∼ | E | should also satisfy the cubic equation (26), the substitution of u = | ζ b ξ b | into the latter yields the condition ν ϕ = 2 | ξ b | ( | ξ b ζ b | + Re { ξ b ζ b } ) (31a)In particular, in the leading order in δ ( k ), the optimal distance 2 h between the two arrays is given by the formula,( h − h b ) = χ c π k z,b ( k b R ) ( ε c − (31b)where as previously, k z,b = p k b − k x .The maximum value σ , max of the conversion ratio σ is the sought-for conversion efficiency . An interesting featureto note is that σ , max = σ ( | ζ b ξ b | ) is independent of the nonlinear susceptibility χ c because the constants C ′ b , ξ b ,and ζ b are fully determined by the position of the bound state ( k b , h b ). In other words, if the distance between thearrays is chosen to satisfy the condition (31a), the conversion efficiency σ , max is the same for a wide range of valuesof the nonlinear susceptibility χ c . This conclusion follows from two assumptions made in the analysis. First, thesubwavelength approximation should be valid for both the fundamental and second harmonics, i.e., the radius ofcylinders should be small enough. Second, the values of h − h b and ν (or χ c ) must be such that the analysis of theexistence and uniqueness of | E | given in Appendix C holds, that is, Eq. (26) should have a unique real solutionunder the condition (31a). The geometrical and physical parameters of the studied system can always be chosen tojustify these two assumptions as illustrated in Fig. 2. FIG. 2. Panel (a): The conversion efficiency is plotted against the cylinder radius R for the critical points ( h b ( n ) , k b ( n )) , n =1 , ,
3, and ε c = 1 . k x = 0 (the normal incidence). The dashed parts of the curves indicate the regions where δ ( k ) > . δ (2 k ) >
1, i.e., the subwavelength approximation becomes inapplicable for the second harmonic.Panel (b): The conversion efficiency is plotted against k x for the critical points ( h b ( n ) , k b ( n )) , n = 1 , ,
3. The curves arerealized for R = 0 .
15 and ε c = 1 . h b (1) ≈ . ν, h ) − plane is defined by the condition τ + + τ − > | E | exists and uniqueas explained in Appendix C. The plot is realized for R = 0 . , ε c = 2 and k x = 0. The parabola-like curve is an actual boundaryof the shadowed region; the top horizontal line represents no restriction. So there is a wide range of the physical and geometricalparameters within the shadowed region which satisfy (31a). The regions of validity for other bound states looks similar. Panels (a) and (b) of Fig. 2 show the conversion efficiency σ , max for the first three symmetric bound states n = 1 , ,
3, as, respectively, a function of the cylinder radius R when k x = 0 and of k x when R = 0 .
15. For all curvespresented in the panels, ε c = 1 .
5. The values of σ , max are evaluated numerically by Eq. (30) where u = | ξ b ζ b | . Thesolid parts of the curves in Panel (a) correspond to the scattering phase δ ( k ) < .
25 with k = k b . Note that thewavelength at which the second harmonic generation is most efficient is the resonant wavelength defined by k = k r ( h )where h satisfies the condition (31a). For a small χ c , the scattering phase at the resonant wavelength can wellbe approximated as δ ( k ) ≈ δ ( k b ). The condition δ ( k b ) < .
25 ensures that the scattering phase for the secondharmonic satisfies the inequality δ (2 k b ) = 4 δ ( k b ) < δ ( k b ) > . h of the system studied here can be as low as a half of the wavelength, i.e., for an infraredincident radiation, 2 h is about a few hundred nanometers. Indeed, as one can see in Fig. 1(c), the first bound stateoccurs at h b (1) ≈ .
259 and k b ≈ π which corresponds to the wavelength λ b = 2 π/k b ≈ δ ( k ): σ , max = σ ( | ζ b ξ b | ) ≈ πδ ( k ) X m op,sh cos ( hk shz,m ) k shz,m (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ( h,k )=( h b ,k b ) (32)Suppose that only one diffraction channel is open for the incident radiation. Then m = 0 , ± σ , max denote the term m = 0, i.e. σ , max is the second harmonic flux inwhich the contribution of the channels with m = ± σ , max < σ , max . One infers from Eq.(32) that, σ , max ≈ πδ ( k b ) k z,b cos (2 h b k z,b )As the pair ( h b , k b ) at which a symmetric bound state is formed satisfies the equation cos( h b k z,b ) = 0, it follows thatcos(2 h b k z,b ) = −
1. Hence, σ , max ≈ k b k z,b πR ( ε c − k b at which the bound states occur lie just below the diffraction threshold 2 π − k x , i.e., k b / π − k x (see Fig. 1(c) and [6] for details). So that in the case of normal incidence ( k x = 0), the above estimate becomes, σ , max ≈ π ( πR )( ε c − δ (2 π ) ≪
1. If, for instance, R = 0 .
15 and ε c = 2, then, σ , max ≈ δ (2 π ) ≈ . Appendices
Appendix A: Estimation of ζ and ξ Here the limit values ζ b and ξ b of the functions ζ and ξ defined in Eq. (24) are estimated as h → h b along theresonance curve C . For ξ b this is immediate. Indeed, in the aforementioned limit, the field ratio µ → ξ b = i πδ ( k ) k z (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) k = k b For ζ b , the estimate follows from that wave numbers k b at which bound states exist are close to the diffraction threshold2 π − k x when only one diffraction channel is open for the fundamental harmonic, i.e., k x < k < π − k x [6]. Indeed,in the first order of δ ( k ), k b ≈ π − k x − π δ (2 π − k x )2 π − k x (A1)This proximity of the wavenumbers k b to the diffraction threshold 2 π − k x allows one to determine the leadingterms in the coefficients a and b defined in Eq.(25), and hence ζ b . To proceed, the coefficients α = α (2 k, k x ) and β = β (2 k, k x , h ) of the matrix H are rewritten by separating explicitly the real and imaginary parts: α + β = ψ + + iS c , α − β = ψ − + iS s (A2a)3where S c and S s are defined by the relations, S c = 16 πδ ( k ) X m op,sh cos ( hk shz,m ) k shz,m , S s = 16 πδ ( k ) X m op,sh sin ( hk shz,m ) k shz,m (A2b)and the index m op,sh indicates that the summations are to be taken over all open diffraction channels for the secondharmonic. Using the estimate (A1), the functions ψ ± are found to obey the estimates, ψ + = 2 + O ( δ ( k b )) , ψ − = O ( δ ( k b ))These expressions are then used to estimate ζ − b = 2( a + b ). In the first order of δ ( k b ) one infers that ζ b ≈ − − πiδ ( k ) X m op,sh cos ( hk shz,m ) k shz,m (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ( h,k )=( h b ,k b ) (A3) Appendix B: Complements on the flux analysis: Flux conservation
For the nonlinear wave equation (2), the Poynting Theorem becomes,18 π ddt (cid:20)Z V (cid:16) εE + B + χ π E (cid:17) d r (cid:21) = − Z ∂V S · d n (B1)where V is a closed region, and ∂V is its boundary. The vector S = E × B is the Poynting vector (for simplicity, it isassumed that ∂V lies in the vacuum so that µ = ε = 1, and χ = 0 in a small neighborhood of ∂V ). In the case of amonochromatic incident wave, the flux measured is the time-average of S over a time interval T → ∞ . By averagingEq.(B1), it then follows that, Z ∂V h S i · d n = , h S i = 1 T Z T S ( t ) dt, T → ∞ This is the flux conservation. In terms of the different harmonics of Eq.(4), the time averaged Poynting vector becomes, h S i = − c πω Im ∞ X l =1 E l ∇ E − l l ! Of interest is the flux of the Poynting vector across the rectangle depicted in Fig. 1(b). By Bloch’s condition (5), thecontributions to the flux from the faces L ± / : x = ± cancel out so that the flux measured is through the verticalfaces L ± = { ( x, z ) |− ≤ x ≤ , z → ±∞} . Note that the vanishing of the flux across the faces L ± / is a consequenceof the fact that the incident wave is uniformly extended over the whole x − axis. For example, consider the normalincidence ( k x = 0) with one diffraction channel open for the incident radiation. Then the Poynting vector of thereflected and transmitted fundamental harmonic is normal to the structure and, hence, carries no flux across L ± / .The second harmonic ( l = 2) has three forward and three backward scattering channels open, m = 0 , ±
1, relative tothe z − axis. The wave with m = 0 propagates in the direction normal to the structure and does not contribute tothe flux across L ± / . Since the incident wave has an infinite front along the x − axis, so do the scattered waves with m = ±
1. The waves with m = 1 and m = − L ± / as the correspondingwave vectors have the same z − components and opposite x − components and, hence, the total flux vanishes. For afinite wave front (but much larger than the structure period), the second harmonic would carry the energy flux in allthe directions parallel to the corresponding wave vectors in each open diffraction channel.If σ l is as defined in Section V, then the flux conservation implies that P ∞ l =1 σ l = 1. Therefore, in the perturbationtheory used, i.e., when the system (6) is truncated to Eqs.(10), the inequality σ + σ ≤ σ is given in Section V. If only one diffraction channel is open for the fundamental harmonic,then the ratio σ of the scattered and incident fluxes of the fundamental harmonic reads, σ = | T | + | R | T and R are the transmission and reflection coefficients which are obtained from the far-field amplitude of E as, E → ( e i r · k + R e i r · k − , z → −∞ (1 + T ) e i r · k , z → + ∞ where k = k x e + k z e is the incident wave vector and k − = k x e − k z e is the wave vector of the reflected fundamentalharmonic. It then follows from Eqs.(8) and (20) that R = i πδ ( k ) k z h ( E + 2 νE E ) e ihk z + ( E − + 2 νE − E − ) e − ihk z i T = i πδ ( k ) k z h ( E + 2 νE E ) e − ihk z + ( E − + 2 νE − E − ) e ihk z i In the vicinity of a critical point ( h b , k b ), the coefficients R and T obey the estimate, R ≈ T ≈ i πδ ( k b ) k z,b ϕE (cid:18) ν | E | ζ b (cid:19) After some algebraic manipulations, it is found that, σ + σ = 1 + 8 πδ ( k b ) k z,b ν | E | (cid:20) A b + 4 πδ ( k b ) k z,b ϕ (cid:18) Re (cid:26) ζ b (cid:27) + ν | E | | ζ b | (cid:19)(cid:21) where A b is the constant defined as, A b = (cid:20) S c | a + b + 1 | − Im (cid:26) ζ (cid:27)(cid:21) ( h,k )=( h b ,k b ) and S c = Im { α + β } is introduced in Eqs.(A2). Expressing a and b defined by (25) via the coefficients α and β of the symmetric matrix H , one also obtains S c = Im (cid:26) a + b a + b (cid:27) Since at the point ( h b , k b ) the value of ζ is ζ b = (2( a + b )) − | ( h,k )=( h b ,k b ) , it follows that, A b = (cid:20) (cid:26) a + b a + b (cid:27) | a + b + 1 | − { a + b } (cid:21) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ( h,k )=( h b ,k b ) For general complex numbers a and b , the expression in square brackets is always zero. Therefore A b = 0, and, σ + σ = 1 + 2 (cid:18) πδ ( k b ) k z,b ϕν | E | (cid:19) (cid:18) Re (cid:26) ζ b (cid:27) + ν | E | | ζ b | (cid:19) (B2)By Eq. (A3), Re { ζ − b } ≈ −
4. In Appendix C it is proved that ν | E | = O ( ϕ / ). Consequently, near the criticalpoint ( h b , k b ), the right hand summand in Eq.(B2) is negative so that σ + σ ≤ Appendix C: Complements on the amplitude E The amplitude of the field E is a root of the cubic polynomial in Eq.(26) which can be solved by Cardano’smethod. Put Y = X + (cid:0) ϕν (cid:1) Re { ζξ } . For the new variable Y , Eq. (26) assumes the standard form, Y + pY + q = 0 (C1)where, p = ϕ ν (cid:0) | ζξ | − { ( ζξ ) } (cid:1) , q = 2 ϕ ν Re { ζξ } (cid:0) { ( ζξ ) } − | ζξ | (cid:1) − ϕ ν | ζ | E is uniquely defined by the system (6), it is therefore expected that the cubic in Eq.(C1) shouldhave a unique real solution in order for the theory to be consistent. The latter holds if and only if the discriminant D = 427 p + q is nonnegative. To prove that D ≥
0, note first that | ζξ | − { ζ ξ } > h b , k b ).This follows from the estimates established in Appendix A. Indeed, in the first order of δ ( k b ), | ζ b ξ b | − { ( ζ b ξ b ) } ≈ π k z δ ( k ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) k = k b > ρ = 427 (cid:20) { ζξ } (cid:18) { ( ζξ ) } − | ζξ | (cid:19) + i ( | ζξ | − { ( ζξ ) } ) (cid:21) The positivity condition (C2) ensures that the coefficient of the complex number i in the expression of ρ is indeedreal. After some algebraic manipulations, it can be shown that, D = ϕ ν (cid:12)(cid:12)(cid:12) | ζ | ν − ϕ ρ (cid:12)(cid:12)(cid:12) Thus D ≥ Y to Eq.(C1) is then, Y = s − q + √ D s − q − √ D | E | = | ϕ | ν p τ + + τ − , τ ± = s (cid:18) ν | ζ | − ϕ Re { ρ } ± (cid:12)(cid:12)(cid:12) | ζ | ν − ϕ ρ (cid:12)(cid:12)(cid:12)(cid:19) − ϕ { ζξ } (C3)provided τ + + τ − ≥
0. The latter condition imposes a limit on the validity of the perturbation theory developed in thepresent study, i.e., the reduction of the system (6) to (10) is justified if τ + + τ − ≥
0. This is to be expected becauseof the lack of analyticity in χ c of the solution to the nonlinear wave equation (2) that can only occur at the criticalpoints ( h b , k b ) at which bound states in the radiation continuum exist. As one gets away from these critical pointsin the ( h, k )-plane, the solution to the nonlinear wave equation becomes analytic in χ c , meaning that all the termsthat were neglected in finding the principal parts of the amplitudes must now also be taken into account to find asolution befitting the series of Eq.(3). The shadowed region depicted in Fig. 2(b) shows the region of the ( h, k )-planein which the condition τ + + τ − ≥ h = h − h b and R , and thephysical parameters, ε c and χ c >
0, which satisfy the conditions (31a) and τ + + τ − ≥ [1] F.J. Garc´ıa de Abajo, Rev. Mod. Phys. 79, 1267 (2007).[2] S. V. Shabanov, Int. J. Mod. Phys. B 23, 5191 (2009).[3] D. C. Marinica, A. G. Borisov and S. V. Shabanov, Phys. Rev. B 76, 085311 (2007).[4] A. Nahata, R. A. Linke, T. Ishi, K. Ohashi, Opt. Lett. 28, 423 (2003).[5] D. C. Marinica, A. G. Borisov and S. V. Shabanov, Phys. Rev. Lett. 100, 183902 (2008).[6] F. R. Ndangali and S. V. Shabanov, J. Math. Phys. 51, 102901 (2010).[7] J. von Neumann and E. Wigner, Phys. Z 30, 465 (1929).[8] F. Capasso, C. Sirtori, J. Faist, D. L. Sivco, S.-N. G. Chu and A. Y. Cho, Nature (London) 358, 565 (1992).[9] H. Friedrich and D. Wintgen, Phys. Rev. A 32, 3231 (1985).[10] J. Okolowicz, M. Ploszajczak and I. Rotter, Phys. Rept. 374, 271 (2003).[11] A. Z. Devdariani, V. N. Ostrovsky and Yu. N. Sebyakin, Sov. Phys. JETP 44, 477 (1976).[12] E.N. Bulgakov and A.F. Sadreev, Phys. Rev. B 78, 075105 (2008).[13] E.N. Bulgakov and A.F. Sadreev, Phys. Rev. B 80, 115308 (2009).6