## Abstract

We obtain band diagrams for a two-dimensional periodic structure consisting of an infinite square array of infinitely thin concentric circles (split rings) with narrow gaps. Our approach exploits the narrowness of the gaps and yields algebraic equations relating the frequency to the Bloch wavenumber and geometric properties of the array. Further asymptotic analysis indicates that the gravest mode has a frequency that scales in an inverse logarithmic fashion with the size of the gap and that exhibits anomalous dispersion. Near the origin of the Brillouin zone this ‘acoustic’ mode is dispersionless. Numerical solution of the eigenvalue problem in the single-gap case confirms these conclusions. The two lowest modes of the split ring can be interpreted as a splitting of the gravest propagating Rayleigh mode.

## 1. Introduction

The study of the propagation of waves through periodic structures has a long and distinguished history going back to Rayleigh (1892) and including classic books such as Brillouin (1953). Interest in the field continues today, given that it underpins modern technological applications that exploit the frequency-dependent filtering and bandgap properties of periodic structures. This rapidly growing area is sometimes referred to as photonics and the devices called photonic crystals. The photonics literature is now vast; we limit ourselves to citing Sakoda (2005) as a helpful introduction and Maynard (2008) for a very elementary description of wave propagation in arrays of scatterers.

Pendry *et al.* (1999) first suggested how to construct ‘left-handed’ media (i.e. with a negative refractive index) by using a split-ring resonator structure. This was the beginning of efforts to design tailored materials exhibiting properties not found in naturally occurring materials. Subsequently, Smith *et al.* (2000) built a split-ring resonator, demonstrating that such materials could in fact be constructed. Subsequently, Movchan & Guenneau (2004) analysed a double ‘C’ resonator design analogous to the split-ring resonator for a two-dimensional cylindrical structure in an elastic medium. Guenneau *et al.* (2007) further examined this structure using a discrete approximation and predicted the lowest stop and dispersion bands.

We study a very similar geometry in which each resonator is made up of two perfectly insulating split rings, each with a hole in it, as illustrated in figure 1. In the single-gap case, we remove the inner ring. We consider just the scalar wave equation. The rings are in a square lattice with periodicity *d*. The gaps are centred at *θ*_{1} and *θ*_{2} and hence in general are not aligned. We combine the asymptotic technique of Davis & Llewellyn Smith (2007), hereafter referred to as DLS, with the method originally introduced by Rayleigh, which is clearly explained in Movchan *et al.* (2002). In the problem studied by Rayleigh, the dispersion relation for plane waves propagating through an infinite array of scatterers is obtained as a function of the Bloch wavenumber **k**_{0}.

The approach of the present paper proceeds by obtaining and matching potentials inside the inner ring, between the rings, and outside the outer ring. We hence calculate the dispersion relation when the gaps are narrow. Our result is the dispersion relation for frequency as a function of the Bloch wavenumber **k**_{0} and of the geometric parameters of the problem. The Rayleigh modes are present as resonances of the system when the gap size vanishes. Guizal & Felbacq (2002) examine scattering by a cylinder with gaps, but the problems are very different since in that work there is a single cylinder, the medium is infinite rather than periodic and the results are purely numerical rather than asymptotic. Our approach is not limited to the specific shape considered (see DLS). More complicated shapes require solving numerically for the interior Green’s function. The present geometry allows explicit results to be given in terms of Bessel functions, which are lengthy but simple to use.

We present the solution in the two interior regions inside the inner ring and between the two rings in §2. The solution in the array is obtained in §3. The two are matched to obtain an eigenvalue problem in §4. Limiting cases of the dispersion relation for small *k* are found in §5, while solutions to the full eigenvalue problem in the single-gap case are obtained in §6. We conclude by summarizing in §7.

## 2. Mathematical model of the cell

We first suppose that, in terms of polar coordinates (*r*,*θ*), the potential *ϕ*(*r*,*θ*) satisfies the Helmholtz equation
2.1
inside and outside the cavities 0≤*r*<*R*_{1}, *R*_{1}<*r*<*R*_{2} that are connected by gaps at *r*=*R*_{1},|*θ*−*θ*_{1}|<*δ*_{1} and *r*=*R*_{2},|*θ*−*θ*_{2}|<*δ*_{2}. We assume that the concentric boundaries are hard and, following DLS, we incorporate the square root edge singularities by writing the boundary conditions at *r*=*R*_{1},*R*_{2} in the form
2.2
and
2.3
where *T*_{N} denotes the Chebyshev polynomial of order *N* and the *a*_{N1} and *a*_{N2} are coefficients that will be obtained as part of the solution.

