Band structure and spatial electromagnetic field distribution in the 2D photonic crystals with nontrivial elementary cell

Periodic systems with the silver-based waveguides are discussed. The elementary cell in such systems consists of several waveguides. We employ the multiple scattering formalism based on solving the Maxwell equations solution for a single cylinder. The method involving the lattice-sum transformation into a fast convergent series is implemented. This enables us not only to increase the computational accuracy but also to reduce the necessary time for performing the numerical simulation. It is shown that under condition of nonzero propagation constant, the guided modes located closely in the frequency region appear. The number of these modes equals the number of waveguides in the elementary cell of the system. As a quasi-wave vector approaches the edge of the Brillouin zone, the frequency of these modes decreases. Increasing the propagation constant affects the modes shifting them towards the high frequency region. It is also shown that it is possible to match simultaneously wave vector and frequency of guided modes in the zigzag shaped systems.


Introduction
The optical waveguides are the components inherent in various optical and optoelectronics devices. In many cases they are important for the optical signal transmission between the different parts of the devices. However, due to tendency to miniaturization, the possible interaction between closely placed different waveguides inside the device becomes significant and leads to emergence of many phenomena. In particular, this concerns the plane arrays of identical long cylindrical waveguides. Such arrays belong to the family of low-dimensional photonic crystals. A material for the components of the photonic crystal can be dielectric (semiconducting) (see, e.g. [1]), metallic (see, e.g. [2]) and even superconducting [3].
It was shown in [4] and [5] that the quality factor for a whispering-gallery mode grows remarkably for the linear and circular arrays of dielectric microspheres. However, the isolated microsphere has no features of this kind. Then, the guided modes in the plane array of cylindrical waveguides are investigated in the numerous publications (see [6] and Refs. therein). In particular, it was discovered in [6] that the interaction between the waveguides inevitably results in the formation of the high quality guided modes. The circular arrays of dielectric rods reveal the specific features similar to that of the plane ones [7].
The optical pulse propagation in the metal nanoparticle chain waveguides is investigated firstly in [8]. Recently, the plasmon mode propagation in the plane array of the closely spaced metallic (silver) cylinders has been investigated in [9], a very specific case for the guided modes being considered. In that paper the propagation constant β, wavevector component parallel to the axes of the cylinders, is considered to be zero. However, the guided modes with β ≠ 0 are also important for various applied problems. In addition, they are relevant for simulating the optical phenomena which have analogues in the solid state physics, like the Bloch oscillations, Anderson localization, Bloch-Zener tunneling etc. [10][11][12][13][14][15].
Below in this paper we study the guided modes in the periodic systems with the elementary cell consisting of two metallic waveguides. These systems are found to possess the effects of anomalous transmission and refraction [16]. The similar effect was found experimentally in the system of silicon nanoposts [17]. However, due to plasmon frequencies located in the near ultraviolet, the metal waveguides are possible to operate at higher frequencies (including the visible part of spectrum) and smaller wavelengths. Thus the metal waveguides can be much smaller than their dielectric analogues. Let us point out that these systems (metallic and dielectric) traditionally were studied for the case of zero propagation constant so that the whole the whole class of guided modes was never considered. We show that the periodic systems with the elementary cell consisting of two metal cylinders reveal two kinds of dispersion curves with the unique properties that depend on the propagation constant value and geometrical parameters of elementary cell. We also determine the propagation constant β min for which all modes are located within the borders of transmission window and β max for which all modes are found in the visible part of spectrum.
The paper is organized as follows. In the next section we describe the system of parallel metallic cylinders. Then, we represent the key formulas of the multiple scattering formalism following Ref. [6]. This formalism is based on the formally exact solution for the scattering of the electromagnetic wave by a single infinite cylinder [18]. We derive a set of equations determining the dispersion law for the guided modes and describe the way of solving this set of equations. Then, we discuss the results obtained.

