## Abstract

We obtain a convergent power series expansion for the first branch of the dispersion relation for sub-wavelength plasmonic crystals consisting of plasmonic rods with frequency-dependent dielectric permittivity embedded in a host medium with unit permittivity. The expansion parameter is *η*=*kd*=2*πd*/*λ*, where *k* is the norm of a fixed wavevector, *d* is the period of the crystal and *λ* is the wavelength, and the plasma frequency scales inversely to *d*, making the dielectric permittivity in the rods large and negative. The expressions for the series coefficients (also called dynamic correctors) and the radius of convergence in *η* are explicitly related to the solutions of higher order cell problems and the geometry of the rods. Within the radius of convergence, we are able to compute the dispersion relation and the fields and define dynamic effective properties in a mathematically rigorous manner. Explicit error estimates show that a good approximation to the true dispersion relation is obtained using only a few terms of the expansion. The convergence proof requires the use of properties of the Catalan numbers to show that the series coefficients are exponentially bounded in the *H*^{1} Sobolev norm.

## 1. Introduction

Sub-wavelength plasmonic crystals are a class of *meta-material* that possesses a microstructure consisting of a periodic array of plasmonic inclusions embedded within a dielectric host. The term ‘sub-wavelength’ refers to the regime in which the period of the crystal is smaller than the wavelength of the electromagnetic radiation travelling inside the crystal. Many recent investigations into the behaviour of meta-materials focus on phenomena associated with the *quasi-static limit* in which the ratio of the period cell size to wavelength tends to zero. Sub-wavelength microstructured composites are known to exhibit effective electromagnetic properties that are not available in naturally occurring materials. Investigations over the past decade have explored a variety of meta-materials, including arrays of microresonators, wires, high-contrast dielectrics and plasmonic components. The first two, especially in combination, have been shown to give rise to unconventional bulk electromagnetic response at microwave frequencies (Pendry *et al.* 1998, 1999; Smith *et al.* 2000) and, more recently, at optical frequencies (Povinelli 2009), including negative effective dielectric permittivity and/or negative effective magnetic permittivity. An essential ingredient in creating this response is local resonance contained within each period owing to extreme properties such as high conductivity and capacitance in split-ring resonators (Pendry *et al.* 1999).

In the case of plasmonic crystals, the dielectric permittivity *ϵ*_{p} of the inclusions is frequency-dependent and negative for frequencies below the plasma frequency *ω*_{p},1.1Shvets & Urzhumov (2004, 2005) have investigated plasmonic crystals in which *ω*_{p} is inversely proportional to the period of the crystal and for which both inclusion and host materials have unit magnetic permeability. They have proposed that simultaneous negative values for both an effective *ϵ* and *μ* arise at sub-wavelength frequencies that are quite far from the quasi-static limit, that is,1.2is not very small, where *d* is the period of the crystal, *k* is the norm of the Bloch wavevector and *λ* is the wavelength. In this work, we present rigorous analysis of this type of plasmonic crystal by establishing the existence of convergent power series in *η* for the electromagnetic fields and the first branch of the associated dispersion relation. The effective permittivity and permeability defined according to Pendry *et al.* (1999) are shown to be positive for all *η* within the radius of convergence *R*, and, in this regime, the extreme property of the plasma produces no resonance in the effective permittivity or permeability. This regime is well distanced from the resonant regime investigated in Shvets & Urzhumov (2004, 2005).

The analysis shows that the radii of convergence of the power series are at least *R*_{m}, which is not too small, as shown in table 1, which contains values of *R*_{m} for circular inclusions of various radii rd. The number *R*_{m} can be put in physical perspective by fixing the cell size and introducing the parameters *λ*_{m} and *k*_{M} such that the power series describes wave propagation for wavelengths above *λ*_{m} and wavenumbers below *k*_{M}. Table 2 contains values of *λ*_{m} and *k*_{M} when *d*=10^{−7} m. The wavelengths lie in the infrared range and the plasma frequency is *ω*_{p}=10^{15} s^{−1}.

We focus on harmonic H-polarized electromagnetic waves in a lossless composite medium consisting of a periodic array of plasmonic rods embedded in a non-magnetic frequency-independent dielectric host material. Each period can contain multiple parallel rods with different cross-sectional shapes; however, the rod–host configurations are restricted to those with rectangular symmetry, i.e. configurations invariant under a 180° rotation about the centre of the unit cell. The regime of interest for this investigation is that in which

— the plasma frequency

*ω*_{p}is high;— the ratio of the cell width to the wavelength is small (

*η*≪1).

From the formula , it is seen that a high plasma frequency *ω*_{p} gives rise to a *large and negative dielectric permittivity* *ϵ*_{p} in the plasmonic inclusions (figure 1). Following Shvets & Urzhumov (2004), the plasma frequency is related to the cell size byThis results in the relationwhere *c* is the speed of light in vacuum. The governing family of differential equations for the magnetic field is the Helmholtz equation with a rapidly oscillating coefficient1.3in which *A* is the matrix defined on the unit period of the crystal by*I* is the identity matrix andThe coefficient *A* is not coercive in the regime *ω*_{p}≥*ω*, as *ϵ*_{p} is negative in this regime, and it is precisely the appearance of negative *ϵ*_{p} that allows us to obtain a *convergent* power series expansion of the electromagnetic field and the frequency *ω*^{2} for a fixed Bloch wavevector , with .

In theorem 5.2, we obtain the following series expansion for the frequency *ω*^{2}:1.4in which is a tensor of degree 2*m*+2 in . This gives rise to a convergent power series for an effective index of refraction defined through1.5The effective property is well defined for all *η* in the radius of convergence and is not phenomenological in origin but instead follows from first principles using the power series expansion. Interpreting the first term of this series as the quasi-static index of refraction *n*^{2}_{qs}, the remaining terms then provide the dynamic correctors of all orders. In §6, we define the effective permeability *μ*_{eff} and prove that both and *μ*_{eff} are positive for *η* in some interval (0,*η*_{0}] and that a mild effective magnetic response emerges for the homogenized composite, even though the component materials are non-magnetic (). Having defined and *μ*_{eff}, the effective electrical permittivity *ϵ*_{eff} can be defined through the equationso that *ϵ*_{eff} is positive whenever both and *μ*_{eff} are positive. Thus, one has a solid basis on which to assert that plasmonic crystals function as materials of positive index of refraction in which both the effective permittivity and permeability are positive. The method developed here can be applied to other types of frequency-dependent dielectric media such as polaratonic crystals. From a physical perspective, this work provides the first explicit description of Bloch wave solutions associated with the first propagation band inside nano-scale plasmonic crystals. In the context of frequency-independent dielectric inclusions, the first two terms of are identified via Rayleigh sums in McPhedran *et al.* (1996).

To emphasize the difference between effective properties defined for meta-material structures where the crystal period *d* is fixed and effective properties defined in the quasi-static limit, i.e. *k* fixed and , we refer to the latter as quasi-static effective properties and denote these with the subscript qs. The situation considered in this paper contrasts with the case in which *ϵ*≈*d*^{−2} in the inclusion and is large and positive, investigated by Bouchitté & Felbacq (2004). In that case for , *μ*_{qs}(*ω*) has poles at Dirichlet eigenvalues of the inclusion and therefore is negative in certain frequency intervals (see also Bouchitté & Felbacq 2004, 2005). In fact, what allows us to prove convergence of the power series in the plasmonic case is precisely the absence, owing to negative *ϵ*_{p}, of these internal Dirichlet resonances.

In the regime where *ϵ*_{p} is negative and large, the perturbation methods used for describing Bloch waves in heterogeneous media developed in Odeh & Keller (1964), Conca *et al.* (2006) and Bensoussan *et al.* (1978) cannot be applied. Our analysis instead makes use of the fact that *ϵ*_{p} is negative and large for sub-wavelength crystals and develops high-contrast power series solutions for the nonlinear eigenvalue problem that describes the propagation of Bloch waves in plasmonic crystals. The convergence analysis takes advantage of the iterative structure appearing in the series expansion and is inspired by a technique that Bruno (1991) developed for series solutions to quasi-static field problems. We prove that the series converges to a solution of the harmonic Maxwell system for ratios of cell size to wavelength that are not too small. Indeed, for typical values of the plasma frequency *ω*_{p}, the analysis delivers convergent series solutions for nano-scale plasmonic rods at infrared wavelengths.

In §6, we compute the first two terms of the dispersion relation for circular inclusions (Shvets & Urzhumov 2004, 2005) and provide explicit bounds on the relative error committed upon replacing the full series with its first two terms. The error is seen to be less than 3 per cent for values of *η* up to 20 per cent of the convergence radius, so that the two-term approximation provides a numerically fast and accurate approximation to the dispersion relation.

The high contrast in *ϵ* gives rise to effective constants *ϵ*_{eff} and *μ*_{eff}. In the bulk relation1.6where *B*_{eff} is the average over the period cell (a flux), whereas *H*_{eff} is the average of *H*_{3} over line segments in the matrix parallel to the rods. Taking the ratio of *B*_{eff}/*H*_{eff} delivers an effective magnetic permeability and one recovers magnetic activity from meta-materials made from non-magnetic materials. This phenomenon was understood by Pendry *et al.* (1999), and has been made rigorous in the quasi-static limit through two-scale analysis in several cases. These include the two-dimensional arrays of inclusions in which *ϵ* scales as *d*^{−2} (Pendry *et al.* 1999; Bouchitté & Felbacq 2004, 2005; Felbacq & Bouchitté 2005*a*,*b*), two-dimensional arrays of ring resonators whose surface conductivity scales as *d*^{−1} (Kohn & Shipman 2008) as well as three-dimensional arrays of split-ring wire resonators in which the conductivity scales as *d*^{−2} (Bouchitté & Schweizer in press). This ‘non-standard’ homogenization has been understood for some decades in problems of porous media and imperfect interface (Cioranescu & Paulin 1979; Auriault & Ene 1994; Lipton 1998; Zhikov 2000; Donato & Monsurró 2004) and recently has given rise to interesting effects in composites of both high contrast and high anisotropy (Cherednichenko *et al.* 2006; Smyshlyaev 2009).

The two-scale analysis in these cases relies on the coercivity of the underlying partial differential equations (PDEs). The problem of plasmonic inclusions, however, is not coercive because *ϵ* is negative in the plasma—but it is precisely this negative index that underlies the convergence of the power series. As we shall see, the uniqueness of the solution of the Dirichlet problem for △*u*−*u*=0 in the plasmonic inclusion gives exponential bounds on the coefficients of the series, which allows us to prove that it converges to a solution of the differential equation (1.3). This result presented here shows that, by considering a finite number of terms in the series, one has an approximation of the true solution, to any desired algebraic order of convergence. In this context, we point out the recent works (Smyshlyaev & Cherednichenko 2000; Kamotski *et al.* 2007; Panasenko 2009) which show that the power series expressed by the formal two-scale expansion of Bakhvalov & Panasenko (1989) is an asymptotic series in certain cases under the hypothesis that the coefficient *A* is coercive.

## 2. Mathematical formulation and background

We introduce the nonlinear eigenvalue problem describing the propagation of Bloch waves inside a plasmonic crystal and provide the context for the power series approach to its solution.

For points ** x**=(

*x*

_{1},

*x*

_{2}) in the

*x*

_{1}

*x*

_{2}-plane, the

*d*-periodic dielectric coefficient of the crystal is denoted by

*ϵ*(

*ω*,

**), whereBoth materials are assumed to have unit magnetic permeability, .**

*x*We assume a Bloch-wave form of the field, where is the unit vector along the direction of the travelling wave and *k*=2*π*/*λ* is the wavenumber for the wavelength *λ*. The magnetic and electric fields are denoted by ** H**=(

*H*

_{1},

*H*

_{2},

*H*

_{3}) and

