Theory for anomalous refraction of visible light by the plane array of metallic waveguides

The effect of anomalous or negative angle refraction in the periodic systems of silver waveguides is investigated. This effect arises solely due to the waveguide interaction and requires the system of specific geometrical configuration. It is also shown that the propagation constant provides a flexible mechanism to manipulate the anomalous refraction. We employ the multiple scattering formalism based on the exact solution of the Maxwell equations for a single waveguide problem. Also the method based on the lattice sum transformation into a fast convergent series is involved. This technique dramatically decreases the time necessary for performing the numerical simulations and increases the accuracy of the computational process. It is shown that all the directions, which the periodic system is possible to radiate, can be obtained beforehand, i.e. analytically, without numerical simulations. The latter ones are needed to determine energy the flux going into each of those directions.


Introduction
The optical waveguides are the most common components in various optoelectronics devices. They are designed for transmitting between the different structural units of the device. However, due to tendency to miniaturization, the possible interaction between the closely placed waveguides inside the device becomes significant and leads to the emergence of many phenomena that have analogues in solid state physics. In particular, these are the Bloch oscillations, Anderson localization, Bloch-Zener tunneling etc. [1][2][3][4][5][6].
The low-dimensional photonic crystal family is considered in many different publications. The material for these photonic crystals varies and can be dielectric (semiconducting) (see, e.g. [7]), metallic (see, e.g. [8]) or even superconducting [9]. This photonic crystal family is commonly represented by the plane arrays of identical long cylindrical waveguides. The guided modes in the plane array of cylindrical waveguides are investigated in the numerous publications (see [10] and Refs. therein). In particular, it was discovered in [10] that the high quality guided modes unavoidably result from the interaction between the waveguides. The circular arrays of dielectric rods reveal the specific features similar to that of the plane ones [11]. It was shown in [12] and [13] that the quality factor for a whispering-gallery mode grows remarkably for the linear and circular arrays of dielectric microspheres. However, an isolated microsphere has no features of this kind. * E-mail: iyppolishchuk@gmail.com It is possible to control the light reflection and refraction using specifically designed surfaces [14]. Below in this paper we consider the effect of anomalous refraction. This phenomenon was experimentally confirmed for the systems made of silica nanoposts [15] and theoretically studied for the metallic waveguides systems [16]. However, in both cases the propagation constant, wave vector component parallel to the axes of the waveguides, was assumed to be zero. Thus, the problem was effectively reduced from the 3D to 2D case. In this paper we also consider nonzero propagation constants, thus we truly solve the 3D problem. We show that it is possible to manipulate the diffraction channels via changing the propagation constant within certain limits: |β| < β cr . It can also be adjusted to achieve the optimal correlation between the incident wave vector and the period of the system. The paper is organized as follows. In the first section we describe the scattering problem for the system of parallel metallic waveguides. We invoke the multiple scattering formalism following Ref. [10]. This formalism is based on the formally exact solution for the scattering of the electromagnetic wave by a single infinite cylinder [17]. Then, we derive a set of equations determining the partial amplitudes for the field expansions and analytically obtain all the possible direction within which the system radiates. Finally, we discuss the results obtained.

