1. INTRODUCTION

This is the second part of the tutorial on Maxwell Garnett approximation. The first, introductory part was published in [1], and we will refer to the equations of [1] in the form (FP#), where # is the equation number.

In the second part, several more advanced topics will be discussed. The overall goal is to justify and generalize the Maxwell Garnett approximation and to place it more solidly in the context of modern classical electrodynamics of continuous media. In particular, we will take into consideration the effects of finite frequencies and electromagnetic interaction of inclusions. We will also develop a rigorous homogenization theory for periodic composites. We will see that the Maxwell Garnett mixing formula is obtained from this result by neglecting a term that is responsible for the interaction of inclusions.

We will start with laying out in Section 2 a general theoretical framework of the macroscopic Maxwell’s equations in continuous media. This old topic has attracted renewed attention recently [2,3]. Traditionally, the derivation of the macroscopic Maxwell’s equations from the microscopic ones was done primarily for pedagogical purposes and the averaging procedure typically involved in such derivations was not used constructively. However, the development of modern computational approaches to the microscopic theory and, especially, the emergence of the density-functional theory (DFT)-based ab initio methods for computing the microscopic electromagnetic quantities have motivated re-examination of some basic assumptions.

One important idea that requires a fresh look is the classical interpretation of the polarization and magnetization fields as the volume densities of the electric and magnetic dipole moments. Ab initio calculations of the microscopic-induced charge density in some solids such as crystalline silicon have shown that this interpretation can be incorrect [4]. In Section 3, we discuss the range of applicability of the classical interpretation and identify two main reasons for its breakdown: charge leakage through the cell boundaries (in crystalline materials) and accumulation of a nonnegligible electromagnetic phase between the neighboring cells. The two causes for the classical interpretation breakdown become important in different physical situations. Charge leakage can be relevant for natural materials in which the elementary cells are microscopic but not for macroscopic composites. On the other hand, the phase shift plays almost no role in natural materials but can come to the fore in the macroscopic homogenization theories.

In Section 4, we consider propagation of waves on three-dimensional cubic lattices of polarizable particles. It is shown that the Clausius–Mossotti relation follows from the law of dispersion for such waves in the long wavelength approximation. The Maxwell Garnett mixing formula is then obtained from the Clausius–Mossotti relation if we assume that the particles are small dielectric spheres and use the appropriate expression for their polarizabilities. We also discuss briefly what happens beyond the long wavelength approximation and, in particular, give a few leads to the extended homogenization theories that include dynamic corrections and effects of nonlocality.

A drawback of the point-particle model is that it does not allow one to account for the shape and physical size of the inclusions in a mathematically consistent manner. In Section 5, we overcome this limitation by considering a general three-dimensional photonic crystal. A rigorous homogenization theory is obtained in the long wavelength limit. The resulting mixing formula is very similar to Maxwell Garnett’s but contains an additional term, which is responsible for the electromagnetic interaction of inclusions. This term also allows one to compute tensorial effective permittivities of anisotropic composites.

In Section 6, we show several numerical illustrations of the rigorous homogenization formula derived in Section 5. Here we consider various two-component composites consisting of the same materials (a conductor and a dielectric) as were used in the first part of this tutorial to illustrate the Wiener and Bergman–Milton bounds. The results are compared to the Maxwell Garnett’s and Bruggeman’s predictions, and it is shown that the Maxwell Garnett approximation is quite robust for isotropic composites, albeit we do not consider strong resonance effects that can render simple mixing formulas grossly inaccurate. However, in the case of anisotropic composites, Maxwell Garnett’s predictions are much less robust, even in the absence of strong resonances.

Consistently with the first part of this tutorial, we will use the Gaussian system of units.

2. MACROSCOPIC MAXWELL’S EQUATIONS IN CONTINUOUS MEDIA

The starting point for our exposition is the set of microscopic Maxwell’s equations for the electric and magnetic fields $\mathbf{E}$ and $\mathbf{B}$, respectively, and the electric current $\mathbf{J}$,

The fields $\mathbf{E}$ and $\mathbf{B}$ are fundamental in the electromagnetic theory because they appear in the expression for the Lorentz force $\mathbf{F}$ acting on a point charged particle, viz.,

Some part of the electric current $\mathbf{J}$ can be viewed as a known and externally controlled source. One example is a radiating antenna, which is fed in a prescribed way. We denote this term by ${\mathbf{J}}_{\text{ext}}$ and refer to it as to the external current. The electromagnetic fields created by ${\mathbf{J}}_{\text{ext}}$ will exert the force given by Eq. (3) and set the charges of any material object in motion. This will create an induced current ${\mathbf{J}}_{\mathrm{ind}}$ in that object. The total current everywhere in space is then given by

The induced current will create additional electromagnetic fields, which will, in turn, act back on the moving charges. Even if we can neglect the back-action of the induced fields on the external source, this interaction will result in a big unsolvable set of nonlinear equations. To complicate things further, it is known that any set of charged particles cannot hold itself together by the forces of the form given by Eq. (3) if we apply the Newtonian laws of motion. Therefore, stability of matter cannot be understood from classical first principles.

Naturally, one wishes to develop an electromagnetic theory of continuous media that avoids these complications. This goal is achieved by the macroscopic Maxwell’s equation. In a typical textbook exposition, these equations are obtained by averaging the microscopic Eq. (1) over a physically small volume. The derivation is based on the observation that $⟨\nabla \times \mathbf{E}⟩=\nabla \times ⟨\mathbf{E}⟩$ and similarly for other fields, where

Unfortunately, there exist several problems with the averaging approach. First, for the method to be mathematically sound, we need to make sure that the integrals of the form given by Eq. (5) exist for all quantities of interest and are independent of the shape of $\mathbb{D}$. However, the field of a truly point charge is not integrable and therefore it cannot be averaged. This problem is not rectified by introducing any number of point charges or by excluding small spherical regions around each charge from integration (this does not suppress the dependence on $\mathbb{D}$). The only solution to the divergence problem is to introduce a continuous density of charge. In particular, the quantum-mechanical DFT and the classical hydrodynamic model operate with continuous densities, which are already in some sense averages. However, these smooth functions are not obtained by actual averaging of some (generally, unknown) microscopic densities; they are simply postulated to exist and are computed according to some approximations.

Second, the microscopic quantities of interest are very difficult if not impossible to compute in continuous media. Of course, there exist some simple models such as the Lorentz or the Drude model. However, these models describe adequately only some special cases and, moreover, they do not employ averaging of the field of point charges. In the first part of this tutorial, we have averaged the field of elementary dipoles (in Section 2.B). This is a somewhat better-defined procedure albeit not without its own mathematical ambiguities, which we have encountered while considering anisotropy in Section 4.

We can conclude that the operation of microscopic field averaging is a purely formal procedure; it is not constructive in the sense that it does not yield any useful predictions. The averaging is simply postulated for pedagogical reasons and not used for anything. Instead, a phenomenological approach to constructing the macroscopic theory is implicitly used. Indeed, all we need to “derive” the macroscopic Maxwell’s equations are the following three postulates:

In the case of local and linear media, we can describe the functional dependence of Postulate (ii) by the equation

Now we need to make a few important comments.

We now substitute Eq. (6) into Eq. (4) and the result into the right-hand side of Eq. (1) and obtain the macroscopic equations

To complete the construction of the theory, we transform the macroscopic Maxwell’s equations to the frequency domain. Let us assume that the source ${\mathbf{J}}_{\text{ext}}(\mathbf{r},t)$ is strictly monochromatic. Then we can write

Equation (14) is the closed set of the macroscopic Maxwell’s equations. If the functions $\epsilon (\mathbf{r})$ and $\mu (\mathbf{r})$ are known and sufficiently “nice” mathematically, we can solve Eq. (14) uniquely. On the other hand, if we have a sample of material that is known to be homogeneous internally with some fixed $\epsilon$ and $\mu$, we can perform a set of experimental measurements with quasi-monochromatic radiation that will determine uniquely these quantities at the working frequency. For example, to determine the properties of a homogeneous isotropic slab, one can perform transmission and reflection measurements at several incidence angles, various states of polarization of the incident beam, and, perhaps, various widths of the slab. It should be kept in mind that the experiments of this kind are usually carried out with partially coherent sources of radiation and one should use the appropriate mathematical description of the polarization states of partially coherent radiation.

Finally, it is often more convenient to consider a given external (incident) field rather than a given external current. This is explained in Section 4.B.

3. PHYSICAL MEANING OF POLARIZATION AND MAGNETIZATION

It is almost universally accepted that the polarization and magnetization fields $\mathbf{P}$ and $\mathbf{M}$ are the differential volume densities of the electric and magnetic dipole moments in continuous media. However, in a typical textbook exposition of the subject, the mathematical arguments in favor of this belief are either missing or flawed. Sometimes it is implied that two functions are pointwise equal if their integrals over a given region are equal [see Eq. (16) below and its discussion]. Sometimes it is remarked that $\mathbf{P}$ and $\mathbf{M}$ cannot be determined uniquely from the induced current ${\mathbf{J}}_{\mathrm{ind}}$ and we must, therefore, apply some additional conditions to define these quantities unambiguously [this argument is discussed in Point (vii) of Section 2 above]. In yet other cases, it is argued that $\mathbf{P}$ and $\mathbf{M}$ must have some “physical meaning,” and then pictures are drawn of little polarized atoms or current loops. The problem with this approach is that it is relevant only to some very special cases.

The classical interpretation of $\mathbf{P}$ and $\mathbf{M}$, as we will call it, has gained wide acceptance for several reasons. First, it can be viewed as pedagogical since it gives a visual qualitative explanation of the polarization and magnetization phenomena. Second, the approach seemed to work fine when it could be applied in practice, that is, for molecular gases or simple composites. In other cases, which include pretty much all condensed matter and composites of complex geometry, there was not enough information to say anything. The microscopic charge density in the condensed phase could not be computed from simple models.

The situation changed in 1990s with the development of the density-functional theory and related computational methods. When the microscopic charge density was computed in some solids, it became obvious that the classical interpretation does not reproduce the experimental results with an acceptable accuracy. This motivated the development of the modern theory of polarization by King-Smith, Resta, and Vanderbilt [9–12]. The main idea of this new theory is that the polarization $\mathbf{P}$ itself (however defined) is not physically important; it is the change of $\mathbf{P}$ that is measurable and can be related to the macroscopic properties of matter. The point of view taken in this tutorial is very similar: the fundamental quantity of interest is the electric current—the time derivative of $\mathbf{P}$ (assuming $\mathbf{M}=0$). However, since our discussion is purely classical, we do not touch here the geometric quantum phase, which is central to the modern theory of polarization. We note that recent introductory-level reviews of this theory are available in Refs. [4,13].

Still, it is hard not to notice that the classical interpretation is surprisingly robust, at least for the electric polarization, and we naturally wish to understand its limits and conditions of applicability. Considering the importance of this subject, we will discuss it in detail in this section, primarily, in relation to the electric polarization (magnetization is also discussed briefly). In particular, we will see that the classical interpretation works fine for the Maxwell Garnett theory as long as the characteristic scale of inhomogeneities is much smaller than the wavelength. However, one should be conscious of the conditions of applicability of the classical interpretation, especially when high-frequency phenomena are considered.

A. Polarization and the Density of the Electric Dipole Moment

Consider a simple connected region $\mathbb{D}$ of volume $D$ and let the closed and regular surface $\mathbb{S}$ be its boundary. The region $\mathbb{D}$ is located arbitrarily with respect to a given sample of continuous medium, which occupies the spatial region $\mathbb{V}$. We wish to compute $\mathbf{d}(t)$—the macroscopic dipole moment of the matter contained in $\mathbb{D}$. To this end, we can start from Eq. (8) and write

If the surface integral in Eq. (15) was zero, we could conclude that the classical interpretation is correct. Indeed, we can take $\mathbb{D}$ to be sufficiently small so that ${\int }_{\mathbb{D}}\mathbf{P}(\mathbf{r},t)d^{3}r\approx \mathbf{P}({\mathbf{r}}_0)D$, where ${\mathbf{r}}_0$ is a point inside $\mathbb{D}$. Then $\mathbf{P}({\mathbf{r}}_0)\approx \mathbf{d}/D$ is, indeed, the dipole moment per unit volume. However, there is no reason to believe that the surface integral in Eq. (15) is zero. An obvious counterexample is a statically polarized piece of metal. Since the macroscopic charge can accumulate in conductors only on the surface (this result holds at arbitrary frequencies if the conductor is internally homogeneous), the dipole moment of any interior volume element is zero while the polarization $\mathbf{P}$ is not. In fact, the macroscopic dipole moment of any internal region of a homogeneous material is zero.

Thus, Eq. (15) written for the general case does not tell us much. To make a more useful statement, we can take $\mathbb{D}\supset \mathbb{V}$ so that the surface $\mathbb{S}$ is in vacuum and completely encloses the object. The surface integral in this case vanishes, and we obtain the following for the total dipole moment of the object:

However, we can use the result given by Eq. (16) constructively if we have a microscopic model or some “more fundamental” expression for the density of charge. As discussed above, this expression can come from one of the theoretical models such as the Lorentz model, the Drude model, the hydrodynamic model; it can result from an ab initio DFT calculation and so on. Let us say that somehow we have computed a fine-scale density of charge $\varrho$. We expect the following equality to hold:

We thus have made an identification of $\mathbf{P}$ as the density of microscopic dipole moment in a very special case. In particular, we have assumed that the geometry and frequency are such that the macroscopic polarization is spatially uniform inside the sample, the sample is dividable into completely identical cells, and the electromagnetic response of all cells is identical, including the cells at the boundary and in the interior of the sample.

All the conditions described above are important. For example, consider the physical situation sketched in Fig. 1. Here the charge distribution is perfectly periodic. However, if we forget about the left and right ends of the chain, there is no unique definition of the unit cell. We can define the unit cell so that its dipole moment takes any value in the interval $[-qh,qh]$. Therefore, if we want to make a meaningful choice, Condition (ii) must be used. Application of this condition depends on the left and right ends of the chain, which are not shown in the figure. If the medium is electrically neutral, there are only two possible terminations: either the left end is negative and the right end is positive or vice versa. In each case, we can break the medium completely into identical elementary cells. In one case, the macroscopic polarization of the medium $P_z$ (along the chain) is $+q/h^2$ and, in the other case, it is $-q/h^2$ (we assume that each cell is cubic); the intermediate values of $P_z$ are not possible.

 

Fig. 1. One-dimensional illustration of the ambiguity of defining the dipole moment of an element of volume. The idea of the figure is borrowed from Ref. [13].

Download Full Size | PDF

The above example suggests that what happens at the medium boundaries is really important. It is not possible to use a microscopic model of a medium without implicit or explicit consideration of the boundaries. Moreover, the familiar correspondence between $\mathbf{P}$ and the dipole moment of a cell can easily break down if the cells adjacent to the boundaries are in some way not equivalent to the cells in the interior.

Here is a simple one-dimensional example illustrating the importance of boundaries. Consider a microscopic charge distribution $\varrho (z)$ of a sample contained between the planes $z=0$ and $z=L=Nh$, which is given by the expression

 

Fig. 2. Illustration of the charge density given by Eq. (18) for $L=10h$. Here Lorentzians of the width $\gamma =0.01h$ are plotted instead of true delta functions.

Download Full Size | PDF

Let us now break the medium into the elementary cells $nh\le z<(n+1)h$, $n=0,1,\dots ,N-1$. It is easy to see that each elementary cell in the interior of the medium (that is, for $1\le n\le N-2$) acquires a dipole moment per unit surface ${\mathrm{d}}_{\mathbb{C}}={\varrho }_0(h^2/2\pi )\sin (2\pi \mathrm{\Delta }/h)$. According to the conventional wisdom, the medium should be assigned the macroscopic polarization

However, the conventional wisdom is wrong in this case. The problem is that the boundary cells are not equivalent to the cells in the interior. Of course, the boundary cells are not electrically neutral and, therefore, we cannot compute unambiguously their dipole moments. However, we can compute the total dipole moment of the whole structure, and it turns out to be equal to zero. It so happens that the dipole moment of the first and last cells (considered together) cancels exactly the dipole moment of all internal cells. If we assign the medium the dielectric permittivity $\epsilon =1+4\pi \chi$, we would get an obviously incorrect result. In fact, and somewhat paradoxically, the correct value of $\epsilon$ is 1; the medium is not macroscopically different from vacuum.

It is interesting to note that the total dipole moment of the sample remains zero if we replace the delta functions in Eq. (18) by Lorentzians of arbitrary width and unit integral weight.

Another way to understand the above result is to consider the microscopic current $j_z(z)$. Let us assume that $\mathrm{\Delta }$ is a function of time with the derivative $\dot{\mathrm{\Delta }}$. Then we can use the continuity equation $\partial \varrho /\partial t+\partial j_z/\partial z=0$ to compute the current. The result is $j_z=-{\varrho }_0\dot{\mathrm{\Delta }}\cos [2\pi (z+\mathrm{\Delta })/h]$. We have ${\int }_0^L{j}_z\mathrm{d}z=0$ and, moreover, the average of $j_z$ computed over an arbitrarily defined cell is also zero, viz., ${\int }_a^{a+h}{j}_z\mathrm{d}z=0$ for $0\le a\le L-h$.

We have just obtained an illustration of one important fact: unlike the dipole moment, the cell average of the microscopic current is independent of the choice of elementary cell. However, this is only natural: polarization and magnetization are introduced as auxiliary quantities while the current is fundamental. Analogously, various multipole moments that are used in the classical interpretation are auxiliary quantities whose introduction can be avoided altogether.

We now need to clarify that the example of Eq. (18) and Fig. 2 illustrates just one particular reason why the classical interpretation can be wrong: the possibility of charge leakage through the cell boundaries. This possibility is the main objection to the classical interpretation that is raised in the modern theory of polarization. We can, however, envisage a physical situation in which every elementary cell contains a localized distribution of charge that cannot leak through the cell boundaries. This situation was referred to as the Clausius–Mossotti picture in Ref. [4]. It is well understood now that the Clausius–Mossotti picture is not a realistic model for polarization of solids. Still, the model may seem to be not so bad for the problem of homogenization of macroscopic composite materials. For example, we can consider isolated metallic inclusions in a dielectric host or even in vacuum. In the latter case, the charge can definitely not leak through the cell boundaries.

However, there is another reason why the classical interpretation can fail. According to Condition (ii) above, we should be able to break a macroscopic sample into identical elementary cells. However, if the size of a unit cell is not very small compared to the free-space wavelength, such partitioning of the medium is not possible. The fields in two adjacent cells will never be identical in this case. Therefore, the assumption that ${\mathbf{d}}_{\mathrm{tot}}=N{\mathbf{d}}_{\mathbb{C}}$ will not hold. One can also say that an elementary cell is not a physically small volume in this case. For the problem of homogenization of macroscopic composites, this objection is much more important than the charge leakage.

We can now make the following conclusions:

B. Magnetization and the Density of the Magnetic Dipole Moment

Although we do not discuss magnetic effects in any detail in this tutorial, we will now look briefly at the vector of magnetization $\mathbf{M}$. The main point that we wish to make is that there is no direct analogy between $\mathbf{P}$ and $\mathbf{M}$ and the arguments used above for $\mathbf{P}$ cannot be applied to $\mathbf{M}$ without additional restrictions.

The classical interpretation of magnetization is even more problematic because an integral relation of the form of Eq. (16) does not generally hold for $\mathbf{M}$. Indeed, the macroscopic-induced magnetic moment $\mathbf{m}$ of the region $\mathbb{D}$ is by definition

Let us now take $\mathbb{D}\supset \mathbb{V}$ so that the surface integral vanishes. Then the total magnetic dipole moment of the object is

The above discussion concerns magnetization caused by the electric currents, which are often referred to as orbital, although they can also include the conductivity current. Of course, the orbital currents obey quantum laws at the microscopic level. However, macroscopically, they can be viewed as classical. Further details of the modern theory of orbital magnetization can be found in [14,15]. In contrast, ferromagnetism is a macroscopic quantum effect, which cannot be fully understood within the theoretical framework of classical physics. Nevertheless, the classical interpretation of $\mathbf{M}$ works fine for ferromagnetism. This is so, first, because the last term in Eq. (26) is negligible at the low frequencies for which the ferromagnetic effect is non-negligible and, second, because the elementary amperian current loops (the quantum-mechanical expectations of the electric current associated with an electron spin) fit well into the Clausius–Mossotti picture outlined above. In other words, the conditions of applicability of the classical interpretation that were discussed above for the case of electric polarization also hold for the elementary current loops in ferromagnetics.

4. WAVES ON LATTICES OF POLARIZED POINT PARTICLES

In the first part of this tutorial, we have derived the Maxwell Garnett approximation for a collection of point noninteracting dipoles at zero frequency. We now wish to construct a more general theory. Specifically, we will consider arbitrary frequencies and allow the dipoles to interact. We will, however, restrict our attention to a very simple geometry: point-like, spherically symmetric particles arranged on a cubic lattice. We will show that the Maxwell Garnett mixing formula follows from the dispersion relation for the electromagnetic waves propagating on such a lattice in the long wavelength limit. We will not use the classical interpretation of polarization as the dipole moment per unit volume for arriving at this conclusion. However, the classical interpretation is valid in the long wavelength limit and can be used to derive the same result.

The model of electromagnetically interacting point-like polarizable particles arranged on a three-dimensional lattice possesses an intuitive physical appeal. Historically, many authors have considered this model [16–21]. Most if not all of these references assume that the particles are embedded in vacuum. We will make an unfortunate compromise and work under the same assumption here, although it could have been more logical to start with the model of small inclusions in a host of permittivity ${\epsilon }_h\ne 1$. The main reason for our choice is simplicity. Introduction of a host medium is not conceptually difficult but will shift attention away from the key points, which we wish to avoid. However, we hope that anyone who has followed through the derivations of this section will be able to repeat them for a more general host. Alternatively, the end result of such a calculation can be obtained from the less general result of this section by the substitution ${\epsilon }_{\text{eff}}\to {\epsilon }_{\text{eff}}/{\epsilon }_h$, ${\epsilon }_i\to {\epsilon }_i/{\epsilon }_h$, where ${\epsilon }_{\text{eff}}$ is the effective permittivity of the structure and ${\epsilon }_i$ is the permittivity of inclusions. Toward the end of this section, we will apply this transformation to derive the Maxwell Garnett approximation for a general two-component mixture.

A. Polarizability of a Small Particle Beyond the Static Limit

Consider a small particle illuminated by an external monochromatic field of the complex amplitude $\mathbf{E}(\mathbf{r})$ and frequency $\omega$. If the particle size $a$ is much smaller than the wavelength of the incident radiation, it can be characterized by the dipole polarizability tensor $\widehat{\alpha }$, which is defined as the linear coefficient in the expression $\mathbf{d}=\widehat{\alpha }\mathbf{E}({\mathbf{r}}_0)$. Here $\mathbf{d}$ is the complex amplitude of the dipole moment, and ${\mathbf{r}}_0$ is a point inside the particle. If the particle has cubic symmetry, we can write $\widehat{\alpha }=\alpha \widehat{I}$, where $\alpha$ is a scalar and $\widehat{I}$ is the identity tensor. For simplicity, we will consider the scalar case in this section. However, the foregoing discussion also applies to all principal components of the tensor $\widehat{\alpha }$, assuming it can be diagonalized by an orthogonal transformation.

The polarizability is not arbitrary but satisfies certain constraints. A well-known result is that the extinction cross section of the particle (if excited by a plane wave) is ${\sigma }_e=4\pi k\mathrm{Im}\alpha$. Since ${\sigma }_e>0$, we conclude that $\mathrm{Im}\alpha >0$, but this is not the most interesting constraint. Another inequality can be obtained by noting that the scattering cross section under the same excitation conditions is ${\sigma }_s=(8\pi /3)k^4{|\alpha |}^2$. We know that ${\sigma }_e\ge {\sigma }_s$, where the equality holds for an ideally nonabsorbing particle. From this, we find that

The apparent contradiction can be resolved by noting that (FP3) is not quite accurate at finite frequencies: it does not account for the leading-order radiative correction (to the inverse polarizability), which is equal to $-2i{k}^3/3$. This correction can be relatively very small. However, it is needed to satisfy the optical theorem. When the electromagnetic interaction of several small particles is considered, not accounting for the radiative correction can yield strange results that contradict energy conservation [22,23].

It turns out that corrections to $\mathrm{Re}(1/\alpha )$ also exist and are even of lower order in $k$. We can refer to these corrections as to nonradiative since they are not associated with scattering. Unfortunately, these corrections are not uniquely defined. Different authors give different expressions for the nonradiative corrections [24]. The reason for this ambiguity is an irremovable singularity of the Green’s tensor for Maxwell’s equations; we will encounter this singularity momentarily. As a result, the nonradiative corrections depend not only on the particle shape (which could be expected) but also on the form of the incident field and on the reference point ${\mathbf{r}}_0$.

On the other hand, the first nonvanishing radiative correction is universal. This is so because the imaginary part of the Green’s tensor does not experience the singularity mentioned above. One can argue that it is the correction to $\mathrm{Im}(1/\alpha )$ that is physically important and should be retained, even in the limit $ka\to 0$. In accord with this argument, we will decompose the inverse polarizability as

One last question that we need to address is the following. If we account for the effects that are of the order of $O(k^2)$, is it correct to neglect all higher-order multipole moments? Strictly speaking, the answer is “no.” The magnetic dipole and electric quadrupole moments of small particles are quantities of the order of $O(k^2)$, the same as the difference $1/\tilde{\alpha }-1/{\alpha }_{\mathrm{LL}}$. We therefore will work in the regime when the nonradiative corrections and the higher-order multipole moments can be neglected.

As for the radiative correction, we will see that, in the particular problem we are solving, it exactly cancels. The cancellation takes place because a plane wave can propagate on an infinite, perfectly ordered lattice without scattering. Formally, the same result could be obtained if we simply forgot about the radiative correction and neglected the terms that are of the order of $O(k^2)$ and higher in all formulas. It may seem therefore that our discussion of the radiative correction was redundant, but this is not truly so. In problems that actually involve scattering, e.g., for finite or disordered sets of dipoles, accounting for the radiative correction can be crucial; its neglect can result in the appearance of nonphysical resonances or divergences and, ultimately, in violations of energy conservation. On the other hand, neglect of the nonradiative corrections does not lead to anything dramatic or unphysical. One can also say that $\mathrm{Re}(1/\alpha )$ and $\mathrm{Im}(1/\alpha )$ are responsible for different physical effects and should be considered separately.

B. Green’s Tensor Beyond the Static Limit

We next need an expression for the electric field radiated by a point dipole. This result is given by the free-space Green’s tensor for Maxwell’s equations. We have already encountered the Green’s tensor in the static limit in the first part of this tutorial [see (FP8)]. Now we need to generalize this result to finite frequencies.

We start from the frequency-domain macroscopic Maxwell’s equations (14). We consider a nonmagnetic medium with $\mu =1$ and let $\epsilon (\mathbf{r})=1+4\pi \chi (\mathbf{r})$, where the susceptibility $\chi (\mathbf{r})$ is zero in vacuum. Upon substitution of $\mathbf{H}$ from the second into the first equation of (14), and by using the definitions in Eqs. (10) and (12) to write $\chi (\mathbf{r})\mathbf{E}(\mathbf{r})=\mathbf{P}(\mathbf{r})$, we obtain

We will refer to the solution of Eq. (31a), ${\mathbf{E}}_{\text{ext}}$, as to the external or incident field. Here ${\mathbf{E}}_{\text{ext}}$ is the field that would have existed everywhere in space in the absence of the sample. Note that

In contrast, the second term in Eq. (30) is the field that is generated or scattered by the sample; we denote this term by ${\mathbf{E}}_{\text{scatt}}$ so that

The physical significance of the Green’s tensor can be understood as follows. Consider scattering of an external field by a small particle occupying the region $\mathbb{V}$ whose characteristic size $a$ is much smaller than the free-space wavelength so that $ka\ll 1$. Let the point ${\mathbf{r}}_0\in \mathbb{V}$ be the “center” of $\mathbb{V}$. Then, for $|\mathbf{r}-{\mathbf{r}}_0|\gg a$, we can write

We will now use Eq. (31b) to compute $\widehat{G}$. Accounting for the translational invariance of free space, we can write

The real-space Green’s tensor can be obtained by computing the Fourier integral in Eq. (35). To this end, we can replace the tensor $\mathbf{p}\otimes \mathbf{p}$ by ${\nabla }_{\mathbf{r}}\otimes {\nabla }_{{\mathbf{r}}^{\prime }}$. Then we can perform the angular integration by using the elementary integration rules and the ingratiation over $p$ by residues. We can use the symmetry of the integrand to extend integration to the whole real axis; also, we should keep in mind that $k$ has an infinitesimal imaginary part according to Eq. (14c). The result of this calculation is

Let us now evaluate the integral of ${\widehat{F}}_R(\mathbf{r})$ over a small ball of radius $a$. We can use the identity ${\int }_{4\pi }(\mathbf{r}\otimes \mathbf{r})\mathrm{d}\mathrm{\Omega }=\frac{4\pi }{3}{r}^2\widehat{I}$, where $d\mathrm{\Omega }$ is the element of solid angle, to compute the angular integrals. We then obtain

On the other hand, if we evaluate the integral in Eq. (44) over a more general nonspherical region, the term $O(r^{-3})$ would not integrate to zero. We will observe a physical manifestation of this fact when we consider waves on lattices of point polarizable particles. We have also seen a similar difficulty when we considered anisotropy in the first part of this tutorial (in Section 4). All we can say now is that the “regularizing” decomposition in Eq. (40) is appropriate when spherical particles and regions are considered.

We can also define the “regular” part of the Green’s tensor in the spatial Fourier domain. That is, we can write

C. Coupled-Dipole Equations on an Infinite Lattice

Consider an assembly of identical, point-like particles of scalar polarizability $\alpha$ located at the positions ${\mathbf{r}}_n$. Let the whole system be illuminated by the external field ${\mathbf{E}}_{\text{ext}}(\mathbf{r})$. Then the dipole moments of the particles ${\mathbf{d}}_n$ are coupled to the external field and to each other by the coupled-dipole equation

Unfortunately, our exposition of the coupled-dipole equation is by necessity brief. We note only that it is important to use the correct definitions for all quantities appearing in Eq. (47). In particular, we should use the polarizability $\alpha$ as it was defined in Section 4.A above and not the “bare polarizability” [25] or the quasi-static polarizability ${\alpha }_{\mathrm{LL}}$. Also, we have used the “regular” part of the Green’s tensor, ${\widehat{G}}_R$, in Eq. (47) rather than $\widehat{G}$. In the real-space representation, this is unimportant as long as the $m=n$ term is not included in the summation. However, upon spatial Fourier transformation, the difference will show up. Additional details about the coupled-dipole equations can be found in Refs. [19,24,26].

Our goal is to find an eigensolution to Eq. (47) on an infinite lattice and to analyze the corresponding dispersion relation. To this end, we assume that ${\mathbf{r}}_n$ are the vertices of an infinite cubic lattice of pitch $h$, set the external field ${\mathbf{E}}_{\text{ext}}$ to zero and seek the solution to Eq. (47) in the form ${\mathbf{d}}_n=\mathbf{d}\exp (i\mathbf{q}·{\mathbf{r}}_n)$. Upon substitution of this ansatz into Eq. (47), we obtain the equation

Lattice sums are ubiquitous in physics [27] but not always convergent. Correspondingly, a lot of effort has been devoted to computing and regularizing diverging lattice sums. Two approaches to this task are common: allowing the point particles to have higher-order multipole moments or allowing the particles to have a finite size and shape. Ultimately, the two approaches are equivalent since they describe the same physical reality. We will resort to the second approach—that is, we will assume that the particles are small spheres.

We now proceed with the calculation. First, we use the Fourier representation in Eq. (35) (applied to ${\widehat{G}}_R$) to rewrite Eq. (49) as

As expected, the expression in Eq. (51) is divergent: the last term in the right-hand side is just ${\widehat{F}}_R(0)$, but ${\widehat{F}}_R(0)$ is not well defined. Thus, we will use the regularizing substitution

We now put everything together. We substitute Eq. (53) into Eq. (51) and Eq. (51) into Eq. (48), use the decomposition in Eq. (28) for $1/\alpha$, and obtain

D. Long Wavelength Approximation

The dispersion relation given by Eq. (54) is still too complicated for analysis. To proceed, we will need to make the long wavelength approximation. Namely, we will assume that $kh\ll 1$ and $|q_p|h\ll 1$, where $p=x,y,z$. We should keep in mind that, in our model, $a\ll h$ so that $ka\ll 1$ also holds. Then, with the precision of the approximation used, we can write $1/\tilde{\alpha }+k^2/a\approx 1/{\alpha }_{\mathrm{LL}}$, where ${\alpha }_{\mathrm{LL}}$ is the quasi-static polarizability of a particle.

We thus have simplified the left-hand side of Eq. (54). In the right-hand side, we will write the term with ${\mathbf{g}}_n=0$ and the rest of the sum separately. We will thus obtain

We now use the above fact and the expression in Eq. (46) for ${\widehat{\mathbf{K}}}_R(\mathbf{q})$ to obtain the following equation:

The only thing left for us to do is to consider the polarization states. The effective permittivity of a cubic lattice of small spherical particles cannot be anisotropic (from symmetry, the principal axes of the tensor ${\widehat{\epsilon }}_{\text{eff}}$ must coincide with the crystallographic axes, and the directions along all three crystallographic axes are equivalent). Therefore, we must assume that the medium is characterizable by a scalar effective permittivity ${\epsilon }_{\text{eff}}$ and $q^2=k^2{\epsilon }_{\text{eff}}$ (we have assumed here that ${\mu }_{\text{eff}}=1$, more on this later). Under the same conditions, the medium can support transverse waves for which $\mathbf{q}·\mathbf{d}=0$. Then it follows from Eq. (59) that the dispersion equation is

Equation (60) is the inverted Clausius–Mossotti relation [compare to (FP16)]. We therefore have derived the Clausius–Mossotti relation as the long wavelength limit of the dispersion relation for electromagnetic waves on a cubic lattice of small spherical particles.

The Maxwell Garnett approximation can be obtained if we use a specific expression for ${\alpha }_{\mathrm{LL}}$. Since we have assumed that the particles are small spheres, we will use (FP3). Substituting everything into Eq. (60) and inverting the equation, we obtain

E. Beyond the Long Wavelength Approximation

We now discuss briefly what can happen beyond the long wavelength approximation. The material of this subsection is not a self-contained exposition but a collection of a few basic ideas and leads for further reading and exploration.

Recall that we have neglected the nonradiative corrections to the inverse polarizability and replaced the term $1/\tilde{\alpha }+k^2/a$ by $1/{\alpha }_{\mathrm{LL}}$. We cannot make a more precise statement about the inverse polarizability without abandoning the model of point-like dipoles. However, we can still use this model and introduce physically meaningful corrections to the long wavelength approximation if $a\ll h$. In this case, the inequality $ka\ll 1$ can hold quite strongly, but the inequalities $kh\ll 1$ and $|q_p|h\ll 1$ cannot be so strong. We then expand the term $K_R(\mathbf{q}+{\mathbf{g}}_n)$ with ${\mathbf{g}}_n\ne 0$ [in Eq. (55)] in powers of $k$ and $\mathbf{q}$. The long wavelength result (56) will give the zeroth order in this expansion. Some of the expansion terms will sum to zero due to symmetry but others will remain. If we make all expressions dimensionless by multiplying the numerator and denominator in the expression for $K_R(\mathbf{q}+{\mathbf{g}}_n)$ by $h$, the corrections will be of the order of $O({(kh)}^n)$ and $O({(q_{p}h)}^n)$, where $n$ is even. This is true for all systems with a center of symmetry.

The terms proportional to ${(kh)}^2$, ${(kh)}^4$, etc., are the dynamic corrections to the dispersion equation. Accounting for these terms does not change the algebraic structure of the equation; it remains quadratic in $\mathbf{q}$ and its solutions remain similar to those in local, isotropic, homogeneous media. Generally, these solutions can be written as $q^2={\epsilon }_{\text{eff}}(k)k^2$. However, if we account for the dynamic corrections, ${\epsilon }_{\text{eff}}(k)$ will no longer satisfy the scaling laws discussed in Section 7 of [1]. In particular, rescaling the composite as $\mathbf{r}\to \beta \mathbf{r}$ will change ${\epsilon }_{\text{eff}}(k)$; the law of unaltered ratios will also not hold.

Keeping the terms proportional to ${(q_{p}h)}^2$, ${(q_{p}h)}^4$, or $q_x^2{q}_y^2{q}_z^2{h}^6$ (anisotropic terms of this form do not contradict the cubic symmetry), etc., will fundamentally change the algebraic form of the dispersion equation: it will become much more complicated than the quadratic equation characteristic of local and isotropic media. The dispersion equation can be written in this case as $q^2={\epsilon }_{\text{eff}}(k,\mathbf{q})k^2$. The explicit and sometimes very complicated dependence of ${\epsilon }_{\text{eff}}$ on $\mathbf{q}$ can be considered as an effect of spatial nonlocality of the constitutive relations.

Many important physical effects of nonlocality can be glimpsed from the dispersion equation. This includes additional waves, birefringence, and optical activity (in systems without a center of symmetry) [6]. Unfortunately, the dispersion equation does not give us the complete information about the medium. Indeed, solving boundary value problems in finite samples requires the knowledge of the influence function $f_e(\mathbf{r},{\mathbf{r}}^{\prime };\tau )$ [the nonlocal generalization of the function $f_e(\mathbf{r},\tau )$ that appears in Eq. (7a)]. However, $f_e(\mathbf{r},{\mathbf{r}}^{\prime };\tau )$ is not the spatial Fourier transform of ${\epsilon }_{\text{eff}}(k,\mathbf{q})$ nor is it a function of the shift $\mathbf{r}-{\mathbf{r}}^{\prime }$; the latter simplification is obtained only sufficiently far from the medium boundaries. Generally, spatial nonlocality and boundary value problems do not come together easily.

Things are further complicated by a widespread belief that weak nonlocality of the dielectric permittivity is physically indistinguishable from nontrivial magnetic response. For the precise statement of this equivalence principle and relevant references, see [28]. In fact, the equivalence holds with respect to the dispersion relation but not the medium impedance. A detailed analysis of the equivalence principle is given in Ref. [29]. Here we only mention an obvious counterexample: an analytical solution can be easily obtained for a plane wave incident on a slab or a sphere with purely local $\epsilon$ and $\mu$. However, the same problem cannot be solved so easily if the object has a nonlocal permittivity and $\mu =1$. If we define some function $f_e(\mathbf{r},{\mathbf{r}}^{\prime };\tau )$ and solve somehow the resultant integro-differential equations, the electromagnetic fields scattered by the object will be not the same as in the first case.

The above discussion underscores the fundamental importance of what happens at the medium boundary. In Section 3.A, we have already seen that disregard of the boundaries can yield an incorrect result, even in very simple physical situations. In the more complicated scenarios considered by the extended homogenization theories, boundaries are crucial.

To illustrate the importance of boundaries one more time, consider an elementary but frequently overlooked example. If we know only the dispersion relation of an isotropic and local medium, there is no good way to define the parameters $\epsilon$ and $\mu$ separately. Indeed, the law of dispersion in this medium is $q^2=n^2{k}^2$, where $n^2=\epsilon \mu$. Assume that we have computed $n^2$. All we can gain from this knowledge is the product of $\epsilon$ and $\mu$, not these two quantities separately. This fact is frequently overlooked when extended homogenization theories are developed to predict both $\epsilon$ and $\mu$ but infinite unbounded media are considered. The main difficulty of the extended homogenization theories is accounting for the lack of translational invariance that occurs near the medium boundaries.

Now we have to backtrack a little and recall that, in the previous subsection, we have assumed without a justification that ${\mu }_{\text{eff}}=1$. This was, however, not a random guess. In the long wavelength limit, one can consider the medium with a planar boundary and compute the reflection and transmission coefficients [30]. Alternatively, one can compute the total magnetic moment of the medium in some simple geometry and apply the results of Section 3.B. Both approaches will indicate that the correct choice is ${\mu }_{\text{eff}}=1$. Also, all experimental evidence tells us that there is no artificial magnetism in composites whose heterogeneity size is much smaller than the wavelength.

Beyond the long wavelength approximation, one enters a gray area in which the scaling laws discussed in Section 7 of [1] do not work, where the magnetic effects can be considered as important and the introduction of a nontrivial effective magnetic permeability is feasible, even in intrinsically nonmagnetic composites. The big question is whether the medium can be reasonably characterized by local tensors $\widehat{\epsilon }$ and $\widehat{\mu }$ or, perhaps, by a more general but still local constitutive tensor. The extended (nonasymptotic) electromagnetic homogenization theories are subjects of active ongoing research [31–38]. The above references and many references therein will give an interested reader useful leads for further exploration of the subject.

5. WAVES IN PHOTONIC CRYSTALS IN THE LOW FREQUENCY LIMIT

We have seen in the previous section that the model of an infinite lattice of point polarizable particles is not well-behaved mathematically unless we assume implicitly or explicitly something about the particle shape and size. All methods of regularizing the divergent lattice sums boil down to making these assumptions, often in a disguised form. It makes sense therefore to abandon the point-particle model and consider the actual physical inclusions arranged periodically. We will see in this section that this approach is ultimately rewarding since it is free from any mathematical difficulties and allows one to account self-consistently for the effects of inclusions interaction, anisotropy, etc. After all, Maxwell’s equations are not self-contradictory and all mathematical difficulties grow from making ill-justified or unphysical approximations.

The mathematical developments shown in this section are not fundamentally different from those used for the case of point polarizable particles. However, some more advanced algebraic manipulations will be used.

A. Bloch Waves in a Photonic Crystal

Consider a three-dimensional periodic medium shown schematically in Fig. 3. Inclusions of arbitrary shape are arranged periodically on an infinite cubic lattice of pitch $h$. We denote the elementary cells of the lattice by ${\mathbb{C}}_n$ and the center of ${\mathbb{C}}_n$ by ${\mathbf{r}}_n=h(n_x,n_y,n_z)$. As above, $n=(n_x,n_y,n_z)$ is a composite index.

 

Fig. 3. A schematic illustration of a two-component photonic crystal. Identical inclusions are periodically arranged in a host medium. For simplicity, we assume that the host medium is vacuum. The pattern is three-dimensional.

Download Full Size | PDF

Each elementary cell contains an inclusion. Let the spatial region occupied by the $n$th inclusion be ${\mathbb{V}}_n\in {\mathbb{C}}_n$. All regions ${\mathbb{V}}_n$ are of exactly the same shape and the same volume $V$ and differ only by translations. In the figure, we have shown ${\mathbb{V}}_n$ as connected, but this assumption is not important for us. Also, ${\mathbb{V}}_n$ can touch the boundary of ${\mathbb{C}}_n$ and the center of the cell, ${\mathbf{r}}_n$, can lie inside or outside of ${\mathbb{V}}_n$.

Further, we assume that the material of inclusions has the permittivity $\epsilon \ne 1$ and the material of the host is vacuum. Later on, we will generalize the consideration to the case of a two-component mixture by the usual substitution $\epsilon \to {\epsilon }_i/{\epsilon }_h$ and ${\epsilon }_{\text{eff}}\to {\epsilon }_{\text{eff}}/{\epsilon }_h$.