Spectral super-resolution in metamaterial composites

We investigate the optical properties of periodic composites containing metamaterial inclusions in a normal material matrix. We consider the case where these inclusions have sharp corners, and following Hetherington and Thorpe, use analytic results to argue that it is then possible to deduce the shape of the corner (its included angle) by measurements of the absorptance of such composites when the scale size of the inclusions and period cell is much finer than the wavelength. These analytic arguments are supported by highly accurate numerical results for the effective permittivity function of such composites as a function of the permittivity ratio of inclusions to matrix. The results show that this function has a continuous spectral component with limits independent of the area fraction of inclusions, and with the same limits for both square and staggered square arrays.


Introduction
This paper links themes evoked in two classic papers, one in mathematics [1] and the other in physics [2]. The first of these poses the question as to whether the spectral content of the radiation from a body can reveal its shape. The second shows that the use of negatively refracting metamaterials in a plane slab can lead to a super-resolving perfect lens, also known as a superlens. We will consider a two-dimensional composite material, composed of polygonal inclusions made of a metamaterial (by which we mean an artificial material with a dielectric constant that has a negative real part and a very small imaginary part) and placed in a positive dielectric matrix material. We will show that, in the spirit of Pendry, the metamaterial makes possible the resolution of an important structural feature of the inclusions, irrespective of how much smaller than the wavelength they are. We will also show that, in the spirit of Kac, this feature relates to the shape of the inclusion, being in fact the corner angle of the polygon, and that it is deduced from spectral measurements on the composite. The fact that a spectral feature could be determined by corner shape, independent of (say) the area fraction of inclusions, was first suggested by Hetherington and Thorpe [3] on the basis of an elegant argument and numerical evidence for dilute composites.
We will base our demonstration, firstly, on analytic results relating to the spectrum of the effective dielectric permittivity function eff of the composite material and, secondly, on remarkably accurate numerical results for this spectrum obtained using a new technique. Note here that we are using the word 'spectrum' in two related, but slightly different senses. In the previous paragraph, its usage meant that the absorption of electromagnetic waves by the composite was being determined as a function of wavelength, ranging over an appropriately wide band. In the first sentence of this paragraph, we referred to the distribution of singularities of the function eff , giving the effective permittivity of a composite having a specified geometry as a function of the ratio σ = 1 / 2 of the permittivities of inclusions and the matrix. The relation between these usages is that, as the wavelength varies so does the ratio σ , so that measurements of (say) optical absorption by a composite over a suitable wavelength interval can reveal details of the singularities of the function eff .
The numerical results for the singularity spectrum of eff reveal that it has a continuous part that runs between the upper and lower limits of σ , which do not vary at all with the area fraction of the inclusions. It is complemented by a discrete spectrum of poles, which does 3 evolve with area fraction. This evolution is in fact necessary, since the continuous spectrum for touching square inclusions in a checkerboard arrangement occupies the entire negative real axis of σ , but for non-touching square inclusions it is confined to the interval −3 σ −1/3. The animations we give show how this transformation is achieved: the discrete spectrum becomes more and more dense as the touching configuration is approached, to supply the required spectral extension, as anticipated by one of the authors [4].
The results we give here are interesting for the insights they give into the connection between metamaterials and super-resolution. They are also important in furthering our understanding of the connection between inclusion shape, geometrical arrangement and spectral properties of the effective permittivity function. This connection helps in the design of structures having enhanced absorption over a wide wavelength range for applications in photothermal or photovoltaic captors [5][6][7][8], or offering strongly enhanced local fields for applications such as sensing or nonlinear optical elements.
The topic of super-resolution using metamaterials is an active one. Two interesting methods for achieving super-resolution are hyperlenses [9,10,12] and perfect imaging without negative refraction [11,13]. However, the proposal we investigate here for morphological super-resolution using the optical absorption spectrum has, the authors believe, not been discussed in the literature previously, beyond the work of Meixner [14] and Hetherington and Thorpe [3]. By morphological super-resolution, we mean that the technique is primarily sensitive to the shape of the particles in the unit cell, but offers little information about their geometrical distribution. This is quite different from inverse problem methods [15][16][17][18] that are capable of inferring information about the geometrical distribution (volume fractions, geometric parameters associated with correlation functions, etc), but offer little information about the shape of the particles.
We give a brief overview in section 2 of some of the important properties of the function eff (σ ), including analytic results relating to the continuous part of the spectrum and some numerical investigations of both the discrete and continuous parts of the spectrum. In section 3, we describe the method that enables accurate calculation of the spectrum of eff and give numerical results illustrating its convergence, even on the negative real axis of σ . In section 4, we discuss the animations that are given in the supplementary material (available from stacks.iop.org/NJP/13/115005/mmedia) of this paper and the physical consequences of the behaviour they show. We present a discussion and the concluding remarks in section 5.

