Leaky waves in planar dielectric waveguide

1 Department of Applied Probability and Informatics Peoples’ Friendship University of Russia (RUDN University) 6, Miklukho-Maklaya St., Moscow 117198, Russian Federation 2 A.M. Prokhorov General Physics Institute Russian Academy of Sciences, Moscow 119991, Russian Federation 3 Bogoliubov Laboratory of Theoretical Physics Joint Institute for Nuclear Research 6, Joliot-Curie St., Dubna, Moscow region 141980, Russian Federation


Introduction
In the books by Marcuse [1], [2], Adams [3], Snyder and Love [4], Tamir [5], and other authors the terms "leaky rays" and "leaky modes" appear when discussing the propagation of polarized light in fiber optical waveguides with the refractive index of the core smaller than that of the cladding, and in planar waveguides with plates of material optically denser than the waveguide layer itself. In this case, Marcuse writes that the outflow of light from such a waveguide is akin to tunneling through a potential barrier in quantum mechanics. The "leaky light", in contrast to the "emitted light", propagates for quite a long time along the axis of the optical fiber. Similarly, in a planar waveguide, the resulting electromagnetic radiation propagates for some time at some distance along the waveguide, in contrast to the emitted light.
At the initial stage of the study of "leaky" modes, T. Tamir et al. [6]- [12], A.W. Snyder et al. [13]- [18], as well as other research teams [19]- [25], investigated the dispersion equations of optical waveguides written in terms of transition matrices from the point of view of choosing one (two) roots of an analytical function. The studies were executed using the theory of residues in Cauchy integrals.
In the papers by V. Shevchenko [24], [25] the behavior of guided modes during the transition of their wave numbers beyond critical values was analyzed, their transformation into leaky modes was shown, and the choice of quadrants to which the wave numbers should shift when passing through critical values was justified. At frequencies below critical, the reflection from the waveguide walls ceases to be complete, so that the waveguide modes continue to propagate experiencing incomplete internal reflection, because of which some radiation from the waveguide occurs. Such (improper) waveguide waves with radiative damping are called leaky waves.
Open waveguides as radiating systems were first investigated by Hansen [19], who proposed an antenna structure implemented using leaky waves. However, there was no understanding of the physical mechanism of the resulting waves. After all, the leaky waveguide mode is characterized by a complex longitudinal wave number with constant attenuation due to radiation losses.
In this case, the longitudinal attenuation leads to an exponential increase in the wave amplitude in the transverse direction. This fact violates the usual radiation condition for guided modes, described by the solutions of self-adjoint problems for the Helmholtz equation. The behavior of the resulting waves that seems non-physical was clarified by Marcuvitz [23] and Oliner [6]- [9].
In the papers by Oliner et al. [6]- [12], a detailed study of the complex roots of the dispersion equation that do not correspond to the guided modes is given. The study begins with the assumption that exponentially damped waves correspond to such roots, the experimental observation of which was earlier reported in Refs. [20]- [22]. First, using ray technique, and then with the help of mode analysis, the authors analyzed the wave solutions corresponding to four different roots of the fourth-power dispersion equation. Two of these roots correspond to solutions that exponentially decrease in the direction of propagation and are located symmetrically with respect to the axes of coordinates and the origin of coordinates. The rest two roots are rejected. Many publications of that time have been devoted to the analysis of the relative position of the variety of roots [6]- [18].
In the first decades of the study of leaky waves, the method of steepest descent was used as the most common method for their numerical search. In this case, the trajectory of the fastest descent comes close to the leakage poles, they begin to make a significant contribution or even dominate in the general directional pattern of an open waveguide. The field distribution of the resulting leaky wave increases in the transverse direction. However, the field amplitude remains finite in a wedge-like region of space that allows leakage.
As shown by Marcuvitz [23], these complex poles can correspond to leaky modes. Although they do not make a direct contribution to the correct spectral solution and can therefore be characterized by non-physical growth towards infinity, they can nevertheless accurately describe the radiation field in limited spatial domains. In Ref. [26], e.g., it is noted that in most publications on leaky modes there are no plots of fields of various types for leaky modes calculated numerically (see, e.g., [4], [5], [27]- [35]). In this case, the authors of some publications (e.g., [34]) propose to replace the leaky modes with radiative ones in limited domains. Our studies have shown that this can lead, firstly, to a large error in the calculation of losses, and secondly, to an inaccurate calculation of the field profiles of leaky modes at distances exceeding several wavelengths (⩾ 2) of the electromagnetic radiation used. The replacement of one wave with another sometimes used requires serious analysis in each specific case. As a consequence, there is an urgent need to develop new algorithms for calculating the fields of both radiative and leaky modes, surpassing the standard methods, e.g., the FDTD method, in count rate and not inferior to them in accuracy.
In quantum physics, such solutions of the stationary Schrödinger equation are called Gamow resonances [36], [37] or Siegert quasi-states [38]. In recent decades, some researchers (see [39], [40]) solve boundary-value problems for the Helmholtz equations with the asymptotic conditions of Siegert leaky waves, obtaining numerical results interpreted by them as leaky waves. We propose to obtain (using a numerical method) the solutions of boundary problems for wave equations with asymptotic conditions of Siegert leaky waves. The numerical solutions obtained using this approach coincide with the solutions of Refs. [27]- [29], [33], [35], [41], but additionally allow description of the phase fronts of leaky waves and "angular outflow cones".
In our opinion, a more rigorous justification of the model of leaky waves of open waveguide systems can be obtained by starting calculations not from the Helmholtz equation, as is traditionally done, but from the wave equation preceding the Helmholtz equation, and most importantly, more adequately reflecting the wave nature of leaky modes.