In the central cavity, *r*<*R*_{1}, we need a Green’s function, *G*_{1}(*r*,*θ*;*r*′,*θ*′) that satisfies equation (2.1) in the primed coordinates and the hard condition ∂*G*_{1}/∂*r*′=0 at *r*′=*R*_{1}. This is given by
2.4
where *ϵ*_{n} denotes Neumann’s symbol. We apply Green’s Theorem to *ϕ*(*r*′,*θ*′) and *G*_{1}(*r*,*θ*;*r*′,*θ*′) in *r*,*r*′<*R*_{1} and obtain
2.5
in which, according to equation (2.4),
2.6
The first expression for *G*_{1} enables equation (2.5) to be obtained directly by inserting a delta-function in (∂*G*_{1}/∂*r*)_{r=R1} (see DLS). The second expression for *G*_{1} displays explicitly the logarithmic singularity.

In the annular cavity, *R*_{1}<*r*<*R*_{2}, we need a Green’s function, *G*_{2}(*r*,*θ*;*r*′,*θ*′) that satisfies equation (2.1) in the primed coordinates and the hard conditions ∂*G*_{2}/∂*r*′=0 at *r*′=*R*_{1},*R*_{2}. We note that its structure must be given by
2.7
where *A*_{n} and *B*_{n} are constants to be determined. The analogues of equation (2.6) are then
2.8
and
2.9
As in equation (2.6), the logarithmic singularity in equation (2.8) is displayed by writing
2.10
We now apply Green’s Theorem to *ϕ*(*r*′,*θ*′) and *G*_{2}(*r*,*θ*;*r*′,*θ*′) in *R*_{1}<*r*′<*R*_{2} and obtain
2.11
Then substitution of equations (2.2) and (2.3) into equations (2.5) and (2.11) gives the following ‘gap’ values of the potential, valid for |*s*|<1:
2.12
2.13
and2.14
Inspection of equations (2.6), (2.8) and (2.9) shows that *πG*_{1}(*R*_{1},*θ*_{1}+*δ*_{1}*s*;*R*_{1},*θ*_{1}+*δ*_{1}*s*′), *πG*_{2}(*R*_{1},*θ*_{1}+*δ*_{1}*s*;*R*_{1},*θ*_{1}+*δ*_{1}*s*′) and *πG*_{2}(*R*_{2},*θ*_{2}+*δ*_{2}*s*;*R*_{2},*θ*_{2}+*δ*_{2}*s*′) have a singular term
2.15
Following DLS, we deduce approximations to the gap potentials in |*s*|<1 that are given, with errors of order *δ*_{1}, *δ*_{2}, by
2.16
2.17
and
2.18
Before matching the ‘gap’ potentials, we must construct the field in *r*>*R*_{2} which has a doubly periodic, square array of the annular boundaries described above.

## 3. Mathematical model of the doubly periodic array

The resonator consists of square cells of side *d*(>2*R*_{2}), centred at
3.1
referred to unit vectors directed along Cartesian axes. Each cell has the structure detailed in the previous section and the potential *ϕ* satisfies equation (2.1) everywhere. The quasi-periodicity of the problem is manifested through the Bloch wavenumber **k**_{0} by imposition of the condition
3.2
For a given **k**_{0}, eigenvalues of *k*^{2} determine the frequencies of propagating modes and singular values of *k*^{2} determine the resonant frequencies. Evidently, equation (3.2) suffices to determine *ϕ* from its values in the cell centred at **R**_{0}=(0,0). In the region *r*>*R*_{2} of this cell, the fundamental quasi-periodic Green’s function is given by
3.3
which satisfies the Bloch conditions
3.4
and
3.5
(One can also define a Green’s function using Hankel functions rather than Bessel functions of the second kind.) With the notation **r**−**r**′=*ξe*^{iγ}, **R**_{p}=*R*_{p}*e*^{iΘp}(*p*≠0), we follow Movchan *et al.* (2002) in writing equation (3.3) as
3.6
where
3.7