Multiple scattering formalism
The multiple scattering formalism is used for the detailed description of the electromagnetic field structure. It is based on the exact solution of the Maxwell equations for a single cylinder 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 constants are ε j and μ j , correspondingly. All the cylinders are parallel to the z-axis. We assume that there are no external current and charge in the system, so Maxwell's equations take form as Here J m (z) is the Bessel function, (ρ j , φ j , z) are the polar coordinates in 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. Unit vectors e j ρ , e j φ , e z form the basis of j-th waveguide's cylindrical framework, M jm and N jm are the vector cylindrical harmonics. Coefficients c jm , d jm are called the partial amplitudes. The prime denotes differentiation with respect to the whole argument of the function. The normalized constant J m (κ j R j ) in the denominator of expressions (2) is introduced to increase the accuracy of numerical simulations.
The electromagnetic field outside all the waveguides is a sum of fields formed by each waveguide. The contribution of j-th waveguide reads as follows     Here H (1) m (z) is the Hankel function of the first kind. The coefficients of expansion a jm , b jm are called the partial amplitudes. The normalized constant H (1) m (κ ex 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 partial amplitudes a jm , b jm , c jm , d jm can be determined via boundary conditions for the Maxwell equations on 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. So one can write: In the framework of the multiple scattering formalism only 4 conditions are used -the continuity of tangent (φ and z) components of fields E and H. The remaining two conditions are fulfilled automatically. We employ 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 (5) and equations (4), one can obtain a set of linear equations: After excluding the partial amplitudes c jm and d jm from expressions (6), we derive the system of equations that contain partial amplitudes a jm и b jm alone: Here we used the notation The relation between the solution of equation (7) and the partial amplitudes c jm and d jm , being excluded from (6), is given below One can determine the electromagnetic field produced by the waveguide system via solving the equation (7) and substituting this solution into formulas (9), (3) and (2).

Periodic waveguide system
Let us consider periodic waveguide system with an arbitrary elementary cell. We use multiindex Jj to encounter the cylinders. The first number J = 0, ±1, ±2 etc. stands for the elementary cell number, while the second j defines the number of the cylinder within the elementary cell and runs from 1 to N. Then the equation (7) Let us point out that 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 on the number j remains, as we consider the system with an arbitrary elementary cell. Next we introduce   Ll Jj l j LJ      a ρ ρ being the distance between waveguides (Jj) and (Ll), Ll Jj  is the angle between vectors a(L -J) +ρ lρ j and a (system's period).
The solution of equation (10) acquires the form of the Bloch wave characterized by the quasi-wave vector k The substitution of expression (11) into the equation (10) results in the equation It is possible to extract the lattice sum in equation (12) and transform it into a fast convergent series. This technique is very useful in the course of the numerical simulations, as the lattice sum arising naturally from Eq. (12) is a slow convergent one. Tens of thousands terms are necessary to be summed up in order to obtain the accuracy required. We employ the method involving the analytical properties of the Hankel functions and their integral representation. For the detailed description, we refer to paper [20]. Here we briefly demonstrate the main steps of that method.
First of all, let us rewrite the lattice sum on the left-hand side of equation (12). The second term in (13) contains the lattice sum in most common form, so we proceed with that one. The first term in (13) can be obtained from the lattice sum in the second one via subtraction the term with L=0 and calculating limit We also choose the coordinate system with the z-axis parallel to the waveguides axes and x-axis along the vector a.  The Hankel function of the first kind can be represented as follows [18]     The integral on the right-hand side of (14) is a convergent one provided the integration contour lies within the shaded region (Fig. 1). This is valid for contour C α . Let us assume that 0 < φ L < π, then we take α = φ L -π / 2. After substitution τ = t + α -π / 2 one can obtain: One can observe that the analytical properties of integral (14) are employed to choose the appropriate contour C α and compensate nontrivial part of phase φ L . Next, we change the integration variable so that sin t = y and move the integration contour to the real axis:  (17) and obtain the lattice sum in the rapidly convergent representation: The same calculations are performed for the case -π < φ L < 0. The only difference is that one should take α = φ L + π / 2 in order to keep the integral in (15) convergent. The expression for the lattice sum valid for all values of φ L reads as follows: Let us return to expression (13) and obtain the second kind of lattice sum: First, we assume that 0 < φ 0 < π, take the integral representation for the Hankel function from (16) and change the integration variable to μ so that y = y μ = (ka + 2πμ)/(κ ex a):   The same calculations are performed for the case -π < φ 0 < 0:  (22) One can observe that the lattice sum has the following symmetry: U -ν =e iππ U ν . This feature is obvious in expression (20). The half-sum of expressions (21) and (22) has the same symmetry. Since we assume that the subscript of the lattice sum is always positive, we may proceed with any of the expressions (21) or (22). The case ν = 0 will be discussed separately as it involves special technique.
Let us split the expression (21) into two parts: The first term in (23) does not need an exponential factor to be convergent. So, we keep the sum unvaried and calculate the integral: Again, the first term in (25) converges without exponential factor. As we did it before, we leave the sum unevaluated, and the calculation of the integral leads us to the result with the absence of factor e -iπν as compared to expression (24). The second term in (25) contains the difference and thus can be written in polynomials of y μ provided the parameter ν takes the integer values.
In order to calculate the limit in the second term in expression (25), we employ the Euler's formula [21] , 0 The same calculations can be performed in the case of odd ν: Finally, we use the relation (-1) n B n (-z)=B n (z)+nz n-1 and obtain the formula for the lattice sum U ν : Expression (29) is valid for the positive integer values of ν. For the negative integer values of ν, one should use the symmetry U -ν =(-1) ν U ν . The expression for the lattice sum with zero subscript can be derived via taking the limit ν→0 in any of the expressions (21) or (22). We omit the detailed description of these calculations and write down the result: Here C = 0.577… -is the Euler-Mascheroni constant. Let us return to the equation (12) that describes the partial amplitudes of the periodic waveguide system with an arbitrary elementary cell. Using the relations above, one can write

Numerical simulations
Let us demonstrate the relations obtained above. For the numerical simulation, we take the periodic system with the elementary cell consisting of two waveguides (Fig. 2). The geometrical parameters of the system are as follows: period a = 280 nm, the distance between the cylinders within each elementary cell b = a, radius of each cylinder R = 100 nm and angle φ = 90°. We choose silver as a material of the cylinders since it has the wide transmission window in the frequency region we are interested in. The silver permittivity is extracted from the experimental data in [22]. The refraction index of the surrounding medium is 1. As it was shown in [16], periodic systems with the similar parameters demonstrate the effect of anomalous refraction, and the investigation of the dispersion law of such systems is of high interest. All the calculations are carried out using MatLab, and we have employed the approximation which takes into account only terms with |m| < 20 in the field expansions (2) and (3). This approximation gives the accuracy of 1% in the frequency region we are interested in.  Figure 3 represents the dispersion law for the guided modes with the different propagation constant: a) -βa/π = 0, b) -βa/π = 0.66 (this corresponds to β min = 7.41 μm -1 ), c) -βa/π = 0.93 (this corresponds to β max = 10.43 μm -1 ). The horizontal dash-dotted line is the lower border of the transmission window. The upper border of the window is located near 3.9 eV and is not shown. The dashed line stands for the light cone that shows the dispersion law for a free photon: From Fig. 3 we see that for nonzero propagation constant, there exist closely located dispersion curves in the frequency region n ex ω < βc. The number of these curves equals the number of cylinders inside the elementary cell ( Fig. 3b and 3c; two curves with the lowest frequency). Unlike the other upper curves, these tend to be dispersionless. Let us point out that the similar effect is found in [23]: in case β≠0 the periodic system with one metallic cylinder in the elementary cell has one guided mode in the frequency region n ex ω < βc.
We have found the propagation constant β min = 0.66π/a for which all dispersion curves belong to the transmission window of the waveguide material ( fig. 3b). As β increases from β min to β max = 0.93π/a all dispersion curves shift towards high frequency region, and for β > β max they are located in the visible part of spectrum ( Fig. 3c; visible spectrum is occupied from 1.59 eV to 3.26 eV). In case β < β min the modes are located below the transmission window, and one can employ the Drude-Lorenz model for the low-frequency region even if no experimental data exist for the permittivity. This model was proven to be correct for silver, gold and copper [23].
In the case β = 0 the number of dispersion curves is also equal to the number of cylinders inside the elementary cell (Fig. 3a).  Let us consider relative shift of cylinder layers. In the terms of the parameters introduced in Fig. 2 this means that b and φ vary simultaneously so that the relative shift δ = b cos φ changes gradually while the distance between layers b sin φ = a remains constant. The period of the system and the radii of cylinders are the same. In the limiting case φ = 0, b = a/2 we have a single-layered system with one cylinder in the elementary cell. The actual period of the system is b, but we would assume it to be a = 2b. This would lead to system's spectrum to be symmetrical with respect to the point |k| = π/(2b). So, if ω(β, k) where -π/b < k < π/b -is the dispersion law of the system with period b and dispersion ω 1 (β, k 1 ) where -π/a < k 1 < π/a -is the dispersion law of the same system which period is considered to be a = 2b, then we have the following transformations: In order to see how the relative shift affects the dispersion law both for the guided modes that exist under condition β ≠ 0 alone and for the modes that still present when β = 0, we take the propagation constant β = β min . We see that it is possible to match simultaneously the quasi-wave vector at the edge of the Brillouin zone and the frequency of two upper dispersion curves (Fig. 4c). The corresponding value of quasi-wave vector is k = π/a, and the frequency ω = 2.02 eV. As β increases, the frequency of the point where curves merge also increases, and one can easily choose the appropriate value of β so that this point is located in the near infrared or visible part of spectrum. For β = β min , as is represented in Fig. 4c, the point belongs to the visible spectrum.