Multiple scattering formalism
The electromagnetic field structure is described within the framework of multiple scattering formalism. This technique uses the exact solution for the Maxwell equations in a single waveguide problem. Let us consider system consisting of N parallel infinite cylinders surrounded by the medium with permittivity ε ex and permeability μ ex . The radius of the j-th cylinder is R j , permittivity and permeability are ε j and μ j , respectively. The z-axis is parallel to the waveguides. Inside each area where the permittivity and permeability are constants, Maxwell equations take the form of In equation (1) the symbols ε and μ stand for ε j and μ j if the region inside j-th cylinder is considered; otherwise they are equal to ε ex and μ ex , correspondingly. The field in the medium surrounding the cylinders is represented as a sum of external plane electromagnetic wave and the fields scattered by each cylinder. The contribution of j-th cylinder reads as follows: , , ρ φ ω ω ω κ ρ κ ρ μ ε κ β κ ρ κ Here H (1) m (z) is the Hankel function of the first kind, (ρ j , φ j , z) are the polar coordinates within the framework of j-th waveguide, and ω is the frequency of the guided mode. The The electromagnetic field inside j-th cylinder, being finite as ρ j → 0, expands as the following (see [10] for details): Here J m (z) is the Bessel function, (ρ j , φ j , z) are the polar coordinates within the framework of j-th waveguide, and ω is the frequency of the guided mode. The propagation constant β is the z-component of plasmon wave vector. The quantities M jm and N jm are the vector cylindrical harmonics (see explicit relations for M jm and N jm in [18]). The coefficients c jm , d jm are called the partial amplitudes and they fully describe the field inside the cylinder. The normalized coefficient J m (κ j R j ) in the denominator of expressions (3) is introduced to increase the accuracy of numerical simulations. All the other notations have the same meaning as in (2).
The external plane electromagnetic wave in the vicinity of j-th cylinder can be written in the following form: Like (2) and (3) the normalized coefficient J m (κ ex R j ) in the denominator of equations (4) is introduced to improve the accuracy of numerical simulations. The coefficients P jm and Q jm are called the partial amplitudes. It follows from expansions (4) that for scattering problem, the propagation constant β in relations (2) and (3) equals the z-component of the wave vector of the incident wave, k 0z .
To determine electromagnetic fields in the system, one has to know all partial amplitudes a jm , b jm , c jm , d jm . So, one writes down the boundary conditions for Maxwell equations at the surface of each cylinder. These conditions express the continuity of the tangent components of vectors E and H and the continuity of normal components of vectors B and D. The multiple scattering formalism we employ utilizes only 4 conditions, namely, the continuity of tangent (φ and z) components of fields E and H. The remaining two conditions are fulfilled automatically. According to (3), the field inside j-th cylinder expands in the framework of j-th waveguide. So, it is convenient to represent the field outside j-th waveguide in the same framework, and we use the Graph's theorem [19] to expand the electromagnetic field in the vicinity of j-th cylinder using its vector harmonics (here ρ lj , φ lj , are the polar coordinates of l-th cylinder with respect to j-th cylinder) Using Graph's theorem (6) and boundary conditions for Maxwell equations, one can obtain a set of linear equations: Eliminating the partial amplitudes c jm and d jm from expressions (7), we derive the system of equations that contain partial amplitudes a jm и b jm alone: The following notations are used in equation (8): The relation between the solution of equation (8) and the partial amplitudes c jm and d jm , eliminated from (7), is given below The solution of equation (8) and the formulae (10), (2), (3) determine the electromagnetic field in the whole cylinder system.

Periodic waveguide system
Next step we implement is to apply the relations obtained above to describing the plane wave scattering problem by the periodic waveguide system with an arbitrary elementary sell structure. Multi-index Jj is used to encounter the cylinders. The first number J = 0, ±1, ±2 etc. stands for the elementary cell number, while the second j determines the number of the cylinder within the elementary cell and runs from 1 to N. Then the equation (8) Matrix S j m does not depend on the number of the elementary cell due to the fact of periodicity. The same is valid for all the waveguide parameters like radius, permittivity and permeability. However, the dependence from the number j remains as we consider the system with an arbitrary elementary cell. Next we introduce The solution of equation (11) acquires the form of the Bloch wave characterized by quasi-wave vector k. In the case of scattering problem the quasi-wave vector is equal to the component of the incident wave vector parallel to the direction of periodicity in the system: After substitution the relations (12) into equation (11) The right-hand side of equation (13) does not depend on the number of the cell J. It follows from relations (4) that the partial amplitudes P Jj m and Q Jj m contain the factor e ikaJ thus cancelling out the dependency on the elementary cell number. One can also extract the so-called lattice sum on the left-hand side of equation (13) and convert it into a fast convergent series as it was shown in [20]. This technique is very useful and provides a great increase in accuracy while optimizing the computational time. Following the notations in [18], we rewrite the equation (13): , , κ κ π π κ κ κ π π κ (14) The lattice sum V ν reads as follows: Here Θ(z) denotes the Heaviside step function. Lattice sum U ν with the positive subscripts has the following form:  (16) For lattice sum U ν with negative subscript, we should use the symmetry U -ν = (-1) ν U ν . Here B n (z) is the Bernoulli polynomial. The expression for the lattice sum with zero subscript U 0 : In expression (17) C = 0.577… stands for the Euler-Mascheroni constant. The other notations used in (17) are the same as in (15).
Let us return to expressions (3) and write down the electromagnetic field scattered by all the cylinders in the Cartesian coordinates, z-axis being parallel to the cylinders axes and xaxis being along the vector a. Using this technique, one is able to obtain all the directions that the system is possible to radiate. We will demonstrate this method for the electric field scattered by the system. So combining (2) and (12), we have:  The lattice sum, extracted in (18), contains only the coordinates in the framework of cylinder j, located in elementary cell with number 0 while all non-coordinate dependencies are collected in the square brackets. This lattice sum can be transformed into a rapidly convergent series involving the same technique we described earlier. So, one has Before transformation (19) the lattice sum W ν was represented as a series over the cylinders with same number j. After that transformation we have W ν decomposed into a sum over plane waves with the wave vectors defined already analytically. Notice that μ in (19) actually enumerates the different plane waves as each wave vector {κ ex y μ , ±κ ex (1y μ 2 ) 1/2 , β} corresponds to unique value of μ (+ stands for reflection and -for refraction). This means we have beforehand all the directions in which the system is possible to radiate , i.e. without numerical calculations. The latter ones are needed to obtain the partial amplitudes and thus to get the coefficients for each plane wave in expansion (18). Parameter y μ depends only on the incident wave vector components and the period of the system. Thus the configuration of the elementary cell does not affect the number of diffraction channels.
Since μ takes integer values, only the finite number of discrete directions exists. One can easily obtain the number of the latter ones by solving the inequality As it follows from (20), κ ex a/π determines the number of integer solutions for inequality (20). That is, if κ ex a/π < 1, then only one integer solution is possible regardless from the value of ka/(2π), and we have only one possible channel for refraction and one possible channel for reflection. One way to change the number of available channels is to adjust the period a of the system. The larger the period a, the more directions are allowed. One also may affect this number by changing the propagation constant β as parameter κ ex incorporates it.