**=(**

*E**E*

_{1},

*E*

_{2},

*E*

_{3}), respectively. For

*H*-polarized time-harmonic waves, the non-vanishing field components are2.1in which the fields

*H*

_{3}(

**),**

*x**E*

_{1}(

**) and**

*x**E*

_{2}(

**) are continuous and**

*x**d*-periodic in both

*x*

_{1}and

*x*

_{2}. The Maxwell equations take the form of the Helmholtz equation (1.3), in which substitution of gives2.2and2.3where

*H*

_{3}satisfies the transmission conditions on the interface between the rods and host material given by2.4Here, the subscripts indicate the side of the interface where the quantities are evaluated and

**is the unit normal vector to the interface pointing into the host material. We denote the unit vector pointing along the**

*n**x*

_{3}direction by

*e*_{3}, and the electric field component of the wave is given by

**=−(**

*E**ic*/

*ωϵ*)

*e*_{3}×∇

*H*

_{3}.

For each value of the wavevector ** k**, equations (2.2)–(2.4) provide a nonlinear eigenvalue problem for the pair

*H*

_{3}and

*ω*

^{2}. One of the main results of this work is to show that this problem is well posed by explicitly constructing solutions using power series expansions. In order to develop the appropriate expansions, we rewrite equations (2.2)–(2.4) in terms of

