Macroscopic optical response and photonic bands

We develop a formalism for the calculation of the macroscopic dielectric response of composite systems made of particles of one material embedded periodically within a matrix of another material, each of which is characterized by a well defined dielectric function. The nature of these dielectric functions is arbitrary, and could correspond to dielectric or conducting, transparent or opaque, absorptive and dispersive materials. The geometry of the particles and the Bravais lattice of the composite are also arbitrary. Our formalism goes beyond the longwavelenght approximation as it fully incorporates retardation effects. We test our formalism through the study the propagation of electromagnetic waves in 2D photonic crystals made of periodic arrays of cylindrical holes in a dispersionless dielectric host. Our macroscopic theory yields a spatially dispersive macroscopic response which allows the calculation of the full photonic band structure of the system, as well as the characterization of its normal modes, upon substitution into the macroscopic field equations. We account approximately for the spatial dispersion through a local magnetic permeability and we analyze the resulting dispersion relation, obtaining a region of left-handedness.


Introduction
The propagation of light within homogeneous materials is completely characterized by their electromagnetic response. In the optical regime, magnetic effects are typically negligible so that knowledge of the dielectric response is sufficient for the study of wave propagation [1]. However, when the system has spatial inhomogeneities, the behaviour of light is not trivial and has captured the attention of many researchers [2,3,4,5,6,7,8,9,10,11]. The possibility of controlling the propagation of photons by producing artificially structured materials has been widely demonstrated in recent years. Novel effects have been obtained, such as negative refraction, inverse Doppler effect, optical invisibility cloaking, optical magnetism etc. [12,13,14,15].
In this paper we obtain the macroscopic optical response using a computationally efficient reformulation of a procedure [5,6] original developed to account for the local field effect of systems with spatial fluctuations. Our main result is that for any system, the macroscopic response to an external excitation is simply the average projection of the corresponding microscopic response. We apply this general result to Maxwell equations in order to obtain an explicit expression for the macroscopic dielectric response of binary composites. We illustrate its use with numerical results for a 2D periodic photonic crystal. As the wavelength of the electromagnetic field becomes comparable to the characteristic microscopic spatial lengthscale of the system, the macroscopic response becomes highly non-local, with a strong dependence on the wavevector besides its usual dependence on frequency. When the full spatial dispersion is taken into account, the photonic band structure of the system can be obtained and analyzed from its macroscopic response.
The paper is organized as follows: In Sec. 2 we formulate a general theory for the macroscopic dielectric response of composite systems. In Sec. 3 we adapt our formulation to periodic systems made of two alternating materials, each characterized by well defined dielectric functions. In Sec. 4 we develop an efficient computational scheme which allows the numerical calculation of the macroscopic response without recourse to explicit operations between large matrices. In Sec. 5 we obtain the nonlocal macroscopic dielectric tensor and use it to calculate macroscopically the full band structure of a 2D photonic crystal made up of non-dispersive dielectric components; we compare it to exact results as well as to results within a local approximation that partially accounts for spatial dispersion through an effective magnetic permeability. Finally, in Sec. 6 we present our conclusions.