Properties of the effective permittivity function for composites
We will now give a concise review of what is known about the properties of the effective dielectric permittivity function eff for composites made of two materials with dielectric permittivities 1 and 2 , with the former corresponding to a phase consisting of isolated particles and the latter to a continuous matrix phase. We assume that the geometry has cubic symmetry (in three dimensions) or square symmetry (in two dimensions) so that eff is scalar valued, i.e. the effective dielectric tensor equals eff I . This review builds on that given by Perrins and McPhedran [19].
The calculation of eff for a given geometry is homogeneous of degree 1 in the variables 1 and 2 , so we can rescale to make eff a function of a single complex variable, the permittivity ratio σ = 1 / 2 : 4 From now on, we let eff be an abbreviated notation for the effective dielectric permittivity function eff ( 1 , 2 ) and we let eff (σ ) denote its scaled counterpart, the effective relative dielectric permittivity function, as defined via (1).
To determine eff (σ ) for a given geometry, we need to solve an electrostatic transport problem repeatedly for various σ . More precisely, we need to solve Laplace's equation for the potential V on a periodic domain with a periodic electric field −∇V (x, y) having a prescribed average value, and boundary conditions of continuity of V and its normal flux ∂ V /∂n at interfaces between materials. The theory of the function eff (σ ) becomes particularly elegant when we deal with two-dimensional problems, in which V = V (x, y) becomes a function in the plane. We may then apply the apparatus of complex-variable theory to the calculation of V and thus to eff (σ ). For the case of a doubly periodic array of inclusions C with unit cell U and square symmetry, the effective permittivity may be defined as [20] where the integral in the numerator includes contributions E 1 from the inclusion region and E 2 from the matrix region. Except for occasional comments on effective permittivity in three dimensions, we will concentrate on two dimensions, which corresponds to arrays of cylinders of arbitrary cross-section, with the average field aligned in the (x, y) plane. The area fractions of the two phases will be denoted by p 1 and p 2 .
Since the geometry has square symmetry, Keller's theorem [21] gives This equation then pairs zeros of eff (σ ) at values σ 0 with poles at values σ p = 1/σ 0 . Of course, from (3), zeros of eff (σ ) require that the contributions E 1 and E 2 add up to zero. Since, with 2 = 1, E 2 is real and positive, this means that E 1 must be real and negative, and so 2 = σ must be real and negative at any zero of eff (σ ) and thus at any pole as well. Bergman [20] proved this property for both two-and three-dimensional composites. Bergman also recognized that even though the transport problem when the average electric field −∇V (x, y) is prescribed does not have a solution at a pole, it should have a solution at a pole when instead the average displacement field − ∇V (x, y) is prescribed. Milton [22] proved that for the value σ = −1 the electrostatic problem of an array of circular cylinders, with either a prescribed value of the average electric field −∇V (x, y) or a prescribed average value of the displacement field − ∇V (x, y), does not have a solution; compare with the discussion of (A.1) below. This suggests that σ = −1 is either an essential singularity or lies on a branch-cut of the function eff (σ ) for arrays of circular cylinders. The fact that branch-cuts cannot be in the upper or lower half-planes, but must lie exactly on the negative real axis of σ was proved by one of the authors [23] using the relationship between composite materials and resistor networks. A rigorous justification of the spectral representation for eff (σ ) was given by Golden and Papanicolaou [24]. Assuming that a certain moment is finite [25], the function eff (σ ) has the representation where a 1 and the spectral measure dµ(τ ) are non-negative. The support of dµ(τ ) is the spectrum. The spectral measure can be recovered from the values the imaginary part of eff (σ ) takes near the negative real σ -axis, since the integral of any smooth test function g(τ ) with respect to the measure dµ(τ ) is given by The discrete spectrum of eff (σ ) is readily exhibited numerically. This has been done for arrays of spheres by Bergman [26] and for arrays of circular cylinders by McPhedran and McKenzie [27], with both studies showing that the poles and zeros of eff (σ ) converge to an essential singularity at σ = −1.
We now focus on what can be said about the continuous spectrum of eff (σ ). One simple geometry for which a result is immediately apparent is the checkerboard, for which Dykhne [28] obtained from Keller's theorem (3) the exact result This then exhibits a branch-cut along the entire negative real axis of σ . An exact result can also be obtained for the polarizability α of a pair of touching cylinders [5]. This can be employed to estimate the effective permittivity of dilute arrays or random arrangements of pairs of touching cylinders, using the Clausius-Mossotti equation [4], and a similar integral representation exists for the polarizability. Using an inversion of coordinates about the contact point, the touching cylinders may be transformed into a slab of matrix material with permittivity 2 = 1 surrounded by two half-planes filled with material with permittivity 1 = σ . Introducing the parameter it is easy to show using the method of images that the polarizability for a pair of touching cylinders of radius a for the case of the applied field parallel to the line connecting cylinder centres is We see that the series in (8) converges provided |λ| < 1, i.e. for real σ , σ > 0. However, we can obtain a meaningful result even when this is not the case by the technique of analytic continuation, since the series in (8) is a known transcendental function, called the dilogarithm, and denoted by Li 2 . Thus, we can replace (8) by The properties of the dilogarithm function are that it has a branch-cut running from λ = 1 to λ = ∞, across which the discontinuity in the imaginary part of Li 2 (λ) is 2π log [ {λ}]. The branch-cut in the plane of relative permittivity runs from σ = −∞ to σ = −1. If the direction of the applied field is perpendicular to the line of centres, the branch-cut runs from σ = −1 to σ = 0. We next consider arrays of square inclusions, for which we have already mentioned the Dykhne result (6). A generalization of this for an array in which the square unit cell was divided into four equal squares with dielectric permittivities 1 , 2 , 3 , and 4 (in clockwise order) was conjectured by Mortola and Steffé [29], and in one direction the effective dielectric constant was proposed to be eff = which reduces to the arithmetic average eff = ( 1 + 3 )/2 when 1 = 2 and 3 = 4 , corresponding to the effective dielectric constant of a laminate of equal portions of phases 1 and 3, in the direction parallel to the layers. This conjecture was proved independently by Milton [30] and Craster and Obnosov [31]. We will consider a particular sub-case of this result, due to Obnosov [32]: an array of square cylinders with area fractions p 1 = 0.25 and p 2 = 0.75, for which (10) gives This formula yields a spectrum consisting solely of a branch-cut running from σ = −3 to σ = −1/3. We compare the result given by this formula for the real part of eff (σ ) with the result of a numerical mode-matching procedure in figure 1. This comparison reveals the difficulty of evincing details of the spectrum using numerical methods: the mode matching technique approximates the branch-cut by a discrete set of poles, which becomes more dense as the number of modes increases. However, it is difficult to distinguish between branch-cuts and sets of poles concentrating around an essential singularity by such methods. Furthermore, the mode-matching method failed to give clear indications of the spectrum for area fractions of cylinders distinct from p 1 = 0.25. The behaviour of fields near corners of inclusions in a matrix of differing dielectric permittivity and magnetic permeability was treated in electromagnetism by Meixner [14]. As well as obtaining the exponents that characterize the singularity of the field components near the edge, Meixner pointed out that exactly the same formulae could be applied in electrostatics 7 and magnetostatics. Hetherington and Thorpe [3] analysed the behaviour of fields near corners in electrostatics and magnetostatics, apparently without knowledge of Meixner's paper. They, however, pointed out the link between field behaviour near corners and the nature of the singularity spectrum in composites containing inclusions with corners. Let us suppose that the electrostatic potential varies with distance r from a corner with an included angle 2ψ as r β and that the permittivity ratio inside the inclusion to that outside is σ . Then [3,14], β is found by solving a transcendental equation: or the same equation with σ replaced by 1/σ . For σ real, the solution β of (12) is either pure real or pure imaginary.
To clarify the behaviour of fields, consider the case of ψ = π/4, corresponding to a 90 • corner. Then (12) gives From this, we see that the right-hand side exceeds one in magnitude for σ lying between −3 and −1/3 and that β is then pure imaginary, corresponding to a solution for the potential that oscillates with r , and the oscillations become ever more rapid as r tends to zero. The electric field is given by the spatial derivative of the potential, and it diverges like 1/r multiplied by the oscillating term as r → 0. The branch-cut region here is that where β is imaginary, i.e. between σ = −3 and σ = −1/3. Hetherington and Thorpe [3] postulate that for an m-sided regular polygon, there will be a branch-cut running between σ = (2 + m)/(2 − m) and σ = (2 − m)/(2 + m). They made this assertion after recognizing that in this interval of σ the surface charge (rather than just the surface charge density) near the corner is infinite, which is unphysical.
Another argument for the position of this branch-cut was put forward by one of the authors [4]. In order that eff have a significant imaginary part when has a very small imaginary part, we see from (2) that the electric field −∇V must be close to losing its square integrability. From the asymptotic form of V near the corner, one finds that this happens when the imaginary part of σ is very small and the real part of σ is between (2 + m)/(2 − m) and (2 − m)/(2 + m). To elucidate this further for the case of a 90 • corner, we solve (13) with β = β + i β and σ = σ + i σ , where β , β , σ and σ are real, −3 < σ < −1/3, and σ and β are both very small and positive. This gives Using polar coordinates (r, θ ) near the corner, the potential V scales as r β as r → 0, and so |∇V | 2 will be close to r 2β −2 |g(θ, σ )| 2 for some function g(θ, σ ). It follows that with 1 = σ and 2 = 1 the imaginary part of E 1 (which is a measure of the power dissipation in the composite) has a contribution near the corner from inside the radius r = r 0 of where the integral over θ is only over those angles in the inclusion. Thus, provided g(θ, σ ) is nonzero, the contribution of the corner to the imaginary part of E 1 remains nonzero even in the limit σ → 0 (and goes to zero when σ < −3 or σ > −1/3, since then σ /β → 0). A corner is not the only geometric feature that can act as a significant energy absorber when the imaginary part of the dielectric constant goes to zero. The centre of a sphere with a dielectric constant λ 1 in the radial direction and dielectric constant λ 2 in the tangential direction acts as an absorber when λ 2 /λ 1 approaches real values less than −1/8: see figure 4 in the paper by Qui and Lukyanchuk [34] (which shows that this energy-absorbing feature extends beyond the quasistatic limit) and see also the related discussion on page 239 of [4] (where there is an error as the limit δ → −1/2 should have been taken rather than the limit δ → 0).