*η*and a dimensionless variable

**in that normalizes a period cell to the unit square**

*y**Q*=[0,1]

^{2},

**=**

*x*

*y**d*=

*y**η*/

*k*. We define the

*Q*-periodic functionand, for convenience of notation, we redefine2.5to arrive at the eigenvalue problem that requires the pair

*h*(

**) and**

*y**ω*

^{2}to be a solution of the master system2.6We prove in theorem 5.2 that this eigenvalue problem can be solved by constructing explicit convergent power series solutions.

The development of the remainder of this paper is as follows. In §3, the power series expansion is introduced and the associated boundary-value problems (BVPs) necessary for determining each term in the series are obtained. The BVPs are given by a strongly coupled infinite system of linear partial differential equations. The existence and uniqueness of the solution to this infinite system is proved under fairly general hypotheses in §4. Because the system is coupled through convolution products, the convergence analysis is delicate. The convolutions are handled through estimates involving sequences of Catalan numbers whose convolution products determine the next element of the sequence. The Catalan numbers and their relevant properties are discussed and used to derive bounds on the Sobolev norm of each term of the series expansion in §5. These bounds are then used to establish the radius of convergence for the power series representations of the field and frequency (the dispersion relation), which solve the nonlinear eigenvalue problem (2.6). Section 6 deals with the computation of error bounds for finite-term approximations of the magnetic field and the dispersion relation.

## 3. Power series expansions

We take *η* to be the expansion parameter for the field *h*(** y**) and the frequency

*ω*

^{2},3.1in which the functions

*h*

_{m}are periodic with period cell

*Q*.