General theory
Consider a non-magnetic inhomogeneous system characterized by its dielectric responseǫ, defined through the constitutive equation D =ǫE, where D and E are the microscopic displacement and electric fields respectively. We designate these fields as microscopic as they have spatial fluctuations due to the inhomogeneities of the material. They are not microscopic in the atomic scale, but rather, on the scale of the spatial inhomogeneities of the system. We use a caret (ˆ) over a symbol to denote its operator nature and we leave implicit the dependence of the dielectric response on position as well as on frequency. Our purpose is to obtain the macroscopic dielectric response of the system, relating the macroscopic displacement and electric fields, from which the fluctuations have been removed.
We start from Maxwell equations for monochromatic microscopic fields. We follow the usual procedure [36] to decouple the magnetic field from the electric field to obtain the second order wave equation is the wave operator, q = ω/c is the wavenumber of light in vacuum, ω is the frequency, c is the speed of light in vacuum and J is the external electric current density. Here we introduced the transverse projectorP T such that the transverse projection of a field F is F T =P T F. For completeness, we also introduce the longitudinal projectorP L such thatP L +P T =1, with1 the identity operator. Notice that Eq. (1) contains both a longitudinal and a transverse part, so it describes not only transverse waves, but also allows for the possible excitation of longitudinal waves such as plasmons [37]. We solve Eq. (1) formally for E to obtain As we are interested on the macroscopic response of the system, we introduce the averageP a and fluctuationP f projectors, such that an arbitrary field F can be written as F = F a + F f with F a =P a F its average and F f =P f F its fluctuating part, and we identify the average field with the macroscopic field F a ≡ F M . Later, we will provide appropriate explicit definitions forP a andP f ; here we remark that they must be idempotentP 2 a =P a ,P 2 f =P f , i.e., the average of the average is the average, and the fluctuations of the fluctuations are the fluctuations [5]. Furthermore, P aPf =P fPa = 0 andP a +P f =1.
As the macroscopic response would be useless unless the external excitations are devoid of microscopic spatial fluctuations, we assume that the external current has no fluctuations, so J f = 0 and J = J M =P a J M , where we used the idempotency ofP a . Thus, acting on both sides of Eq. (3) withP a we obtain Here we defineÔ αβ ≡P αÔPβ , with α, β = a, f for any operatorÔ. As Eq. (4) relates the macroscopic external electric current to the macroscopic electric field we may identify the macroscopic inverse wave operatorŴ −1 M , We summarize this results stating that the macroscopic inverse wave operator is simply the average of the microscopic inverse wave operator. This is a particular case of a more general result: the macroscopic response to an external excitation is simply the average of the corresponding microscopic response. From the macroscopic Maxwell equations we may further relate the macroscopic wave operator in Eqs. (5) and (6) to the macroscopic dielectric response of the system ǫ M throughŴ in analogy to Eq. (2). Thus, we have to invert the wave operator, average it, and invert it again to finally identify the macroscopic dielectric response. Notice that we have made no approximation whatsoever. We remark that asǫ relates two fields, E and D, that have spatial fluctuations, we may not simply averageǫ to obtain its macroscopic counterpart, i.e.ǫ M =ǫ aa . The difference constitutes the local field effect [5].
Our result may easily be generalizable to other situations and other response functions. The procedure consists on first identifying the response (in our caseŴ −1 ) to the external perturbation (i.e. J), and then averaging it to yield the corresponding macroscopic response (i.e.Ŵ −1 M =Ŵ −1 aa ), which may further be related to the desired macroscopic response operator (i.e.ǫ M ).