Statement of the problem of modeling leaky modes of symmetric waveguides
Consider ( Figure 1) a symmetric three-layer planar waveguide consisting of a dielectric film having the height ℎ with real refractive index , surrounded by a cladding layer with real refractive index < . The propagation of radiation in such structures is described by the Maxwell equations, material equations [3]- [5], and boundary conditions that distinguish the class of solutions interesting for the researcher -in the present case, the leaky modes [6]- [12], [26], [42]. The generally accepted model of the electromagnetic field in a planar (infinitely extended along the -axis) are fields that are independent of the variable . In this case, Maxwell's equations are considerably simplified, since / = / ≡ 0 for any = , , , and they are divided into two independent subsystems -the subsystem for the so-called TE-modes and for the TM-modes. The subsystem for the TE-modes can be represented as a single wave equation for the master component with the boundary conditions and the initial conditions where ( ) = √ 2 − 2 , and is the electrodynamic constant, ( ) is the variable refractive index of the considered three-layer waveguide, defined below. The subsystem for the TE-mode also includes two equations for the connection of components , with the master component . As a model of leaky modes propagation, we will consider Eq. (1), i.e., in other words, we will consider the propagation of leaky modes in terms of a wave process. As asymptotic boundary conditions, we will consider the conditions of leaky modes corresponding to the Gamow-Siegert model [36]- [40].
In the case under consideration, the function describing the refractive index depends only on , which makes it possible to separate the variables in Eq. (1).
As a result, we obtain solutions corresponding to leaky modes propagating in the positive direction of the -axis: where is the frequency, and are the eigenvalues of the non-self-adjoint Sturm-Liouville problem with boundary conditions that extract the leaky modes [3], [34], [41]: The eigenfunctions of the problem (7) are defined as general solution of the ordinary differential equation subject to the boundary conditions of this problem, that is, they have the form and the constants , , , are determined from the field joining conditions at the boundaries of the waveguide layer = 0 and = ℎ, which with Eq. (8) taken into account constitute a homogeneous system of linear algebraic equations: The homogeneous system of Eqs. (9) has a nontrivial solution if and only if the determinant of the matrix of the system (9) is zero [43]. The equality to zero of the determinant of the matrix of the system (9) can be achieved for some values of the spectral parameter , which, in turn, determine the eigenvalues of the problem (7).
In each subdomain > ℎ, 0 < < ℎ, and < 0 the solution of the wave equation corresponding to the leaky modes is representable as a wave with a complex wave vector. In the case of a field corresponding to leaky modes running in the positive direction of the -axis for > ℎ and > ℎ due to the symmetry of the waveguide, the wave vector is determined as and in the waveguide layer (0 < < ℎ) there are two waves with wave vectors the modules of the wave vectors being equal to the corresponding wave numbers: | ⃗ ± | = 0 and | ⃗ ± | = 0 . We formulate the problem of finding solutions corresponding to the leaky modes as an eigenvalue problem for a differential operator with non-self-adjoint boundary conditions (7), which we will further solve numerically.