Inserting equation (3.1) into equation (2.6) and identifying coefficients of like powers of *η* on the right- and left-hand sides yields the equationsfor *m*=0,1,2,…, in which *h*_{m}≡0 and for *m*<0 and the terms involving the subscript ℓ are convolutions written according to the following summation conventions:where [*n*/2] denotes the largest integer less than or equal to *n*/2. The BVP satisfied by *h*_{0} in isso that this function is necessarily a constant in . We denote this constant value by . It will be convenient to use the dimensionless parameters *ψ*_{m} and defined through3.2In terms of *ψ*_{m} and , the above equations for the functions *h*_{m} become3.3for *m*=0,1,2,…, in which *ψ*_{m}≡0 and for *m*<0. We thus have an infinite system that the sequences and must satisfy. The system is written in terms of a Poisson equation in with Neumann boundary data and a Helmholtz equation in *P* with Dirichlet boundary data, namely3.4and3.5in which3.6where and *ψ*_{m}|_{p} indicate the trace of *ψ*_{m} on the and *P* sides of the interface ∂*P* separating the two materials. The Neumann BVP (3.4) is subject to the standard *solvability condition* given by3.7Here, the area integral over the domains and *P* are denoted by and 〈·〉_{P}, while the line integral over the interface ∂*P* separating the two materials is denoted by 〈·〉_{∂P}.

The iterative algorithm for solving the system is as follows. First note that, from the definition of *ψ*_{0}, it follows that *ψ*_{0}=1 for ** y** in . The function

*ψ*

_{0}is determined inside

*P*by solving equation (3.5) with Dirichlet boundary data on ∂

*P*. Then,

*ψ*

_{1}on is the solution of equation (3.4) with Neumann data on ∂

*P*. The process then continues with the boundary values on ∂

*P*of

*ψ*

_{m}in providing the Dirichlet data for

*ψ*

_{m}in

*P*, which, in turn, provides the Neumann data for

*ψ*

_{m+1}in , up to an additive constant. The term is determined by the consistency condition (3.7), and an inductive argument can be used to show that it is a monomial of degree

*m*in . The equations satisfied by

*ψ*

_{0},…,

*ψ*

_{4}inside ,

*P*and ∂

*P*are listed in table 3.

Note that in table 3 (meaning for ). In the next section, we identify a large class of shapes for the plasmonic rod cross sections for which the sequences and satisfy the infinite system (3.4)–(3.7) and , *m*=1,2,…. The mean zero property of *ψ*_{m} on provides a tractable scenario for proving the convergence of the resulting power series. This topic is discussed further in §4.

In what follows, we will make use of the equivalent weak form of the infinite system. To introduce the weak form, we introduce the space of complex-valued square integrable functions with square integrable derivatives *H*^{1}(*Q*). For *u* and *v* in *H*^{1}(*Q*), the inner product is defined by and the norm is given by . The *H*^{1} inner products and norms over *P* and are defined similarly.

The weak form of the infinite system is given in terms of the space *H*^{1}_{per}(*Q*) of functions in *H*^{1}(*Q*) that take the same boundary values on opposite faces of *Q*. The weak form of the system (3.4)–(3.7) is given by3.8for all , where and . The equivalence between equations (3.4)–(3.6) and the weak form follows from integration by parts and the solvability condition (3.7) follows from equation (3.8) on choosing the test function *v*=1 in equation (3.8).

## 4. Solutions of the infinite system for plasmonic domains with rectangular symmetry

The goal here is to identify solutions of the infinite system for which one can prove convergence of the associated power series with a minimum of effort. Looking ahead, we note that the convergence proof is expedited when one can apply the Poincaré inequality to the restriction of *ψ*_{m} on for *m* greater than some fixed value. To this end, we seek a solution (, ) such that for *m*≥1 one has and the sequences satisfy equations (3.4)–(3.7) or equivalently satisfy equation (3.8). We show that we can find such solutions for the class of plasmonic domains *P* with rectangular symmetry. Here, we suppose that the unit period cell is centred at the origin and the class of rectangular symmetric domains is given by the set of all shapes invariant under 180° rotations about the origin. This class includes simply connected domains such as rectangles and ellipses as well as multiply connected domains. For these geometries and for each *m*=1,2,3…, it is demonstrated that one can add an arbitrary constant to the restriction of the function *ψ*_{m} on without affecting the solvability condition (3.7). Under the assumption of rectangular symmetry for the inclusion *P*, we will show that there exists a pair (, ) satisfying equation (3.8) with the functions *ψ*_{m} in the subspace of real-valued functions with zero average in .

We now record the symmetries necessarily satisfied by any solution to equations (3.4)–(3.6) for plasmonic domains with rectangular symmetry. We denote the dependence of *ψ*_{m} on the unit vector by writing so that

;

.

Statement (i) is true for inclusions of arbitrary shape, while statement (ii) is true only for inclusions with rectangular symmetry. Taken together, these statements imply that so that is even or odd in *Q* accordingly as the index *m* is even or odd. From its definition, *ψ*_{0}≡1 in and trivially satisfies the solvability condition (3.7). The solvability of *ψ*_{m} when *m*≥1 is proved by induction on *m* using the weak form (3.8). We have the following theorem.