Periodic binary systems
In this section we use Eq. (6) to obtain the optical properties of an artificial binary crystal made of two materials A and B with dielectric functions ǫ A and ǫ B . We assume that both media are local and isotropic so that ǫ A and ǫ B are simply complex functions of the frequency. For convenience, we will further assume that ǫ A is real, though this assumption is easily relaxed [38].
We introduce the characteristic function B(r) of the inclusions, such that B(r) ≡ 1 whenever r is on the region B occupied by the inclusions, and B(r) ≡ 0 otherwise. Thus, we may write the microscopic dielectric response as where we defined the spectral variable u ≡ 1/(1 − ǫ B /ǫ A ) [17]. The microscopic wave operator of Eq. (2) may be written aŝ where the characteristic operatorB corresponds to multiplication by B(r) in real space, and we defined which, as shown below, plays the role of a metric. Using Eq. (9) we write Eq. (6) aŝ where we have taken advantage of the fact thatĝ is unrelated to the texture of the crystal, so that it does not couple average to fluctuating fields. Using Bloch's theorem [39], we consider an electric field of the form [19] where k is a given wavevector, {G} is the reciprocal lattice of our crystal and the coefficients E G represent the field in reciprocal space. In this representation, all operators may be written as matrices with vector index pairs G, G ′ , besides other possible indices, such as Cartesian ones. Thus, we represent the longitudinal projector as the matrix with δ GG ′ the Kronecker's delta, so the transverse projector becomes P T GG ′ = 1δ GG ′ − P L GG ′ with 1 the Cartesian identity matrix. The Laplacian operator in reciprocal space is represented by and we can define the average projector as a cutoff in the reciprocal space [6] ( so that average fields simply keep the contributions with vector k while all other wavevectors are filtered out. As shown by Eq. (11), we only require to obtain the macroscopic inverse wave operator, where the subindices 0 denote the projection onto the subspace with G = 0.

Recursive method
The calculation of Eq. (16) is analogous to that of a projected Green's function [40,41] G aa (ε) = a|Ĝ(ε)|a , onto a given state |a , wherê is the Green's operator corresponding to some Hermitian HamiltonianĤ and ε is a complex energy. In Ref. [26] a similar result was obtained, where the HamiltonianĤ was identified with the longitudinal projection of the characteristic functionB LL , the energy ε with the spectral variable u, and |a with a slowly varying longitudinal wave, and it was shown that the corresponding projected Green's function was proportional to the inverse of the longitudinal macroscopic dielectric function. Haydock's method may be applied to obtain the projected Green's function (17) in a very efficient way [42,43], and it has been adapted previously [26] to the calculation of the optical response of nanostructured systems in the long-wavelength local limit. In this section we generalize the approach of Ref. [26] to arbitrary frequencies and wavevectors. Eqs. (17) and (18) are similar to Eq. (16), although the operator product H =Bĝ that plays the role of Hamiltonian does not correspond to a symmetric matrix. Nevertheless, it would correspond to a Hermitean operator if we useĝ as a metric, that is, if we use ψ|ĝ|φ instead of ψ|φ as the scalar product of two arbitrary states |ψ and |φ . Then, it is easy to verify the Hermiticity condition [ ψ|ĝ(Ĥ|φ )] * = φ|ĝ(Ĥ|ψ ). Notice, however, thatĝ is not positive definite.
According to Eq. (11), we need the average (u −Bĝ) −1 aa → a|(u −Bĝ) −1 |a , where we projected onto an average state |a = b 0 |0 consisting of a plane wave with a given wavevector k, frequency ω and polarization e. Here, |0 is normalized according to the metricĝ, i.e., 0|ĝ|0 = g 0 = ±1, and the coefficient b 0 is a real positive number chosen to guarantee that |a is normalized in the usual sense a|a = 1. Now, we define | − 1 ≡ 0 and we obtain new states by means of the recursion relation | n + 1 ≡Bĝ|n = b n+1 |n + 1 + a n |n + g n−1 g n b n |n − 1 , where all the states |n are orthonormalized according to the metricĝ, that is with g n = ±1 and δ nm the Kronecker's delta function. The requirement of orthonormality yields the generalized Haydock coefficients a n , b n+1 , and g n+1 given the previous coefficients b n , g n and g n−1 . Thus, a n are obtained from n|ĝ| n + 1 = a n g n , and b n+1 from where we choose the sign g n+1 = ±1 so that b 2 n+1 is positive and we may choose b n+1 as a real positive number. In the basis {|n },Bĝ is represented by a tridiagonal matrix with a n along the main diagonal, b n along the subdiagonal and g n−1 g n b n along the supradiagonal, so that Notice that we may represent the states |n through the corresponding Cartesian field components for each reciprocal vector G or for each position r within a unit cell. Thus, the action ofĝ is a trivial multiplication in reciprocal space, while that ofB is a trivial multiplication in real space, and we may repeatedly computeBĝ|n and Haydock's coefficients without performing any actual matrix multiplication, by alternatingly fast-Fourier transforming our representation between real and reciprocal space. According to Eq. (11), we do not require the full inverse of the matrix (23), but only the element in the first row and first column. Following Ref. [26], we obtain that element as a continued fraction, which substituted into Eq. (11) yieldŝ The right hand side of Eq. (24) denotes the macroscopic response projected onto a state with the given wavevector k, frequency ω and polarization e. Due to the translational invariance of the homogenized system, the inverse wave operator for a given k corresponds simply to a tensor with Cartesian components (W −1 M ) ij . Thus, so far we have calculated its inner products with e, as shown by the second term of Eq. (24). We may repeat the calculation above for different orientations of the (possibly complex) unit polarization vector e and solve the resulting equations for all the individual Cartesian components (W −1 M ) ij . Finally, we perform a matrix inversion of this macroscopic tensor and we use Eq. 7 and 6 to compute the macroscopic dielectric tensor Eqs. (24) and (25) constitute the main result of our formalism. Notice that in general, this dielectric tensor would depend on both ω and k, that is, it is a temporal and spatially dispersive response function.

Results
To test our formalism, we first calculate the local limit of the longitudinal macroscopic dielectric function ǫ M L (ω) ≡ ǫ M xx (ω, k → 0x) of a 2D system made up of a square array of empty cylindrical holes (ǫ B = 1) placed with lattice parameter a within a dielectric medium with permittivity ǫ A = 12. The holes are of circular cross section with radius ρ = 0.45a. We assume that the axes of the cylinders are parallel to the z axis. This system has been frequently been used as a testground for calculations of photonic crystal properties [44]. In this calculation we took the limit k → 0 along the x direction, and we introduced the Cartesian unit vectors x, y and z. We remark that it is important to assign a direction to k even in this limit, as the transverse T and longitudinal L projection of W M differ, although, for this system, the resulting ǫ M (ω) = ǫ M (ω, k → 0) is isotropic within the xy plane.
In Fig. 1 we show ǫ M L (ω) calculated for two different lattice parameters a = 40nm and a = 120nm. Our calculations were performed using the Perl Data Language [45] ‡ . The figure shows that our formalism yields the same results as a full straightforward matricial calculation [30] that solves Maxwell equations in a plane wave basis. Nevertheless, our calculation is much more economical in memory usage, as we don't have to store the full matrices that represent the dynamics of the system, and is at least four orders of magnitude faster as we do not perform any matrix manipulation with our scheme. To illustrate the effects of retardation, in the figure we also show the result of using the non-retarded calculation of Ref. [26] (NR), together with those of the well-know 2D Maxwell-Garnet (MG) effective medium formula [19]. As the dielectric functions of the components are independent of the frequency, since the constituents were taken to be non-dispersive, in this case the frequency dependence of the result arises solely from to the finite ratio of the free space wavelength to the lengthscale of the system, namely, the lattice parameter a. Thus, in the ω → 0 limit our results approach those of the NR calculation. The difference between the NR T (ω, k ∆ ) of the same system as in Fig. 1 for a given wavevector k ∆ along the ∆ line of the first Brillouin zone (solid lines). The line (k ∆ /q) 2 (long dashes) intersects (circles) ǫ M T (ω, k ∆ ) at the frequencies ωα(k ∆ ) of the TE normal modes, with α =I, II, IV, V. The intersections (diamonds) of ǫ M T (ω, k ′ ∆ ) (short dashes) with (k ′ ∆ /q) 2 (long dashes) corresponding to a shifted wavevector k ′ ∆ = k ∆ + (2π/a)y yields the modes α =III, IV (see text). and MG calculations is due to the high filling fraction f ≈ 0.64 of the cylinders, which interact non-dipolarly with their nearest neighbors which they almost touch, while MG contains only dipolar interactions. Naturally, retardation effects become stronger as a increases. Furthermore, ǫ M (ω, k) depends only on the products qa and ka so that the two curves shown in Fig. 1 could be collapsed into one curve if appropriately plotted. We use this in Fig. 2, which shows the transverse response ǫ M T (ω, k ∆ ) = ǫ M zz (ω, k ∆ ) corresponding to out-of-plane or TE polarization as a function of the normalized frequency qa/2π. Here, we have chosen a finite in-plane wavevector k ∆ = (π/2a, 0, 0), corresponding to the midpoint of the ∆ line between Γ and X in the first Brillouin zone (BZ) [39]. This figure covers a larger frequency range than Fig. 1, and therefore it displays a series of poles related to resonant multiple coherent reflections at the interfaces between the A and B regions.
The dispersion relation of transverse waves with TE polarization propagating along the xy plane can be simply written as [46] Therefore, in Fig. 2 we also plotted the curve (k ∆ /q) 2 . Those points where it intersects the curves for ǫ M T correspond to normal TE modes of the system with frequencies ω α (k = k ∆ ), α =I, II, . . . Here, we use roman numerals to number the modes in increasing frequency order.
In Fig. 3 we show the photonic bands (circles) of the TE modes of our system obtained by solving the dispersion relation (26)  plot the TE modes obtained by solving the eigenvalue equation [47] 1 ǫ(r) through an implicitly restarted Arnoldi iteration [48] , applying 1/ ǫ(r) in real space and ∇ → i(k + G) in reciprocal space [38]. We further verified our results by recalculating them with the freely available package MPB [50] and comparing it to previous results [44] for the same system. The agreement between these calculations shows that it is feasible to use the macroscopic dielectric response of the system to obtain its photonic band structure. However, to succeed, we have to account fully for the spatial dispersion of ǫ M T (ω, k). We remark that although we obtained the correct photonic bands even for large wavevectors, going all the way around the first BZ, there are some regions where we failed to obtain all of the normal modes. For example, our procedure did not produce the third band (labeled III in the figure) along the ∆ line, and neither the fourth band (IV) along the Σ line, between M and Γ. Furthermore, it did not yield the second band (II) which is almost degenerate with the third band (III) along Σ. It is easily shown that the microscopic field corresponding to the missing band along ∆ is antisymmetric with respect to the reflection G y ↔ −G y about the ∆ line. Similarly, the missing bands along the region Σ are antisymmetric with respect to the reflection G x ↔ G y about the Σ line. Thus, the reciprocal vector 0 does not contribute to the electric field for those bands [51], i.e., the missing modes have no macroscopic field components and thus, apparently they may not be obtained from the roots of the macroscopic dispersion relation (26).
Nevertheless, within our formalism, the wavevector k is actually not restricted to lie within the first BZ. Thus, we may calculate the macroscopic dielectric function for wavevectors beyond the first BZ, and obtain the modes in the extended scheme. In Fig.  2 we show part of the macroscopic dielectric function ǫ M T (ω, k ′ ∆ ) with a wavevector that has been shifted out of the first BZ from k ∆ along the y direction by the reciprocal vector (2π/a)y, i.e., k ′ ∆ = k ∆ + (2π/a)y. We remark that ǫ M (ω, k) is not a periodic function of k, as it corresponds to a specific plane wave, not to a microscopic Bloch's wave. The intersection of the curve (k ′ ∆ /q) 2 with ǫ M T (ω, k ′ ∆ ) yields the corresponding macroscopic modes. For example, the diamonds in Fig. 2 illustrate two of the normal modes that could be obtained using this shifting procedure. One of these modes is identical to the mode labelled IV that was obtained previously using the unshifted response ǫ M T (ω, k ∆ ). Nevertheless, another mode, labelled III appears just below qa/2π = 0.4 and it corresponds to one of the previously missing modes. We may shift it back into the first BZ in order to display it in the reduced zone scheme. Proceeding in this fashion, we have obtained all of the previously missing bands, shown by diamonds in Fig. 3. Thus, we have shown that we can calculate the full photonic band structure for the TE modes from the macroscopic non-local dielectric function of the composite system.
A similar procedure to that discussed above can be employed also for the TM modes of the system, with the electric field within the xy plane. In Fig. 4 we show the in-plane macroscopic dielectric functions ǫ M xx (ω, k ∆ ) and ǫ M yy (ω, k ∆ ) of the system for the same wavevector k ∆ as in Fig. 2. This wavevector and the system are symmetric under y ↔ −y reflections. Therefore, ǫ M xy = 0 and we may identify the transverse and longitudinal responses as ǫ M T (ω, k ∆ ) = ǫ M yy (ω, k ∆ ) and ǫ M L (ω, k ∆ ) = ǫ M xx (ω, k ∆ ) respectively. The normal modes are then obtained from the transverse dispersion relation (26) and the longitudinal dispersion relation [46] ǫ M L (ω, k) = 0.
These are identified with circles and diamonds in the figure.
Along the Z line, from X to M , there is no crystal symmetry operation that leaves invariant the wavevector, and thus transverse and longitudinal fields mix among themselves, i.e., there are no longitudinal nor transverse modes. Nevertheless, we may obtain the frequencies of the mixed modes through the singularities of the wave operator matrix as illustrated in Fig. 6 for a wavevector k Z = (π/a, π/2a) along the Z line. From calculations such as those illustrated by Figs. 4, 5 and 6, or more generally, from the zeroes of det W M (ω, k) for arbitrary wavevectors k, we can obtain the full TM photonic band structure, as shown by Fig. 7. In order to test our theory, we also plot the TM modes obtained by solving the eigenvalue equation [47] ∇ · 1 ǫ(r) using similar methods [48,38] as for the TE case and comparing them successfully to previous results [44] for the same system. The agreement between these calculations shows that it is also feasible to use the macroscopic dielectric response of the system to obtain its photonic band structure in the TM case.
An approximate band structure could be obtained for transverse waves in the region of small wavevectors k → 0 by expanding the LHS of Eq. (26) up to second order in k and solving for k 2 . The result is a local dispersion relation where the non-locality of ǫ M T (ω, k) is partially accounted for through an effective magnetic permeability [23] µ(ω) = 1 In Fig. 8 we show the resulting approximate bands for the TE case. We can see that the acoustic band within the local approximation is indistinguishable from the exact result for small wavevectors. This local acoustic band extends up to the frequency qa/2π ≈ 0.2 for which µ has a pole. The local approximation also reproduces the upper band of the exact calculation (labelled V in Fig. 3), but only in the immediate vicinity of Γ, close to the zero qa/2π ≈ 0.43 of ǫ M (ω). For larger wavevectors it acquires a larger positive dispersion. As discussed above, the band labelled III in Fig. 3 doesn't couple to macroscopic fields for k along the ∆ line, while the band labelled IV doesn't couple along the Σ line. Thus, neither band couples to macroscopic fields at Γ and as a consequence, these bands are not reproduced at all by the local approximation (31).
In the region qa/2π ≈ 0.36 − 0.40 the local approximation predicts a curious band that displays backbending. Its lower part begins at the pole of ǫ M (ω), which corresponds to a double zero of µ(ω) [38]. This part of the local dispersion relation is spurious, as the second order Taylor expansion of the dielectric function that yields Eqs. (32) and (31) starting from Eq. (26) would be meaningless at a pole. Consequently, this part of the local band doesn't correspond to any band within the exact calculation. On the other hand, the upper part of the backbending band starts at a simple zero of µ which arises from a different kind of singular behavior, consisting of a pole in ǫ M T (ω, k) at qa/2π ≈ 0.40 which is suppressed at k = 0 as its weight is approximately proportional to k 2 . This part of the local dispersion relation agrees with band II of Fig. 3 at Γ and displays a negative dispersion as does band II. Furthermore, for this band, both ǫ M and µ are negative. Thus, our photonic crystal, made up of holes within a dispersionless dielectric, behaves for some frequencies as a left-handed metamaterial [52,53]. Nevertheless, caution should be exercised when using Eqs. (32) and (31) close to singularities of ǫ M T (ω, k).

Conclusions
We have shown that the macroscopic inverse wave operator of composite systems with arbitrary spatial fluctuations is simply given by the average projection of the corresponding microscopic operator. These operators may be related to the macroscopic and macroscopic dielectric functions, thus yielding a general procedure to incorporate the effects of spatial fluctuations in the calculation of all the components of the macroscopic dielectric tensor of the system. We have extended Haydock's recursive scheme in order to calculate very efficiently the response of periodic two-component systems in terms of a continued fraction whose coefficients may be obtained without recourse to operations with the large matrices that characterize the microscopic response and the field equations. We illustrated and tested our results through the calculation of the wavevector and frequency dependent macroscopic dielectric tensor of a 2D dielectric system with non-dispersive non-dissipative components. We identified the transverse and longitudinal components of the response for wavevectors along symmetry lines and we obtained all the components of the tensorial response for the case of reduced symmetry. The macroscopic response has a series of poles, related to resonant multiple coherent reflections and by accounting for its spatial dispersion we obtained the full photonic band structure of the system from the dispersion relations corresponding to the homogenized material. Besides yielding the correct bands from a macroscopic approach, our scheme allowed us to classify the polarization of each mode as either transverse, longitudinal or mixed. Even though the macroscopic approach might fail to yield some modes, namely, those which have no coupling to the macroscopic field due to symmetry, they may be recovered by extending the notion of macroscopic state, allowing its wavevector to lie beyond the first Brillouin zone. The non-locality of the dielectric response may be partially accounted for through a magnetic permeability which can then be employed in the calculation of the optical properties of the system. We compared the band structure obtained through this local approximation to the exact results and we obtained partial agreement close to Γ for those bands that do couple to long-wavelentgh fields. In particular, we showed that this approach can yield a negative dispersion in frequency regions where both the permitivity and permeability are negative. Nevertheless, we discussed some shortcommings of the local approach. Although here we presented results for dispersionless transparent dielectrics, our calculation does not require that the materials that make up the system be non-dispersive nor non-dissipative, so that calculations for real dielectrics and metals may be performed with the same low computational costs.