We now apply Green’s Theorem to *ϕ*(*r*′,*θ*′) and *G*(*r*,*θ*;*r*′,*θ*′) in the region *r*′>*R*_{2} of the cell centred at **R**_{0}=(0,0). The quasi-periodicity (3.2) and conjugate quasi-periodicity (3.4) conditions, as functions of **r**′, enable us to obtain
3.8
The boundary values of *ϕ* are now introduced and expressed as Fourier series by writing
3.9
The first of these is essentially a restatement of equation (2.3), while the second one introduces a new set of unknown coefficients *D*_{n}. We need a double Fourier series for *G* in order to deduce equations for the coefficients *D*_{n} from equation (3.8). Applying the Graf summation formula (e.g. Abramowitz & Stegun 1969, eqn 9.1.79) to equation (3.6) yields
3.10
Substitution of equations (3.9) and (3.10) into equation (3.8) now gives
3.11
which satisfies equation (3.9) provided
3.12
for . Since *R*_{−p}=*R*_{p},*Θ*_{−p}=*Θ*_{p}+*π*, we deduce from equation (3.7) that
3.13
Hence the matrix of coefficients in equation (3.12) is Hermitian. Thus, the resonant frequencies are real and coincide with the propagating frequencies in the absence of the outer gap (Rayleigh 1892), that is, when *δ*_{2}=0. Further, the formula (3.10) is real-valued when *θ*=*θ*′ because the double sum is then a Hermitian quadratic form.

From equations (2.3) and (3.9),
3.14
which suffices for the asymptotic estimate
3.15
in equations (3.9) and (3.12). We then deduce from equation (3.9) that the ‘gap’ value of the potential is approximated for |*s*|<1, with errors of order *δ*_{2}, by
3.16
To show that this estimate is real, we multiply the conjugate of equation (3.12) by *D*_{n}*J*_{n}′(*kR*_{2}) and sum over *n*. Using equation (3.13), this gives, without consideration of convergence,
3.17
which is real-valued because the double sum is a Hermitian quadratic form.

## 4. The eigenvalue problem

Admissible values of *k*^{2}, and hence the frequency at given wavenumber **k**_{0}, are determined by matching values of the potential in the two gaps. At *r*=*R*_{1}, equations (2.16) and (2.17) give
4.1
while, at *r*=*R*_{2}, equations (2.18) and (3.16) give
4.2
with errors of order *δ*_{1} and *δ*_{2}.

Since the interior gap orientation *θ*_{1} appears only in |*θ*_{1}−*θ*_{2}|, we have verification that in the annular region only the separation of the gaps is significant in this asymptotic framework. In the periodic array, we note a symmetry with respect to direction reversal of either **k**_{0} or **R**_{p} and prefer to consider the latter. We replace *n* by −*n* in the conjugate of equation (3.12) and use equation (3.13) to obtain
4.3
which allows us to deduce from equation (3.12) that the coefficients are related by
4.4
When equation (4.4) is used to evaluate equation (3.17) at *π*+*θ*_{2}, we find that
4.5
with the last identity a consequence of the real-valued sum. Thus the last term in equation (4.2) is unchanged when *θ*_{2} is replaced by *π*+*θ*_{2}. So we need consider only the range −*π*/2<*θ*_{2}≤*π*/2 in the eigenvalue problem. But this reflective property of *θ*_{2} does not change the sufficient range |*θ*_{1}−*θ*_{2}|≤*π* of the gap offset.

A simpler calculation suffices if we delete the annular region and have a single gap by setting *R*_{1}=*R*_{2}=*R*, *θ*_{1}=*θ*_{2}, *δ*_{1}=*δ*_{2}=*δ*. Then equations (3.13) and (3.15) are unchanged except for *R*_{2}=*R* and matching of equations (2.16) and (3.16), with *a*_{01}=*a*_{02}=*a*_{0}, gives
4.6
with errors of order *δ*. The symmetry property (4.4) remains applicable.