### Theorem 4.1

*For each* *, there exists a sequence of functions* *,* *, and a sequence of real numbers* *, with* *, solving the weak form (*3.8*) for each integer m.*

### Proof.

The proof is divided into the base case (*m*=1 and 2) and the inductive step.

*Base case.* The solvability for *ψ*_{1} and *ψ*_{2} can be established without the need to restrict to rectangular symmetric inclusions. This restriction will be necessary only in the inductive step. Setting *m*=1 and *v*≡1 in equation (3.8), we see that the left-hand side of equation (3.8) vanishes. This establishes the solvability for *ψ*_{1}. If we then take , we have a solution . Setting *m*=2 and *v*≡1 in equation (3.8), we obtainAs 〈*ψ*_{0}〉_{Q}>0 (see appendix, electronic supplementary material) and , this is one equation in one unknown . Solving for then gives Choosing this value for and also taking , we have a solution .

*Inductive step.* Let 2*n* be an even positive integer and assume that equation (3.8) has solutions for *m*=1,2,…,2*n*, with and . Then equation (3.8) has solutions for *m*=2*n*+1,2*n*+2 with and .

The solvability condition for *ψ*_{2n+1} is obtained by setting *v*=1 and *m*=2*n*+1 in the weak form, namelyThe hypothesis , odd≤2*n*−2, will imply that the convolutions *σ*_{m}, *m*≤2*n*−2, have the same even/odd property as the functions *ψ*_{m}. Indeed, writing out *σ*_{2n−3}, we have and as 2*n*−3−ℓ is odd when ℓ is even, it follows that *σ*′_{2n−3} is a linear combination of odd functions and is, therefore, an odd function. The same reasoning applies to all the other convolutions of index less than or equal to 2*n*−2. Moreover, is an odd function, as *σ*′_{2n−2} is even. Thus, all integrals in the consistency condition above vanish (for the integrals in , we can also use the fact that all functions belong to ), except thatAs 〈*ψ*_{0}〉_{Q}>0 (see appendix, electronic supplementary material), the solvability condition for *ψ*_{2n+1} is simply . We thus take to establish the existence of *ψ*_{2n+1}=0. Moreover, as *ψ*_{m} and *ξ*_{m−2} are real by the induction hypothesis, 0≤*m*≤2*n*, it follows that *ψ*_{2n+1} is real-valued. Thus, taking , we have a solution . Also, *ψ*_{2n+1} is an odd function as its index is odd. We now proceed to the solvability of *ψ*_{2n+2}, namelyAll terms in the above equation are real numbers, as we assumed *ψ*_{m} and real for 0≤*m*≤2*n*, with , and we just took and *ψ*_{2n+1} is real-valued. Thus, this equation contains the only one undetermined term . Thus, we have one real equation with one real variable, so that, taking to be such as to solve this equation and also taking , we complete the proof of the inductive step. ■

## 5. Proof of convergence

In this section, we show that the power series , and , where and *p*_{m}=∥*ψ*_{m}∥_{H1(P)}, converge and provide lower bounds on their radius of convergence. This will then be used to show that the pair and is a solution to equation (2.6). In §5*a*, we present the Catalan bound, which is used to provide a lower bound on the radius of convergence of the power series. In §5*b*, we derive inequalities that bound , *p*_{m} and in terms of lower index terms. In §5*c*, we present the properties of the Catalan numbers relevant for bounding convolutions of the kind appearing in equations (3.4) and (3.5). In §5*d*, we use an inductive argument on the inequalities of §5*b* to prove the Catalan bound. Finally, in §5*e*, we prove that the pair *h*_{η} and is a solution to the eigenvalue problem (2.6).

### (a) The Catalan bound

The following theorem is one of the central results of this paper.

### Theorem 5.1 (Catalan bound)

*For every integer m, we have that*5.1*in which C*_{m} *is the mth Catalan number,* *and* *, where the numbers J*_{1} *and J*_{2} *are determined as follows: J*_{1} *is the smallest value of J such that equation (5.1) holds for m≤4 and J*_{2} *is the smallest value of J for which the following polynomials Q***,R***,S** *in the variable J*^{−1} *are all less than unity:**The constants A,* *, β,* *,* *,* *,* *and p*_{2} *are determined by the particular choice of inclusion, while E(4)=16C*_{2}*/C*_{5}*≤0.7619.*

All bounds obtained here are expressed in terms of the Catalan numbers, area fractions and geometric parameters that appear in the Poincaré inequality and in an extension operator inequality. We start by listing these parameters and give the background for their description. It is known by Nečas (1967) that any function *ϕ* can be extended into *P* as an *H*^{1}(*Q*) function *E*(*ϕ*) such that *E*(*ϕ*)=*ϕ* for ** y** in and5.2where

*A*is a non-negative constant and is independent of

*ϕ*depending only on

*P*. For general shapes,