Description of the algorithm for numerical solution of the leaky mode problem
The spectral problem for a differential operator with non-self-adjoint boundary conditions (7) is formulated numerically as a problem of approximate determination of complex solutions of the equation where ( ) denotes the coefficient matrix of Eqs. (9) [42]. In the case of a problem similar to (7), but with self-adjoint boundary conditions, any classical method of finding the real roots of the equation can be applied (see, e.g., [44]- [46]). Our problem (7) is not self-adjoint, therefore, the eigenvalues of this problem are generally complex and the standard methods for root search can no longer be applied.
The problem (12) can be reformulated as a problem of finding the minimum of a function of two variables as follows. The desired quantity = ′ + ″ is a complex number. Any solutioñ=̃′ +̃″ of Eq. (12) will be also a local minimum of the non-negative function The eigenvalues corresponding to the leaky modes are localized in the first quadrant of the complex plane Re ( ) > 0, Im ( ) > 0; moreover 0 < Re ( ) < (see Refs. [34], [47], as well as our papers [42]).
To find all the local minima of function (11) in this region, it is proposed to introduce a mesh in the region 0 < Re ( ) < , 0 < Im ( ) < , where is a constant that defines the boundary of the search for eigenvalues [42]. The nodes of the introduced mesh are used as initial approximations for the numerical method of finding the minimum of the function of two variables (13). In Refs. [42] the Hook-Jeeves method [48] was used, but there are also other efficient numerical methods for zero-order multidimensional minimization [48].

Analysis of leaky modes in terms of inhomogeneous plane waves
Consider the complex values (4)-(6), explicitly distinguishing their real and imaginary parts, which will allow us to reformulate Eqs. (4)-(6) in terms of inhomogeneous waves, whose amplitude is also a function of coordinates and : Consider the expression (14) in the form of an inhomogeneous wave where ( , ) is the amplitude of the inhomogeneous wave defined as Consider in more detail the inhomogeneous wave in the form (17)   The region of existence of leaky modes corresponding to non-uniform waves with non-increasing amplitude is shown in Figure 2 (the cone between two dashed lines). Namely, if the wave vector of the non-uniform wave is located in the cone between two dashed lines, then this leaky mode has a non-increasing amplitude and can propagate in the positive direction of the -axis.
Let us consider the inhomogeneous wave (14) in more detail. As shown in Appendix (see (22)), in the region 0 ″ + 0 ″ < 0 of a non-uniform wave the amplitude will increase indefinitely, therefore we will consider the wave (14) in the region 0 ″ + 0 ″ ⩾ 0, see Figure 2. This exponential growth is real within a limited transverse distance surrounding the origin [49]. Using conservation of energy flux, one can show [49] that any mode that decreases exponentially as it propagates must increase exponentially transverse to the direction of propagation. However, it is evident that exponential growth of the field (and mode energy) that extends to infinity is unphysical since we have a finite energy source. A more detailed analysis of this problem is beyond the scope of our paper.