Numerical method
We now describe a numerical method stable and accurate enough to verify the conjecture of Hetherington and Thorpe for the case m = 4 and to show the spectral evolution as a function of area fraction for composites with square inclusions.
Laplace's equation is to be solved on a doubly periodic domain C. The boundary conditions on the positively oriented interface between the inclusion phase and the matrix phase are given in section 2. An average electric field E 0 of unit strength is applied. The permittivity of the matrix phase is set to 2 = 1 so that the effective permittivity is equal to the effective relative permittivity. From the repeated solution to this problem for various 1 = σ , we obtain eff (σ ). Two types of domains are investigated: the 'square array of square cylinders' and the 'staggered array of square cylinders' (a square array of cylinders with diamond-shaped cross-sections, as studied for example in [35,36]); see figure 2.
Boundary value problems on domains involving sharp corners may require extreme resolution close to corner vertices, even when the demands for overall accuracy are moderate. One has to be selective with the choice of numerical method. We chose an integral equationbased scheme. Such schemes have the advantage that they can also retain stability in very difficult situations.
Our particular choice of integral equation is standard-a single-layer equation [37]. For its solution we use a novel numerical method called recursive compressed inverse preconditioning.  Conceptually this is a local multilevel technique that makes a change of basis and expresses the non-smooth solution to the single-layer equation in terms of a piecewise smooth transformed layer density that can be cheaply resolved by polynomials. Discretization leads to a block diagonal transformation matrix R (an inverse preconditioner) where the columns of a particular block can be interpreted as special basis functions for the original density in the vicinity of a corner vertex multiplied by suitable quadrature weights. The blocks of R are constructed in a fast recursion, i = 1, . . . , n, where step i inverts and compresses contributions to R involving the outermost quadrature panels on level i of a locally n-ply refined mesh. We emphasize that the method is strictly numerical and fully automatic. There is no separation of variables or eigenvalue analysis involved. The recursive compressed inverse preconditioning method was originally described in [38] and further developed in [39][40][41]. The appendix below highlights some of the method's features, relevant to the domain C. A fuller description will be included in a forthcoming paper [42].