*A*can be calculated via numerical solution of a suitable eigenvalue problem. Constants of this type appear in Bruno (1991) for high-contrast expansions of the DC fields inside frequency-independent dielectric media. The second constant is the Poincaré constant given by the reciprocal of the first non-zero Neumann eigenvalue of and we have that . The last two geometric constants appearing in the bounds are the volume fractions

*θ*

_{P}and of the regions

*P*and . Using that

*C*

_{m}≤4

^{m}(§5

*c*), theorem 5.1 shows that , and are convergent for

*η*≤1/4

*J*, so that one may prove the following theorem.

### Theorem 5.2 (solution of the eigenvalue problem)

*Let R=1/4J, where J is the number prescribed by theorem 5.1. Then,* *converges as a series of real numbers and* *converges in the H*^{1}*(Q) Sobolev norm for η≤R. Moreover,* *and* *satisfy the eigenvalue problem given by the master system (*2.6*) (or by equation (*1.3*)).*

### (b) The , *p*_{m} and inequalities—stability estimates

We now derive the inequalities that bound , *p*_{m} and in terms of lower index terms. These inequalities follow from stability estimates for equations (3.4)–(3.6).

### Theorem 5.3

*Let m≥0 be an integer. Then*5.3*where the* *inequality holds for m≥2 only.*

Here, we have introduced the notation

### Proof.

We start by proving the *p*_{m} inequality. Recalling that equation (3.5) is satisfied by *ψ*_{m} in *P* giveswhere . Write the orthogonal decomposition *ψ*_{m}=*u*_{m}+*v*_{m}, where5.4andWe then have by the triangle inequality that5.5

The term ∥*u*_{m}∥_{H1(P)} is bounded using equation (5.2) and5.6to obtain5.7Here, equation (5.6) follows from the fact that the solution of equation (5.4) minimizes the *H*^{1}(*P*) norm over all functions with the same trace on ∂*P*. The term ∥*v*_{m}∥_{H1(P)} can be bounded using a direct integration by parts on the BVP for *v*_{m}5.8Now,where *p*_{m}=∥*ψ*_{m}∥_{H1(P)}. Using equations (5.7) and (5.8) in equation (5.5) gives5.9and the *p*_{m} inequality is established. We now prove the inequality. In the weak form (3.8), set *v*=*ψ*_{m} in and *v*=*u*_{m} in *P* to obtainWe now use the Cauchy–Schwarz inequality on the product of integrals appearing in each individual term. For the convolutions, we obtainwhere we used that . For the double convolutions, we obtainProceeding similarly with the other terms, we obtain5.10As the functions *ψ*_{m} have zero average in , we have the Poincaré inequality5.11where the constant can be computed from the Rayleigh quotient characterization of the first positive eigenvalue for the free membrane problem in . A simple computation using equation (5.11) then gives where . Using this inequality (5.10) gives5.12It will turn out to be to our advantage to apply equation (5.12) to the last term in equation (5.12) so as to replace it with5.13Using equation (5.13) in equation (5.12) yields the inequality in equation (5.3), valid for *m*≥2 (for *m*=1, use equation (5.12)).

Finally, we establish the inequality. Setting *v*=1 in the weak form (3.8), we obtain5.14(Recall that, for *m* odd, each term on the left-hand side of the above equation vanishes individually.) Solving for , we then obtain5.15where . We shall be using this equality for *m*≥5 only, so that and . Moreover, using that and where *θ*_{P} and denote the volume fractions of the regions *P* and , we have that and . Thus, proceeding with equation (5.15) as we did in the previous stability estimates, we obtainAs the iteration scheme at each step involves *p*_{m} and and , we adjust subscripts in the above inequality to obtain the inequality5.16 ■

### (c) The Catalan numbers

We briefly present some facts about the Catalan numbers that will be used in the sequel. The Catalan numbers *C*_{m} are defined algebraically through the recursion5.17These numbers arise in many combinatorial contexts (Lando 2003) as well as in the study of fluctuations in coin tossing and random walks (Feller 1968). It can be shown that , and a simple computation then gives the ratio of successive Catalan numbers5.18The first inequality above provides the exponential bound5.19It will be convenient to introduce the notation . From equation (5.18), it is clear that is decreasing in both *m* and *k*. In §5*d*, we shall make use of table 4, in which values of are listed.

#### (i) The even part of the Catalan convolution

The fact that needs to be taken into account in order to provide a suitable upper estimate on the incomplete convolution term *q*^{′*}_{m−1} appearing in the inequality in equation (5.3). Thus, we consider the convolution *C*_{n−ℓ}*C*_{ℓ} with the odd values of the index ℓ omitted and denote it by . We then define the even part *E*(*n*) by5.20The following lemma gives the estimate *E*(*n*)≤*E*(4), *n*≥4.

### Lemma 5.4

*The following two statements are true for all m*≥0: (i) *E*(2*m*) *is a decreasing sequence and* (ii) *E*(2*m*+1)=1/2. *Thus, for all m*≥4, *we have that* .

### Proof.