Numerical analysis of leaky modes of symmetric three-layer waveguides
Let us proceed to numerical analysis of the obtained representation of the leaky modes (14)- (16). Since the structure of the waveguide under consideration is symmetric, it is sufficient to consider only Eq. (14). We give in Figure 3 the calculated values of the complex phase deceleration coefficient calculated for a waveguide with =1.47, =1.565, = 0.55 and ℎ = 1.1 .  As seen from Figure 4, the field is concentrated in a cone formed by lines of constant amplitude, and the maximum intensity is observed at the boundaries of the cone where the amplitude is maximal. Outside the cone, there is an area of infinite growth of the amplitude of inhomogeneous waves represented by Eqs. (14)- (16); in this area the inhomogeneous waves characterizing the leaky mode in cladding layers cannot exist.
Inside the guiding layer, on the contrary, the field represented by Eqs. (14)- (16) attenuates rather rapidly and becomes almost completely damped at a distance of several wavelengths. The fields in the coating layer and the substrate are inhomogeneous waves, whose amplitudes decay exponentially the stronger, the smaller the distance to the waveguide layer. Due to the rapid attenuation of the field in the waveguide layer and the gradual removal (escape) of inhomogeneous waves, which characterize the behavior of the leaky mode in the cladding layers, one can observe a virtual "separation" of the leaky mode from the waveguide layer. A similar "separation" is also characteristic of higher leaky modes (see Figures 5, 6). We also note that experimental data on the propagation of leaky modes, which qualitatively agree with the results obtained by us, are given in [50] (see Figures 4(b), 5(b), 6(b)). In the experimental studies given in [50], leaky modes also propagate in a cone, and are also characterized by the presence of "separation" of the leaky mode from the waveguide layer.
It is commonly assumed that for the leaky modes the amplitude increases with the distance from the waveguide along the vertical axis (at a fixed longitudinal distance and in the absence of losses in the waveguide). However, as it propagates along the axis , this mode decays due to permanent energy losses from the waveguide layer to the environment. Functionally, the fields of leaky modes (vertical profile) are identical to the fields of ordinary guided modes; however, since unlike normal guided (homogeneous) modes, the leaky modes are inhomogeneous waves. In this regard, the representation of leaky waves of planar waveguides using the solutions of the wave equation seems to be preferable for us, compared to the traditionally used representation using the solutions of the Helmholtz equation.
At the same time, some features were revealed that we plan to analyze in our further publications. It is important to emphasize that the region of existence of leaky modes corresponding to inhomogeneous waves with non-increasing amplitude is detected (the cone between two dashed lines in Figure 2). Moreover, if the wave vector of an inhomogeneous wave is located in the region of the cone between two dashed lines, then such a leaky mode has a non-increasing amplitude and can propagate for a sufficiently long distance in the waveguide without absorption.

Conclusion
As is well known, conventional guided modes that exist when the waveguide layer thickness is above the critical value are considered in the optical beam representation as plane waves propagating in a regular waveguide due to the total internal reflection of waves at the interfaces between the waveguide media. From this point of view, the leaky waves propagate due to the effect of disturbed total internal reflection: during each act of disturbed total reflection at the interfaces of the media forming the waveguide, some of the power of this guided mode is radiated, i.e. "flows out" into the space surrounding the waveguide.
It is important to emphasize that the number of leaky modes with a gradual leakage is limited, unlike the continuum of radiative modes. The resulting gradual leakage waves form a discrete spectrum and are plane inhomogeneous waves. On the contrary, the radiative modes form a continuum (their spectrum is continuous) and are plane homogeneous waves. As a result, the replacement of one kind of these waves with another kind requires serious analysis in each specific case. From this point of view, the methods developed by us are undoubtedly useful for theoretical and numerical studies of dielectric and, in particular, optical waveguides supporting leaky modes, for example, when used as basic elements in the development of advanced sensors or various interface elements in integrated optical processors.
Lines of equal phase will be given by equations ′ + ′ = Const. Lines of equal amplitude will be given by equations ″ + ″ = Const. Lines of equal phase and equal amplitude will be orthogonal to each other in nonabsorbing media due to the fact that 2 + 2 = 2 , where is the wave number corresponding to the medium in which the wave propagates.

Transition to phase-ray coordinates
We introduce new coordinates attached to the lines of constant phase and constant amplitude: By virtue of the previously established orthogonality property (20), the introduced coordinate system will be orthogonal. In the coordinates , the form of the considered inhomogeneous plane wave (19) is considerably simplified: In the new variables, according to (22), the amplitude of the inhomogeneous wave under consideration decreases along the positive direction of the axis ,. For each fixed = > 0 the non-uniform wave is characterized by the amplitude ⋅ − and the behavior harmonic along (see Figure 7).
For each fixed = , the amplitude of the non-uniform wave decreases in the positive direction of the axis (see Figure 7). In the half-plane < 0, the amplitude of the non-uniform wave will increase indefinitely, therefore from the physical point of view the half-plane < 0 corresponds the so-called shadow region of the non-uniform wave under consideration.