Achievable accuracy
We first compute eff (σ ) for the square array of squares at p 1 = 0.25 in the limit of σ approaching the negative real axis from the upper half-plane H, as in the example of figure 1. The exact result (11) is used as a benchmark. Figure 3 shows that the relative error is close to machine epsilon (the upper bound due to rounding in floating point arithmetic) except for a neighbourhood of three points where it is higher: the ends of the branch-cuts at σ = −3 and σ = −1/3 and at the singularity of the integral equation at σ = −1. This demonstrates that the problem of computing eff (σ ) for arrays of square cylinders is well conditioned in general and that our scheme is stable.  The staggered array of square cylinders at p 1 = 0.499 999 999 95 is a more challenging geometry than the square array of square cylinders at p 1 = 0.25: • there are more length scales involved, • eff (σ ) varies more rapidly and has more poles and zeros, • there is no exact result to compare with. The first problem is the least difficult. The multilevel property of our numerical method should enable the resolution of almost arbitrarily small separation distances between corner vertices. The second problem is more severe. For eff (σ ) close to zero, one can expect cancellation in (A.2) and the relative accuracy should suffer. Furthermore, it is harder to resolve wildly varying functions in floating point arithmetic than slowly varying ones. The third problem is solved by using the extent to which (3) is satisfied as an indicator of the relative error. For this, since σ and 1/σ lie on different sides of the real axis and our numerical method takes limits from H, we use (3) in the equivalent form where the ' * ' symbol denotes complex conjugation. Figure 4 suggests that despite the difficulties, typically, only a few digits are lost compared to the square array at p 1 = 0.25. The numerics seem to give a relative precision of at least 10 −8 even for the most extreme values of eff (σ ).
For the interpretation of various limits it might also be of interest to study eff (σ ) for σ some finite distance into H. Figure 5 shows again the staggered array of square cylinders at p 1 = 0.499 999 999 95, but unlike in figure 4 we have here interrupted the limit process at σ a relative distance of 10 −5 away from the real axis.