Alternative expressions for the lattice sum are derived in Movchan *et al.* (2002). For later estimates, we need eqn (3.101) of Movchan *et al.* (2002):
4.7
which is obtained by expressing the function in wavenumber space. The reciprocal lattice is defined by **Q**_{h}=**k**_{0}+**K**_{h}=*Q*_{h}*e*^{iΓh} and **K**_{h}=2*π*(*h*_{1}**e**_{x}+*h*_{2}**e**_{y})/*d*, with .

## 5. Small eigenvalue asymptotics

When we seek an estimate of the gravest eigenvalue for the single-gap case by inserting small-*k* expansions into equation (4.6), we deduce from equation (3.12) that the second term gives an *O*[1/(*kR*)^{2}] contribution, as in the cases discussed by DLS. This leads to the estimate
5.1
with errors of order *δ* and valid for |**k**_{0}|≫*k*. Large values of the coefficients in equation (4.6), and hence eigenvalues, also occur when *k* is close to either a Rayleigh mode or an interior mode in the absence of the gap. This is analogous to the case of an outer concentric circular boundary, with equations (3.12) and (4.6) together displaying the features observed in equation (3.17) of DLS.

In seeking similar approximations to the two-gap matching formulas, we first note that the coefficients {*D*_{n}} in equation (4.2) are still determined by equation (3.12), because the array geometry is a common property. In particular, *D*_{0} gives an *O*[1/(*kR*_{2})^{2}] contribution, which must be taken into account in a small-*k* estimate of the gravest eigenvalue. Arranged for the dominant terms to cancel in their determinant, the generic form of equations (4.1) and (4.2) is
5.2
and
5.3
where *A*=(*R*_{2}/*R*_{1})^{2}>1 is the area ratio of the split ring. Evidently *B*_{1} must include −(*kR*_{1})^{−2}, which therefore is at most logarithmically large. On setting the determinant to zero, the surviving leading terms in the determinant yield the identity *B*_{1}+*B*_{2}=2*B*_{3} which enables us to deduce, after much algebra, the small-*k* estimate
5.4
with errors of order *δ*_{1} and *δ*_{2}. Again, eigenvalues also occur when *k* is close to either a Rayleigh mode or an interior mode in the absence of the gaps, with the extra equation accounting for the annular modes.

### (a) *k*_{0}*d*=*O*(1)

*k*

The validity of equations (5.1) and (5.4) depends on the existence of the stated limit. We rearrange equation (3.12) to obtain
5.5
which defines *g*_{n} and in which *β*_{n} is estimated by equation (3.15). Since equation (3.7) implies that if *l*≠0, we deduce that, if *n*≠0 in equation (5.5), *g*_{n}=*O*(1) and only coefficients with *m*≠0 and *mn*<0 are *O*(1). Equation (4.7) shows that
5.6
and the equations (5.5) reduce to
5.7
Since *D*_{0} is absent from the *n*≠0 equations, the middle equation is an *O*(1) estimate for *D*_{0} that depends on the other, already computed, limit values of the coefficients.

### (b) *k*_{0}=|*k*_{0}|=*O*(*k*)

*k*

This is the acoustic limit, and we can obtain the asymptotic behaviour by truncating the full system (3.12) to three terms:
5.8
The required estimates are (Movchan *et al.* 2002, eqns 3.132–134, 3.155–156)
5.9
Solving equation (5.8) shows, as expected, that
5.10
and hence the estimate (5.1) is inapplicable here. Instead, to leading order, the eigenvalue equation (4.6) for the single-gap case gives from equations (5.8) and (5.9)
5.11
where *f*=*π*(*R*/*d*)^{2} is the area fraction inside the outer ring. This estimate has errors of order *δ* and valid for |**k**_{0}|=*O*(*k*)≪*d*^{−1}. As *k*→0, both sides of equation (5.11) vanish, leading to the acoustic dispersion relation
5.12
The matrix in equation (5.8) is singular when *k*=(1+*f*)^{−1/2}*k*_{0}, which is the Rayleigh dispersion relation (Movchan *et al.* 2002, eqn 3.158) corresponding to modes propagating through a periodic array of circular cylinders.

## 6. Numerical approach and results