Statement (ii) is actually just an observation, as one can see by writing out the sum *C*_{n−ℓ}*C*_{ℓ}. Statement (i) can be deduced from the identity *C*_{2m−2ℓ}*C*_{2ℓ}=4^{m}*C*_{m} (Koshy 2008). Indeed, dividing both sides of this identity by *C*_{2m+1}, we obtainFrom equation (5.18), each of the above fractions 4*C*_{m}/*C*_{m+ℓ}, ℓ=1,2,…,*m*+1 is a decreasing sequence in *m* so that their product is also decreasing in *m*. This completes the proof. ■

### (d) Proof of the Catalan bound

### Proof.

(Catalan bound, theorem 5.1) Fix the values of the geometric parameters *A*, , *θ*_{P} and in equation (5.3). Starting with the initial estimates and , the inequalities (5.3) can be used recursively for *m*= 1, 2, 3, 4 to determine a number *J*_{1} such thatWe now proceed inductively; assume that5.21where *m*≥5. We then get for the single convolutionswhere and lemma 5.1 was used to introduce the factor *E*(4). Similarly, for double convolutions, we getwhere the factor *E*(4)^{2} comes from using lemma 5.1 twice. For the non-convolution terms, we get The same bounds hold for the terms , and , so that we have5.22The proof now essentially consists of applying these bounds to all terms in inequalities (5.3). The factor *J*^{−k} appearing on the right-hand side of each inequality is the workhorse of the proof: by taking *J* sufficiently large, it will allow us to close the induction argument. The incomplete convolution term *q*^{′*}_{m−1} presents special difficulties, as attempting a bound of the type (5.22) for this term does not produce a factor of *J*^{−k} (actually, it produces *J*^{0}=1).

Recall the inequality from equation (5.3)

Using equation (5.22) on this inequality gives5.23where *Q*_{m} is the following polynomial in *J*^{−1}:As we shall be using this inequality for *m*≥5 only, table 4 can be used to bound the numbers , so that we may write *Q*_{m}≤*Q**, where5.24The strategy now is to determine similar polynomials *R*_{m} and *S*_{m−1} for the other two inequalities, that is, *p*_{m}≤*R*_{m} *βJ*^{m}*C*_{m} and , and then take *J* large enough that all three polynomials are less than unity, allowing us to complete the induction argument. Having obtained *Q*_{m}, it is straightforward to obtain *R*_{m}. Indeed, using equations (5.22) and (5.23), the *p*_{m} inequality in equation (5.3) yields *p*_{m}≤*R*_{m} *βJ*^{m}*C*_{m}, whereThus, *R*_{m}≤*R**, where5.25

The inequality requires a little more care owing to the presence of the incomplete convolution term *q*^{′*}_{m−1}. For the remaining terms, we proceed as we did with the previous inequalities,As this is an upper bound on , we must replace the term *βJ*^{m}*C*_{m} with *βJ*^{m−1}*C*_{m−1} as follows:Using this replacement and the bounds on the numbers , we obtain the upper bound5.26It remains to deal with *q*^{′*}_{m−1}. To do this, we first write *q*^{′*}_{m−1} as follows:5.27The non-convolution terms then give5.28as *C*_{m−3}/*C*_{m−1}≤*C*_{2}/*C*_{4}=1/7, if *m*≥5. The remaining term is treated in a completely different manner,5.29

Thus, adding equations (5.26), (5.28) and (5.29), we set5.30Thus, taking *J*_{2} such that *Q**(*J*_{2}), *R**(*J*_{2}), *S**(*J*_{2})≤1 and , we have shown that the induction hypothesis (5.21) implies5.31so that in fact equation (5.31) holds for every integer *m*. ■

### (e) Proof of theorem 5.2: solution of the eigenvalue problem

### Proof.

The weak form of the master system is5.32

Using that in *P* and *ϵ*_{η}=1 in , and multiplying by , gives the equivalent systemin whichThis form can be expanded in powers of *η*,5.33in which the *a*_{m} are real formsDefine the partial sumsFor *η*<*R*, the sequence converges to a number and the sequence converges in to a function *h*_{η}; thus,Therefore, , *j*=1,…,4, has a convergent series representation in powers of *η*, in which the *m*th coefficient is related to the coefficients *ξ*_{ℓ} and *h*_{ℓ} byFrom these, one obtains the *m*th coefficient of (equation (5.33)), which, by means of the relations *h*_{m}=*h*_{0}*i*^{m}*ψ*_{m}, and , is seen to be equal to the −*i*^{m}*h*_{0} times the right-hand side of equation (3.8). All these coefficients are therefore equal to zero, and we conclude that . This proves that the function *h*_{η}, together with the frequency , solves the weak form (5.32) of the master system. ■

## 6. Effective properties, error bounds and the dispersion relation

In this section, we start by identifying an effective property directly from the dispersion relation. We then discuss the relation between effective properties and quasi-static properties. Next, we provide explicit error bounds for finite-term approximations to the first branch of the dispersion relation for non-zero values of *η*. The error bounds show that numerical computation of the first two terms of the power series delivers an accurate and inexpensive numerical method for calculating dispersion relations for sub-wavelength plasmonic crystals.

### (a) The effective index of refraction—quasi-static properties and homogenization

The identification of an effective index of refraction valid for *η*>0 follows directly from the dispersion relation given by the series for . Indeed, the effective refractive index is defined by expressing the dispersion relation as6.1and it then follows from the expansion for that the effective refractive index has the convergent power series expansion6.2

We now discuss the relationship between the effective index of refraction and the quasi-static effective properties seen in the limit with *k* fixed. The effective refraction index *n*_{eff} can be rewritten in the equivalent form by the equation . By setting *v*=*h*_{η} in the weak form of the master system (5.32), it is easily seen that for all *η* within the radius of convergence, so that for those values of *η*. Following Pendry *et al.* (1999), see also Kohn & Shipman (2008), we define the effective permeability by6.3and we then define *ϵ*_{eff} through the equation6.4The quasi-static effective properties are recovered by passing to the limitsA simple computation shows that *μ*_{qs}=〈*ψ*_{0}〉_{Q}>0 (see appendix, electronic supplementary material). Hence, we have that *μ*_{eff}>0 for *η* in a neighbourhood of the origin, so that *ϵ*_{eff}>0 for these values of *η*, as for all *η* in the radius of convergence. Thus, one has a solid basis on which to assert that plasmonic crystals function as materials of positive index of refraction in which both the effective permittivity and permeability are positive.

For circular inclusions, we have used the program COMSOL to compute that 〈*ψ*_{0}〉_{Q}≈0.98, so that only a mild effective magnetic permeability arises.

Having established that is the solution to the unit cell problem, we can undo the change of variable ** y**=

*k*

**/**