Discussion
The guided modes in the periodic systems of metallic cylinders with an arbitrary elementary cell have been studied. The effect of the propagation constant variation and relative shift of cylinder layers on the dispersion law in such systems is considered. It is shown that these waveguide systems have the guided modes of two kinds. The first one is the modes that exist under condition of nonzero propagation constant alone. They form a group of the closely located dispersion curves, and the number of these curves is equal to the number of cylinders within one elementary cell. The modes are practically dispersionless and their frequency gradually decreases as the quasi-wave vector runs over the Brillouin zone and approaches its edge. We also found the propagation constant β min = 7.41 μm -1 for which the modes enter the transmission window of the waveguide material.
The other kind of guided modes is formed by the curves that exist regardless from the value of the propagation constant; in particular, they exist even for β = 0. In the case β = 0 their number also equals the number of cylinders inside one cell. Unlike the modes of the first kind, they show dispersion and their frequency increases as the quasi-wave vector runs from zero to the edge of the Brillouin zone. The common feature of these modes is that they are shifted towards the high frequency region as the propagation constant increases and enters the visible part of spectrum for β > β max = 10.43 μm -1 . This fact can be used for realizing the mechanism of the mode selection in order to transmit at the desired frequency using the variation of characteristic mode frequency.
The relative shift of cylinder layers results in the fact that the dispersion curves of the second type merge when the quasi-wave vector is located at the edge of the Brillouin zone. The corresponding system parameters are b = a/2 0.5 , φ = 45 o (see Fig. 2), so the system has a zigzag shape. This feature can be used for the optical switching in the near infrared or visible part of spectrum (provided the propagation constant is carefully adjusted).
The paper is accomplished in the framework of Russian Foundation for Basic Research, grants 19-02-00433, 18-32-00204.