We limit ourselves to the single-gap case in the numerical calculations. Our task is to solve equation (4.6), that is find the zeroes of a function of the frequency *k*, the wavenumber **k**_{0}, or rather its angle *θ* and the geometric parameters of the array: *R*, *δ* and *D*. We view equation (4.6) as an equation for *k* and solve for *k* on the boundary of the first irreducible Brillouin zone KGM. To compute the sum in equation (4.6), we first need to solve for the *D*_{n} in equation (3.12). We truncate the infinite series at *N*=8. This is more than enough, given the spectral accuracy of the series. Tests showed negligible differences between *N*=4 and 8. The lattice sums are computed using the approach of McPhedran *et al.* (2000) and Yasumoto & Yoshitomi (1999). Tests matched the digits in these references. While our interest was focused on the lowest three modes, the possibility of crossing of higher modes made the computation of some of the dispersion curves delicate. To deal with this issue, we used the continuation code PITCON66 (a descendant of the algorithm in Rheinboldt & Burkardt 1983) to track the curves. Rather than computing the sum in equation (5.1) at *k*=0 as in equation (5.6), we use small values of *k* and use a parabolic fit to obtain the values at *k*=0.

As default values we take *d*=1, *R*=0.35, *δ*=0.05 and *θ*=0. This produces the three lowest modes in the band diagram of figure 2*a*. We see that the gravest mode (mode 1) has the typical V-shape of the acoustic mode and then weak variation along the rest of the boundary of the Brillouin zone. Hence the group velocity is low, corresponding to anomalous dispersion. The next two modes are typical of such period arrays, and in fact the whole diagram is very similar to figure 2 of Guenneau *et al.* (2007). There are band gaps between all three modes.

The other panels in figure 2 correspond to varying one parameter of the default values. In figure 2*b*, *θ*=*π*/3; mode 3 has lower group velocity along most of the boundary of the Brillouin zone, but there are no qualitative changes. In figure 2*c*, *δ*=0.001 while in figure 2*d*, *δ*=0.2. The small-*k*_{0} asymptotic limit decreases as *δ* decreases and the group velocity is even smaller, while the opposite happens as *δ* increases. Mode 2 shows smaller changes than mode 3. In figure 2*e* *R* is reduced to 0.2. Mode 2 now has a frequency larger than 6 and is no longer visible. The small-*k*_{0} limit is larger since it is inversely proportional to *R*. In figure 2*f* *R*=0.48, which has the outer rings almost touching. The small-*k*_{0} limit is lower and shows more structure. Both modes 2 and 3 show weak dependence on *k*_{0} with small group velocities.

The features seen in figure 2 can be interpreted as a splitting of the lowest propagating Rayleigh mode. In the split ring case, the lowest mode has a universal shape, with essentially the acoustic mode in the large-wavelength, low-frequency limit, asymptoting to an anomalous trapped mode along the rest of the boundary of the Brillouin zone. Mode 1 is then a slightly distorted version of the lowest Rayleigh mode with a smooth parabolic form near *G*. This is shown in figure 3 for the default parameter values of figure 2*a*. The different curves show the split ring dispersion curves, the Rayleigh dispersion curves and the asymptotic limits. As *δ* decreases, the gravest split ring mode becomes flatter with lower frequency, and mode 2 resembles the gravest Rayleigh mode even more closely. Mode 2 and Rayleigh mode 1 do not seem to have any particular resemblance.

## 7. Summary

We have computed the band diagram for a periodic structure of split ring resonators. An analysis valid when the gaps are narrow gives an algebraic dispersion relation. Asymptotic estimates are obtained for low frequencies, showing anomalous dispersion. For the single gap we find the acoustic mode with *k*∼[(1−*f*)/(1+*f*)]^{1/2}*k*_{0}. Numerical solutions for the single gap show that these asymptotic predictions are robust. Band gaps are also found between the lowest three modes. The lowest two modes can be interpreted as a splitting of the gravest propagating Rayleigh mode.

The only extra complication in the numerical solution for the double-ring case is the solution of the two linear equations (4.1) and (4.2) rather than the single equation (4.6). The asymptotic behaviour for the double-ring case is obtained and is no more difficult to implement than the single-ring case, since both require finding the asymptotic behaviour of the coefficients *D*_{n}.

This work shows how asymptotic methods can be used to demonstrate how band diagrams for periodic structures are changed by the introduction of narrow gaps in the geometry.

## Acknowledgements

We are grateful to Richard Craster for helpful comments.

## Footnotes

- Received January 29, 2010.
- Accepted April 16, 2010.

- © 2010 The Royal Society