*x**η*to see that the function6.5where is the

*Q*-periodic extension of

*h*

_{η}to all of , is a solution of6.6for every

*η*in the radius of convergence.

We directly investigate the quasistatic limit using the power series (6.5). Here, we wish to describe the average field as . To do this, we introduce the three-dimensional period cell for the crystal [0,*d*]^{3}. The base of the cell in the *x*_{1}*x*_{2}-plane is denoted by *Q*_{d}=[0,*d*]^{2} and is the period of the crystal in the plane transverse to the rods. We apply the definition of *B*_{eff} and *H*_{eff} given in Pendry *et al.* (1999), which in our context is6.7and6.8Taking limits for *k* fixed and in equation (6.5) givesin which . These are the same average fields that would be seen in a quasi-static magnetically active effective medium with index of refraction *n*_{qs} and *μ*_{qs} that supports the plane waveswhere *μ*_{qs}=〈*ψ*_{0}〉_{Q}. It is evident that these fields are solutions of the homogenized equation . This quasi-static interpretation provides further motivation for the definition of *n*_{eff} for non-zero *η* given by equation (6.2).

Now, we apply the definition of effective permeability *μ*_{eff} given in Pendry *et al.* (1999), together with *n*_{eff} to define an effective permeability *ϵ*_{eff} for *η*>0. The relationships between the effective properties and quasi-static effective properties are used to show that plasmonic crystals function as meta-materials of positive index of refraction in which both the effective permittivity and permeability are positive for *η*>0.

### (b) Absolute error bounds

The Catalan bound provides simple estimates on the size of the tails for the series and *h*_{η},We have established convergence for *η*≤1/4*J*, so that we may write *η*=*α*/4*J*, 0≤*α*≤1. Then, using that , *C*_{2m}≤4^{2m} and 4*Jη*=*α*, we have6.9Similarly, for *h*_{η} we have that6.10

### (c) Relative error bounds

In this section, we use the absolute error bound (6.10) with *m*_{0}=1 to obtain a relative error bound for the particular case of a circular inclusion of radius *r*=0.45 (Shvets & Urzhumov 2004). The first term approximation to *h*_{η} is6.11For a circular inclusion of radius *r*=0.45, we have *J*≤85, *β*≤0.79, ∥*ψ*_{0}∥=0.97 and ∥*ψ*_{1}∥=0.02, where ∥⋅∥=∥⋅∥_{H1(Q)}. Thus, using bound (6.10), the relative error *R*_{1,h} is bounded byso that, for *α*≤0.2, the relative error is less than 8.2 per cent. The graphs of *ψ*_{0} and *ψ*_{1} can be found in the appendix in the electronic supplementary material. The first term approximation to is6.12In the appendix in the electronic supplementary material, we indicate how the tensors may be computed. For an inclusion of radius *r*=0.45, we have and . Thus, using bound (6.9), the relative error *R*_{1,ξ} is bounded byso that, for *α*≤0.3, the relative error is less than 2 per cent (figures 2 and 3).

## Acknowledgements

R.L. was supported by NSF grant DMS-0807265 and AFOSR grant FA9550-05-0008, and S.S. was supported by NSF grant DMS-0807325. The PhD work of S.F. was also supported by these grants. This work was inspired by the IMA ‘Hot Topics’ Workshop on Negative Index Materials in October 2006.

## Footnotes

- Received October 24, 2009.
- Accepted January 11, 2010.

- © 2010 The Royal Society