Animations of spectral evolution
We now consider the evolution of the variation of eff (σ ) with permittivity ratio σ , as the area fraction of square inclusions p 1 in the square array ranges from zero to unity. This evolution   is shown in animation 1, from which a typical frame is given in figure 6. The most important feature of animation 1 is quite clear: for σ real and for all values of p 1 , except at the poles, non-zero values of { eff (σ )} are confined to the interval −3 σ −1/3, in agreement with the suggestion by Hetherington and Thorpe [3]. Of course, the value of this imaginary part is always positive if we restrict ourselves to composites without gain (for which the imaginary part would be always negative).
Below the area fraction of 0.25, the real part of eff (σ ) is positive, and it develops its first pole at this area fraction. It is interesting to compare frames from animation 1 and the left image of figure 3 with figure 1; the results of the mode-matching method clearly correspond to those of the new method, but are capable of a resolution limited by the number of terms employed in field expansions.
For p 1 > 0.25, the real part is negative between σ = −3 and σ = −1/3, while the pole migrates to more negative values of σ . For area fractions near p 1 = 0.75 (see figure 6), 'features' that we call quasipoles develop from near σ = −1, and one moves towards σ = −3, while the other moves towards σ = −1/3. When they reach these values and then are not muted by the absorbing nature of the corners, they transform into actual poles, which move towards σ = −∞ and σ = 0, respectively. At higher values of area fraction, more quasipoles evolve from σ = −1 and give rise to additional actual poles when they emerge from the branch-cut region. { eff (σ )} becomes small as p 1 → 1, while { eff (σ )} tends towards σ , apart from the increasingly numerous but increasingly narrow pole regions.
For the staggered array of square prisms, animation 2 illustrates the behaviour of eff (σ ) as a function of σ for area fractions ranging from zero to 0.5 (with the behaviour for p 1 in the range 0.5-1.0 following from that in the lower range using Keller's theorem (3)). As has been commented in section 2, the interesting question is how the branch-cut location (from σ = −∞ to σ = 0) for p 1 = 1/2 of equation (6) can be reconciled with that (from σ = −3 to σ = −1/3) predicted by Hetherington and Thorpe [3] for p 1 arbitrarily near 1/2. A mechanism for this reconciliation was provided by one of the present authors [4]: discrete sets of poles in −∞ σ −3 and 1/3 σ 0 were predicted to become denser and denser as p 1 approached 1/2, thus extending the branch-cut in the limit to that required by (6). The accuracy of this prediction is evident in Animation 2: poles develop from the quasipoles generated at σ = −1 and move left and right into the embryonic branch-cut regions −∞ σ −3 and 1/3 σ 0. The left frame in figure 4 shows a stage in this evolution where p 1 is very close to 1/2. Animation 2 makes the 'nursery role' of the region around σ = −1 in the development of the spectrum much more evident (due to larger amplitudes of the quasipoles) than does animation 1.
The sensitivity of the spectral details for the staggered array near the checkerboard configuration is very evident in animation 2. Even when the inclusions are extremely close to touching, the poles outside the branch-cut are still fairly well spaced apart, and so the absorption spectrum acts like a Vernier scale for determining the inclusion separation. It would be interesting to obtain an asymptotic formula for how the distance between poles depends on the inclusion separation distance.
As we have commented in section 2, the asymptotics of fields near corners are the same in electromagnetism as in electrostatics [14]. Thus, attempts such as that in [33] to model the transition from electromagnetically reflecting structures to electromagnetically transmitting structures as p 1 moves through 1/2 would require an adaptive and recursive method such as that described in section 3 to be able to achieve sufficient accuracy.

Discussion and conclusions
In this paper, we have brought together rigorous mathematical results with numerical investigations of unprecedented accuracy. The latter have revealed the generality of the former, and have substantiated a conjecture of Hetherington and Thorpe [3] in a striking and conclusive way.
We conclude by commenting further on how the arguments and results we have presented can be implemented in a practical demonstration of morphological super-resolution, uniting the ideas of Kac [1] and Pendry [2]. Such a demonstration would require the fabrication of a set of parallel cylinders with a square cross section (or polygonal cross section). The cylinders do not have to be arranged in a geometrically perfect array, and they do not have to be densely packed. They have to be made of a material that is essentially non-absorbing and with a negative permittivity or permeability.
These requirements suggest the set of cylinders be made of metal, of size 10 µm or larger and be probed with wavelengths far greater than the cylinder size in the far-infrared or longer. Such cylinders are large by today's lithographic standards, and so it should be possible to accurately form their corners to achieve sub-wavelength accuracy. Going into the far-infrared region diminishes metallic loss from its value in the visible and near-infrared [45]. It is crucial that the metallic loss be very low, since the experimental signature we suggest be probed is enhanced absorption by a set of such cylinders over a wavelength interval in which the metal's permittivity ranges from say −1/3 to −3 (scaled relative to the permittivity of a host dielectric in which the cylinders are embedded). Note that the enhanced absorption of incident radiation detected will increase as more lines of cylinders are added to the set.
The experimental result that would indicate morphological super-resolution is an enhanced absorption for wavelengths far in excess of the cylinder size, switching on and off at geometrically determined limits described above, independent of the arrangement and area fraction of the cylinders. We stress, however, that such a demonstration would indicate the physical relevance of the ideas we have described for a particular system. The mathematical results we have described are of course rigorous, and the numerical examples of them we have given are highly accurate, so our demonstration of super-resolution for wavelengths arbitrarily larger than the size of the particles probed does not rely for its validity on experimental support. They may be applicable to governing equations other than the Helmholtz equation, for which the ideas of metamaterials and their applications are currently being explored [46].
Here n z is the outward unit normal of at z, 0 denotes the restriction of to U and λ is as in (7). Note that, as σ → −1, we have λ → ±∞ and (A.1) is no longer an equation of the second kind, but an equation of the first kind whose (unique) solvability is by no means guaranteed. Therefore one can say that σ = −1 is a singularity of (A.1).
Once (A.1) is solved for ρ(z) and under the assumption that the inclusions do not overlap the unit cell boundary, the effective relative permittivity in the direction of the applied electric field can be computed from Depending on how the unit cell is chosen, the squares in the staggered array may overlap the unit cell boundary. With the choice in figure 2, they certainly do. But since ρ(z) is a periodic function and identical on all squares, one can circumvent this problem by modifying (A.2) so that it integrates ρ(z) twice on the square at the centre of the unit cell and ignores ρ(z) on the other squares.

A.2. Discretization
We discretize (A.1) and (A.2) using a Nyström method based on composite 16-point polynomial interpolatory quadrature and a parameterization z(t) of . The parameter t is real. See [44] for a review of Nyström methods including error analysis.
An initial coarse mesh that resolves the kernel of the integral operator in (A.1) away from the corner vertices is constructed on ; see figure 2. The coarse mesh is refined by subdividing those panels that neighbour corner vertices. The subdivision is done n times in a direction towards the vertices. On quadrature panels that neighbour corner vertices we choose quadrature nodes according to the zeros of certain Jacobi polynomials. On the remaining panels we choose quadrature nodes according to the zeros of Legendre polynomials. Upon discretization on the refined mesh, (A.1) assumes the form where I fine and K fine are square matrices and ρ fine and g fine are column vectors. The vector g fine corresponds to the discretization of the piecewise smooth right-hand side. This right-preconditioned equation corresponds to the discretization of a Fredholm equation of the second kind with compact operators. The solutionρ fine is the discretization of a piecewise smooth function.

A.3. Compression
The matrix K • fine , the densityρ fine and the right-hand side g fine in (A.6) can be evaluated on the coarse mesh without loss of precision. Only (I fine + K fine ) −1 needs the refined mesh for its accurate evaluation. This enables a compression of (A.6). We introduce the compressed weighted inverse R = P T W (I fine + K fine ) −1 P. (A.7) Here P is an unweighted prolongation operator that performs panelwise 15th-degree polynomial interpolation in the parameter t (which as we recall parameterizes through z(t)) from points on the coarse mesh to points on the fine mesh when acting on column vectors from the left. P W is a weighted prolongation operator. See section 5 of [40]. Substitution of (A.7) into (A.6) and the use of some relations between prolongation operators make (A.6) assume the form This equation, defined solely on the coarse mesh, is used in our computations.

A.4. Fast recursion for R
The construction of R from its definition (A.7) may be a costly and unstable operation when the refined mesh has many panels. The number of subdivisions n needed to reach a given accuracy may grow without bounds due to the singularities in ρ(z) that arise as σ approaches certain values. Fortunately, the construction of each block of R, associated with a corner of the square array or with a corner-meet of the staggered array, can be greatly sped up and also stabilized via a recursion. This recursion uses matrices K on local meshes centred around corners or cornermeets. It would be going too far to describe all the fine details of this procedure, but sections 3.2 and 3.3 of [41] give a fairly good idea of how the recursion is set up in the present context. A key step is the (partial) conversion of the recursion into a nonlinear matrix equation. This equation is solved using a variant of Newton's method relying on numerical homotopy to approach purely negative σ from the upper half-plane H. The ratio σ , which enters into K, is initially multiplied by a constant q = 1 − 0.01i. The imaginary part of q is reduced by a factor of ten after each of the first 14 Newton iterations. Then q is set to unity and the iterations are continued until either a sharp convergence criterion is met or a total of 30 iterations is reached. A full description is given in [42].