Numerical simulations
Let us apply the relations obtained above to the periodic system shown in Fig. 1. The incident wavelength is λ = 440 nm and incident angle is θ = 45° (angle between k 0 and yaxis). The wave is TE-polarized, i.e. magnetic field is parallel to the cylinder axis and electric field is orthogonal to the cylinder axis. The period of system is a = λ/2 1/2 = 311 nm, gap between the waveguides inside one cell is b = 2 nm and radius of each cylinder is R = 53 nm. All cylinders are made of silver. The permittivity data for silver is extracted from tables in [21]. The refraction index of the surrounding medium is 1. As it follows from (20), the system with the parameters mentioned above has 4 possible channels to radiate. These are schematically shown in Fig. 1b: 1 -anomalous reflection, 2 -normal reflection, 3normal refraction, 4 -anomalous refraction. Note that μ = -1 corresponds to both anomalous reflection and anomalous refraction while μ = 0 determines the normal reflection and normal refraction.
All calculations are carried out using MatLab and we use the approximation which takes into account only the terms with |m| < 20 in the field expansions (2) and (3). This approximation gives the accuracy of 1 %.  Figure 2 represents the spatial distribution of the z-component of magnetic field for the case β = 0: (a) incident wave, (b) response of the system and (c) -total field. We observe the response of the system to be predominantly located in the space below the cylinders. Thus, the total field is practically the same as incident field in the space above the cylinders while they differ dramatically in the space below the cylinders. Figure 3c shows the refracted wave to form negative angle with respect to the normal, i.e. the effect of anomalous refraction appears.
As it follows from (19) and (20), the wave vector of anomalously refracted wave is k 0 /2 1/2 {-1, -1, 0}. So, the refracted wave propagates at angle -45° to the normal. We have already mentioned that it is possible to manipulate the number of available diffraction channels via changing the period of the system or adjusting the propagation constant β. However, one cannot choose β randomly. If the propagation constant is sufficiently large, the channels with μ = -1 disappear. One can write the components of the incident wave vector in the spherical coordinates   Figure 3 represents the case when the anomalous channels are forbidden. We take φ = 50° so that the propagation constant exceeds the threshold: β ≈ 0.54k 0 > β cr . The other parameters are the same. As it follows from Fig. 3b, the system responds with the wave reflected in direction k 0 /2 1/2 {cos 50°, 1, sin 50°}, and the transmitted wave propagates in the direction of the incident wave in counterphase. This results in the normal reflection while the transmitted wave is greatly attenuated (Fig. 3c).
Let us point out that the anomalous refraction emerges when the incident wave vector and the period of the system correlate. That is, the incident wave covers the system so that the neighboring elementary cells are excited in counterphase. The corresponding relation reads