On the Theory of Light Propagation in Crystalline Dielectrics

A synoptic view on the long-established theory of light propagation in crystalline dielectrics is presented, providing a new exact solution for the microscopic local electromagnetic field thus disclosing the role of the divergence-free (transversal) and curl-free (longitudinal) parts of the electromagnetic field inside a material as a function of the density of polarizable atoms. Our results enable fast and efficient calculation of the photonic bandstructure and also the (non-local) dielectric tensor, solely with the crystalline symmetry and atom-individual polarizabilities as input.

small sized 3M × 3M matrix eigenvalue problem, with M denoting the number of polarizable atoms (ions) in the unit cell. Identifying the macroscopic electric field inside the sample with the spatially low-pass filtered microscopic local electric field, the dielectric 3 × 3-tensor εΛ (q, ω) of crystal optics emerges, relating the accordingly low-pass filtered microscopic polarization field to the macroscopic electric field, solely with the individual microscopic polarizabilities α (R, ω) of atoms (molecules, ions) at a site R and the crystalline symmetry as input into the theory. Decomposing the microscopic local electric field into longitudinal and transversal parts, an effective wave equation determining the radiative part of the macroscopic field in terms of the transverse dielectric tensor ε (T ) ab (q, ω) is deduced from the exact solution to the field-integral equations. The Taylor expansion ε (T ) ab (q, ω) = ε (T ) ab (ω) + i c γ abc (ω) qc + c,d α abcd (ω) qcq d + ... around q = 0 provides then insight into various optical phenomena connected to retardation and non-locality of the dielectric tensor, in full agreement with the phenomenological reasoning of Agranovich and Ginzburg in "Crystal Optics with Spatial Dispersion, and Excitons" (Springer Berlin Heidelberg, 1984): the eigenvalues of the tensor ε (T ) ab (ω) describing chromatic dispersion of the index of refraction and birefringence, the first order term γ abc (ω) specifying rotary power (natural optical activity), the second order term α abcd (ω) shaping the effects of a spatial-dispersion-induced birefringence, a critical parameter for the design of lenses made from CaF2 and BaF2 for optical lithography systems in the ultraviolet. In the static limit an exact expression for εΛ is deduced, that conforms with general thermodynamic stability criteria and reduces for cubic symmetry to the Clausius-Mossotti relation. Considering various dielectric crystals comprising atoms with known polarizabilities from the literature, in all cases the calculated indices of refraction, the rotary power and the spatial-dispersion-induced birefringence coincide well with the experimental data, thus illustrating the utility of the theory. For ionic crystals, exemplarily for CsI and RbCl, a satisfactory agreement between theory and the measured chromatic dispersion of the index of refraction is shown over a wide frequency interval, ranging from ultraviolet to far infra-red, accomplishing this with an appreciably smaller number of adjustable parameters compared to the well known Sellmeier fit.

I. INTRODUCTION
When optical signals traverse a transparent dielectric, for example a fused quartz (silica) prism, the light travels at different speed depending on frequency f = ω 2π , so that the shape of a wave packet, say composed of mixed frequencies |f − f c | ≤ ∆fc 2 around a carrier frequency f c in a frequency interval of width ∆f c , tends to spread out. This is the well known chromatic dispersion effect resulting from the frequency dependence of the refractive index n = n (ω). Microscopic considerations based on first principles reveal, that the frequency dependence of the refractive index n (ω) is directly connected to the retarded response of the polarizable constituents of matter, the latter distinguishing themselves as atoms, molecules or ions. The chromatic dispersion effect is further supplemented by the effects of crystalline anisotropy and also by the effects of spatial dispersion brought about by the non local dependence of this response, that is a charge at point r recollects the action exerted on it at another position r ′ [1, 2].
Both fundamental features of the electromagnetic response, retardation and non locality, can be described jointly by a dielectric (tensor) function ε Λ (q, ω) depending on (circular) frequency ω and on wave vector q. While the dependence of the dielectric function on ω explains the chromatic dispersion effect, and its anisotropy for q = 0 describes birefringence, the dependence on wave vector q is directly connected to phenomena like natural optical activity and spatial dispersion induced birefringence. Furthermore, under the influence of a static magnetic induction field B (0) , respectively a static electric field E (0) , the additional dependence of the dielectric function on those static fields gives rise to many other magneto-optic and electro-optic phenomena, for instance the Faraday effect, the Pockels effect and also the Kerr effect [3].
When electromagnetic waves propagate inside a dielectric material, the microscopic local electric field convoyed by those waves exerts a small supplemental force onto charge carriers inside the atoms (ions, molecules) comprising that material, pulling apart inside each atom the positively charged nucleus and the negatively charged bound electrons by a (small) shift in opposite direction. Of course the position of the total mass of a charge neutral atom, considering valence electrons and nucleus together, remains then unchanged. Since the atomic nucleus is certainly much heavier than the bound electrons, the resulting shift of the barycenter of the bound electrons by far surpasses the associated tiny shift of the position of the atomic nucleus. This is the effect of induced electronic polarization.
A fully microscopic theory of the optical properties of a material certainly requires to consider the positive charged atomic nuclei and the electrons, taking into account the laws of quantum statistical physics and low energy quantum electrodynamics, for example [1,4]. In the ensuing discussion we shall accept a phenomenological (semiclassical) picture of the matter light interaction based on the fundamental fact, that matter is stable, i.e. the energy of an atom located at a site R is, in the absence of external fields, a minimum against any (small) displacement of its bound electrons. Accordingly, in reaction to the presence of a (weak) local time dependent electric field E (R, t) the electrons bound to an atom redistribute themselves, so that considered from outside, an atom comprising a number Z of electrons, acquires an induced electric dipole moment. This is the basic idea of the phenomenological classical model of atomic polarizability due to Lorentz, who described the induced dipole moment of electrons, tied to a heavy (immobile) nucleus by a harmonic spring, solving a harmonic oscillator problem with that local electric field acting as a drive. Incidentally, the Lorentz model predicts a frequency dependent Fourier amplitude [5] of the induced electronic dipole moment of an atom at site R that coincides with a full quantum mechanical calculation, see supplemental material [6]: a, a ′ ∈ {x, y, z} Here α a,a ′ (R, ω) denotes the atom-individual electric polarizability given by the summation index ν running here over all eigenstates |Φ ν (R) of the multi-electron configration of that atom except the groundstate |Φ 0 (R) . Indeed, this expression looks like it was derived from an ensemble of (classical) harmonic oscillators with resonance frequencies ω ν (R), each oscillator being driven by the local electric field with a factor of proportionality f a,a ′ (ν; R) > 0, the so called oscillator strength. However, the physical meaning of the oscillator strength is only revealed by quantum mechanics. Withd a = − |e| Z n=1r (n) a denoting the dipole operator of the count Z of electrons tied to an atom at position R, then according to quantum mechanics f a,a ′ (ν; R) = 1 |e| 2 2m ω ν (R) Φ 0 (R) |d a |Φ ν (R) Φ ν (R) |d a ′ |Φ 0 (R) , explaining why the oscillator strength is a measure for how much a bound electron contributes to the electric polarizability of an atom, say under a transition from the multi-electron ground state |Φ 0 (R) of the atom Hamiltonian with eigenvalue Ω 0 (R) to an excited multi-electron eigenstate |Φ ν (R) with eigenvalue Ω ν (R). The differences ω ν (R) ≡ Ω ν (R) − Ω 0 (R) in (2) are the optical transition frequencies, the life-time parameter τ ν (R) > 0 describes spontaneous emission as reasoned by quantum electrodynamics and thus being always present if the atom was excited to an eigenstate state ν = 0. Expanded details how the result (2) can be derived within first order time dependent perturbation theory in response to a weak time dependent electric field, including a discussion of the f -sum rule, see supplemental material [6]. De facto, the spectrum of atoms with Z > 1 cannot be calculated ab initio with sufficient precision, the exact taking into account of electronic correlations being (alas) an unsolved problem. In what follows we therefore conceive the optical transition frequencies ω ν R (j) , the life-time parameter τ ν R (j) and the oscillator strengths f a,a ′ ν, R (j) of each atom species positioned at a site R (j) = R + η (j) , with R ∈ Λ a lattice vector and η (j) indicating a position of an atom (ion, molecule) inside a unit cell C Λ of the crystal Λ, as fitting parameters, so that the optical properties as calculated from the dielectric tensor ε Λ (q, ω) of the crystal coincide with experiment. How this objective can be accomplished, and in particular how ε Λ (q, ω) depends on the individual atom polarizabilities α a,a ′ η (j) , ω , and thus via (2) on the atom individual multi-electron spectrum, we shall elaborate in what follows.
A time dependent external electromagnetic signal E ext (r, t) incident upon a material probe polarizes the atoms inside, and for weak amplitude of the incident field the induced polarization at the atom sites R will then be proportional to that amplitude. However the microscopic local electric field E (r, t) polarizing an atom (ion, molecule) positioned at site R is not known a priory, because all atoms in the sample will get polarized and hence all act back with a (retarded) reaction in response to the primary external field E ext (r, t). Then everywhere inside (and outside) of the probe there holds the secondary induced electric field E ind (r, t) being a superposition of the individually radiated and retarded electric fields emitted by all the atoms inside the sample, that have been polarized in turn by that field E (r, t). With a suitable model of the polarizability of individual atoms thus an integral equation determining the microscopic local electric field E (r, t) emerges. So far, everything said is well known from the original (early) literature [7][8][9][10][11] and from highly cited textbooks on crystal optics, for example [3, [12][13][14]. For a concise summary of the pioneering works on crystal optics of Ewald and v. Laue (and later authors) see [15]. Nevertheless, we believe the approach we present in what follows truely discerns from traditional presentations of the subject. For example, nowhere do we make use of the Ewald-Oseen extinction theorem to recover the correct index of refraction n of a material, it being customary practice in the so called rigorous theory of dispersion [12,13] to consider the wave incident from free space to be propagated with vacuum light velocity c and the signal induced in the sample to be propagated with velocity c/n, the Lorentz-Lorenz formula connecting the index of refraction n with the polarizability of individual atoms thus emerging as a solvability condition for the field-integral equation stated (implicitely) in (4).

Outline
In Sec. II we first establish on the basis of the fundamental Maxwell equations an exact integral equation for the microscopic local electric field E (r, t) in (4) with an explicit formula for the integral kernel that derives directly from the atom polarizability (2). Posing boundary conditions for the components of the electromagnetic field at the boundary of a material probe is then redundant. Moreover, the frequency dependence of longitudinal and transversal parts of the electromagnetic field are treated consistently, thus making everywhere in our calculations the correct static limit ω → 0 accessible.
In Sec. III we solve the field-integral equation for crystalline dielectric materials exactly making use of a set of non-standard Bloch functions, not constructed from plane waves but designed from eigenfunctions of the position operator, representing a complete orthonormal system of eigenfunctions of the translation operator T R under a shift by a lattice vector R ∈ Λ. Accordingly, instead of expanding the kernel of the field-integral equation in the well known basis of plane waves e i(q+G)r , thus requiring to handle for each wave vector q in the Brillouin zone C Λ −1 of a lattice Λ then (infinite) matrices labelled by reciprocal lattice vectors G, G ′ ∈ Λ −1 , our choice of eigenfunctions of the translation operator T R sidesteps the inversion (and truncation) of such large matrices, thus easing notably the determination of the photonic bandstructure of a crystal. Also we show, if the incident electromagnetic wave was purely transversal, yet the microscopic local electric field features both, a transversal and a longitudinal component, see Fig.(4), the strength of the longitudinal component being strongly dependent on the density of polarizable atoms (ions, molecules) in the crystal.
Thereafter we present, exemplifying our calculation method, results for the photonic bandstructure of diamond (M = 2), that have been calculated previously with other (computationally more time-consuming) methods. Our findings for the photonic bandstructures for various monoatomic Bravais lattices (M = 1) we present and discuss in the supplemental material [6]. In comparison to well known (phenomenological) work on photonic bandstructures [16,17], within which the (macroscopic) Maxwell's equations in a superlattice are solved assuming a spatially repetetive varying index of refraction, it turns out that in our approach based on the field-integral equations the need to eliminate unphysical "longitudinal" modes doesn't arise.
In section IV we introduce the notion of a macroscopic electric field E (r, t) inside a material, conceiving it with regard to spatial variations as a low pass filtered signal, with the solution to the field-integral equation (33) as input. Relating the macroscopic polarization to that macroscopic electric field, thereafter the macroscopic dielectric 3 × 3−tensor ε Λ (q, ω) of a crystalline dielectric material emerges, with chromatic and spatial dispersion fully taken into account. The exact expression for ε Λ (q, ω) thus obtained is solely dependent on the symmetry of the lattice Λ under consideration and on the polarizability α R (j) , ω of individual atoms (ions, molecules) positioned at their equilibrium sites R (j) in the crystal. We also confirm, exemplarily for the ionic crystal CsI, that our formula for the dielectric tensor with regard to its frequency dependence indeed obeys to the Lyddane-Sachs-Teller relation. With a view to the key role of locality claimed by macroscopic electrodynamics [18] we caution the reader not to discard the q-dependence of the dielectric tensor ε Λ (q, ω). As emphasized by Ginzburg and Agranovich [3], a rich variety of optical effects, including rotary power and spatial dispersion induced birefringence, manifestly proof the importance of the non local nature of the dielectric response.
We next derive directly from the microscopic field-integral equations a set of differential equations that determine the spatial variation of the transversal and longitudinal parts of the macroscopic electric field E (r, t), without prior knowledge of the microscopic local field E (r, t).
If the external field was purely transversal, this set of coupled differential equations reduces to a wave equation determining the radiative part of the macroscopic field. In this way the parts of the dielectric tensor ε Λ (q, ω) are indentified that determine the transversal dielectric tensor ε (T ) (q, ω) comprising the optical properties of a dielectric crystal. Assuming spatial dispersion is a weak effect, which assumption applies for many transparent media, the Taylor expansion of ε (T ) ab (q, ω) around q = 0 then provides (yet in an implicit manner) access to various optical phenomena featuring the propagation of light in dielectric crystals [3], for instance chromatic dispersion and birefringence, rotary power (natural optical activity) and also the (weak) effects of a spatial-dispersion-induced birefringence, the latter being a critical problem for the design of lens elements made from crystalline materials like CaF 2 and BaF 2 widely used in optical lithograpy systems in the ultraviolet [19,20]. Further we summarize in section IV, see Table I, Fig.10 and also Fig.11, to what large extend our theory of the dielectric tensor for crystalline dielectrics agrees with measurements over a wide range of optical frequencies for a series of well known crystalline materials, including for example Bi 12 TiO 20 and also Bi 12 SiO 20 , both crystals featuring a large number of basis atoms (M = 66) in the unit cell, thus demonstrating the utility of our approach.
If on the other hand the external field was purely longitudinal, the differential equations derived from the fieldintegral equations reduce to a Poisson type equation for a scalar potential function determining the macroscopic electric field, thus identifying the parts of the dielectric tensor ε Λ (q, ω) featuring electric-field screening, like in electrostatics. Also we deduce an exact analytic formula for the static dielectric tensor ε Λ , that conforms with general thermodynamic stability criteria [21] and applies for all 14 monoatomic Bravais lattices (M = 1). For the special case of cubic symmetry the long-known Clausius-Mossotti relation is recovered.

II. THE FIELD-INTEGRAL EQUATIONS
Consider a fixed inertial frame, with a polarizable dielectric material at rest occupying a volume Ω in that frame. Without loss of generality let the microscopic charge density inside Ω be given the representation with P (r, t) denoting the vector field of electric polarization [1]. We find it then convenient to decompose the associated microscopic current density j (r, t) = ∂ t P (r, t) flowing inside Ω into longitudinal and transversal parts. To serve this purpose we introduce integral kernels with labels a, a ′ , a ′′ ∈ {x, y, z} specifying Cartesian components, and Denoting the convolution of two kernels A aa ′ (r, r ′ ) and B aa ′ (r, r ′ ) with the symbol and writing for the kernel representing unity then the validity of all the relations distinctive for projection operators are readily confirmed: In position space, for |r| = 0, the kernel of the longitudinal and the transverse projection operator both correspond to a dipole field, see [22] for a concise derivation: It follows from what has been said that the longitudinal part (L) , respectively the transversal part (T) of the current distribution j a (r, t) conforms with the convolution integrals It should be underlined, this link between the original vector field j(r, t) and the associated longitudinal (transversal) part j (L,T ) (r, t) is non local. There holds by construction According to the Helmholtz-Hodge theorem such a decomposition of a vector field j (r, t) is unique. A recent thourough discussion and compilation of the literature on the subject can be found in [23]. Because Maxwell's equations are linear, particular monochromatic solutions a (r, ω) +j (L) a (r, ω) = 0 , the second group involves the Cartesian components of the Fourier amplitudes of the transversal electromagnetic field, which obey to the following six inhomogenous scalar Helmholtz equations, the respective source terms being provided by the Fourier amplitudes of the transversal current density, with c = 1 √ ε0µ0 the speed of light in free space: We are interested to solve (16) and (17) considering now a geometry consisting of two disjoint material bodies at rest, see Fig.1, i.e. Ω = Ω S ∪ Ω P and Ω S ∩ Ω P = ∅.
Let us refer to the body Ω S as the source, and to the body Ω P as the probe. Accordingly, we split the current densityj a (r, ω) into an externally controlled partial currentj ext,a (r, ω) flowing solely inside Ω S , and into an induced partial currentj ind,a (r, ω) flowing solely inside the probe Ω P : If the source Ω S was a cavity producing a laser beam, and if the surface of the probe Ω P would be reflecting (parts) of the radiation incident back into that cavity, certainly there would exist a backaction from the probe to the source, leading then to the existence of an induced partial current flowing also inside the source region Ω S . We preclude here and in the following any such backaction effects. Provided transfer of charge between Ω S and Ω P is prohibited, there holds charge conservation separately (!) inside Ω S and inside Ω P . Consequently, inside the domain Ω S there holds iωρ ext (r, ω) = ∇ ·j ext (r, ω).
Of courseρ ext (r, ω) ≡ 0 andj ext,a (r, ω) ≡ 0 for r / ∈ Ω S . It follows from (16) that the longitudinal electric fieldẼ (L) a (r, ω) is already determined by the longitudinal part of the current distribution (18), a seemingly simple result. Indeed from (16) and (18) theñ ext,a (r, ω) + where the external longitudinal field is derived from a scalar potentialφ like in electrostatics. Even thoughj ext,a (r, ω) was restricted to be non vanishing solely inside the source domain Ω S , its longitudinally projected partj ext,a (r, ω) also exists outside of Ω S , the non locality of that projection thus becoming manifest.
The electromagnetic fields radiated by the transversal current densityj (17) can be readily determined introducing the retarded Green's functiong Im (ω) → 0 + , withg (r − r ′ , ω) the solution of the three-dimensional inhomogeneous scalar Helmholtz equation in free space R 3 with a point source (of strength unity) at position r ′ : As the distance |r − r ′ | to the source position at r ′ increases theng (r − r ′ , ω) behaves like an outgoing spherical wave. Based on Green's identity applied to the domain Ω = Ω S ∪ Ω P , and again applied to the complementary domain R 3 \ Ω, there follow directly from (17) and (24) for all points r / ∈ Ω S the integral representations andB a (r, ω) =B ext,a (r, ω) + µ 0ˆd For details of the derivation of (25) and (26), see supplemental material [6]. According to (6) the Fourier amplitude of the induced current density flowing inside the probe volume Ω P is directly proportional to the Fourier amplitude of the microscopic electric polarization: Combining now the respective longitudinal and transverse parts by adding the integral representations (20) and (25) there followsẼ withG denoting the electromagnetic kernel. The Fourier transformed kernelG aa ′ (r, ω) in the wave vector domain we denote asḠ = Assuming a small amplitude of the perturbing external fieldẼ ext,a (r, ω) there results inside the probe at a position r ∈ Ω P the (total) microscopic polarizatioñ Within a fully microscopic approach the response kernelχ ext,aa ′ (r, r ′ , ω) is calculated from Kubo's formula in reaction to the presence of the external fieldẼ ext,a ′ (r ′ , ω) . However, what we are really interested in here, is not the response kernelχ ext,aa ′ (r, r ′ , ω) connecting the polarizationP a (r, ω) inside the probe with the external fieldẼ ext,a ′ (r ′ , ω), but the dielectric susceptibility kernelχ aa ′ (r, r ′ , ω) connectingP a (r, ω) with the microscopic local electric fieldẼ a ′ (r ′ , ω) that acts on each atom (ion, molecule): As emphasized by Keldysh [1], there is no need to consider dielectric and magnetic susceptibilities separately, the latter being already incorporated in the non locality of the dielectric kernel. The fundamental field-integral equation determining the microscopic local electric field is then given by an inhomogenous integral equation comprising the dielectric kernelχ: The link between the dielectric susceptibility kernelχ and the response kernelχ ext is readily identified [1], combining (31) with (32):χ For a more detailed explanation, see supplemental material [6]. So, ifχ ext was known, say from a full microscopic calculation with Kubo's formula, thenχ follows from (34). In case the probe volume Ω P was of finite size, solving the integral equation (34) for the microscopic dielectric susceptibility kernelχ poses a formidable (numerical) problem, the non-locality radius ofχ ext being substantially enhanced up to the macroscopic scale by the presence of a boundary [1]. In the next section III we circumvent this problem and introduce a phenomenological model forχ that proves a posteriori to be appropriate to describe in detail the propagation of light in dielectric crystals.

III. MICROSCOPIC LOCAL ELECTRIC FIELD IN CRYSTALLINE DIELECTRICS
Crystalline order assumes each equilibrium position of atoms (ions, molecules) inside a material corresponds to a site vector R (j) = R + η (j) , where R denotes a lattice vector in a Bravais lattice Λ, and η (j) with j ∈ {1, 2, ..., M } is a set of position vectors indicating equilibrium positions of the atoms (molecules, ions) inside the unit cell C Λ . The entire probe volume Ω P can be thought of to be filled translating a number |Λ P | of Wigner-Seitz unit cells C Λ by lattice vectors R ∈ Λ. This subset of all such Bravais lattice vectors R ∈ Λ we denote as Λ P ⊂ Λ. Under the action of the external electric fieldẼ ext,a ′ (r ′ , ω) each atom (ion, molecule) gets polarized proportional to the local electric fieldẼ a R (j) , ω acting at position R (j) of that atom. Then a simple phenomenological model for the microscopic dielectric susceptibility kernel in crystalline insulators emerges assembling first all individual atom contributions located at positions η (j) inside the unit cell positioned around the origin R = 0, and then sum over all such cells filling the probe volume Ω P of the crystal lattice Λ under consideration: Essentially, this susceptibility kernel describes a lattice periodic arrangement of point dipoles [10,11,14]. The diagonal terms j = j ′ refer to the afore mentioned effects of induced electronic polarization of single atoms (ions, molecules) at the positions η (j) inside the unit cell C Λ around the lattice vector R = 0. Based on the functional form of the frequency dependence of the individual microscopic polarizability (2) of a single atom (ion, molecule), in actual fact being equivalent to a Lorentz-oscillator model, we write now a simple phenomenological ansatz with fitting parameters α (j) 0 and ω (j) 0 (and possibly also including a small life-time parameter τ (j) = 1/γ (j) > 0 representing spontaneous emission, see supplemental material [6]): the atom polarizability is well approximated by its static value α (j) 0 . For most frequencies, except if ω approaches a transition frequency ω (j) 0 near to an absorption band, we may assume γ (j) → 0 + . In case there are two or more basis atoms in the Wigner-Seitz cell qualitatively new effects need to be considered. Off diagonal terms j = j ′ exist for M ≥ 2 and designate mutual influences of atoms (ions, molecules) positioned at different sites η (j) and η (j ′ ) inside the unit cell. On one hand, in reaction to the presence of a propagating electromagnetic wave the overlap integral(s) determining the sharing of electron pairs between neighbouring atoms undergo (slight) changes. On the other hand, in crystals like N aCl, CsI, RbCl etc., the formation of ions needs to be taken into account. Besides the effect of induced electronic polarization of single atoms (ions), attributed to a shift of the barycenter of the electrons bound to individual atoms (ions) under the action of the local field on-site R (j) , an additional shift of the position of the positive ions relative to the position of the negative ions concurs, thus leading to ionic displacement polarizability. This effect is typically noticeable in the electromagnetic response to radiation with frequency ω of order of characteristic lattice vibration frequencies ω ph , de facto being mainly of concern for radiation in the infrared.
Conceiving also the off-diagonal polarizabilities j = j ′ not as input from microscopic theory, but in the phenomenological guise of a Lorentz-oscillator model for oppositely charged ion pairs [24], with fitting parameters α (j,j ′ ) 0 , ω (j,j ′ ) 0 and small damping γ (j,j ′ ) , we find the well known chromatic dispersion of the index of refraction n (ω) is nicely reproduced by our theory of the dielectric tensor ε Λ (q, ω) for a variety of dielectric crystals, see Table I in section IV. While the electronic polarizabilities (37) of single atoms are often well approximated by their static value, the ionic displacement polarizabilities (38) have a characteristic dependence on frequency in the infrared. Note, that we refrain adding the contributions of the electronic polarizabilities of single atoms to the contributions of the ionic displacement polarizabilities, there being no justification for this [24].
For example, in the case of CsI with M = 2 ions in the unit cell, our phenomenological approach in (35) brings into altogether 6 parameters. These are the two transition frequencies ω . With the model (35) we then successfully reproduced experimental data for the chromatic dispersion of the refraction index n (ω), exemplarily for CsI and RbCl, over a wide frequency interval ranging from ultraviolet to infrared, utilizing only these 6 parameters instead of 17 parameters as required by a Sellmeier fit, see Fig.10 and Table II. Restricting to optical frequencies well above the infrared range (but always well below atomic excitation energies), then mostly the electrons bound around individual ions (atoms) will react to the electromagnetic fields. In this case the effect of induced electric polarization is predominant and the effects of ionic polarisation, as represented by the off diagonal contributions j = j ′ in (36), can be considered as small, so that It should be noted that retaining in (39) the possibility of non zero off diagonal Cartesian terms a = a ′ in the Lorentzoscillator model of atomic polarizabilitiesα aa ′ η (j) , ω then enables the study of crystals composed of anisotropic polarizable ionic or molecular subunits in the elementary cell. For instance, the study of the influence of an external static magnetic induction field B 0 or static electric field E 0 on the propagation of light in crystals also requires to retain non zero off diagonal Cartesian components a = a ′ in (39). As important examples we mention the magneto-optical Faraday effect and the electro-optical Pockels effect [3].
Obviously, under a translation by a lattice vector R ∈ Λ the dielectric kernel (35) remains invariant: The dielectric kernel considered at fixed position r ∈ Ω P as a function of r ′ − r will typically undergo already discernible variations traversing a short route of order of the interatomic distance a. In what follows we show, that the microscopic local electric field amplitudeẼ (r, ω) then also displays discernible spatial variations on that same short length scale a, even though the primary incident light signalẼ ext (r, ω) was discernibly varying only on the much longer length scale set by the wavelength λ ≫ a of light in free space.

Non-Standard Bloch Functions
Because of the periodicity (40) it is an obvious choice to expand the microscopic local electric field in a complete basis of Bloch eigenstates of the translation operatorT R shifting the argument of any function f (r) according tô with R ∈ Λ a Bravais lattice vector. Routinely in problems with a lattice periodic potential an expansion provided by the orthonormal and complete basis system of plane waves constructed from eigenfunctions of the momentum operator is deployed Here Λ −1 denotes the reciprocal lattice conjugate to the lattice Λ and C −1 Λ is the Brillouin zone. Making use of e iG·R = 1 for R ∈ Λ and any reciprocal lattice vector G ∈ Λ −1 one readily confirmŝ The eigenvalues e ik·R associated with the eigenfunctions e i(k+G)·r of the translation operatorT R being highly degenerate, any function of the form (Fourier coefficients denoted as u G ), is an eigenfunction of the translation operator(s)T R : Based on a rigorous theory of the microscopic electromagnetic response kernel, Dolgov and Maximov [2] presented for crystalline systems a profound analysis of the dielectric function, representing the microscopic kernels in (34) as (infinite) matrices with respect to the basis system (42). Adversely, the integral kernelG •χ in (33), when represented in the basis of plane waves e i(k+G)·r k∈C −1 Λ ,G∈Λ −1 , turns out to be a full up matrix G + k| G •χ aa ′ |G ′ + k , labeled by a wavevector k ∈ C Λ −1 and an infinite set of reciprocal lattice vectors G, G ′ ∈ Λ −1 . In the (numerical) calculations thus the handling of large matrices is required.
An alternative set of eigenfunctions of the translation operator(s)T R , so that the dielectric kernel when represented in the new basis appears as a sparse matrix, is highly desirable. Observing, that any lattice periodic function u (r) can be generated from a fragment u (0) (s) that equals to u (s) inside the Wigner-Seitz cell C Λ and is zero outside, we may as well represent any such fragment u (0) (s) ≡ s|u (0) as a linear combination with the complete and orthonormal set of eigenstates {|s } s∈CΛ of the position operatorr a obeying to the well known relationsr So instead of expanding into the basis system of plane waves e i(k+G)·r k∈C −1 Λ ,G∈Λ −1 there emerges as an alternative expanding into the following system of eigenfunctions of the translation operatorT R : By construction:T The set of states {w (r; s, k)} k∈C Λ −1 ,s∈CΛ is for one thing labeled by a wavevector k ∈ C Λ −1 , that ranges through the count |Λ P | of values of wave vectors k in the first Brillouin zone C Λ −1 , consistent with the Born-von Karman periodic boundary conditions, and for another thing it is labeled by position vectors s ∈ C Λ that range within the Wigner-Seitz cell C Λ of the lattice Λ. The system {w (r; s, k)} k∈C Λ −1 , s∈CΛ spans a complete and orthonormal basis system of eigenfunctions of the translation operator(s)T R , for a proof see supplemental material [6].

Solution of the Field-Integral Equations for a Dielectric Crystal
Representing now the microscopic local electric field in the complete and orthonormal basis system {w (r; s, k)} k∈C Λ −1 , s∈CΛ we writeẼ Conversely, the expansion coefficients representing that field arẽ the field-integral equation (33) is transformed into an equivalent integral equation determining the expansion coefficientsẽ a (s, k, ω): The matrix elements of the dielectric kernel in the basis (50) are readily evaluated, see supplemental material [6]: Here we introduced notation such that andζ Λ denotes a 3 × 3 matrix of lattice sums formed with the electromagnetic kernel (29): For any Bravais lattice vector R ∈ Λ there holds The definition of the lattice sums (58) by cases is to ensure no atom can get polarized by its self-generated electromagnetic field. It is important to realize that Instead there holds A fast and precise numerical method for the calculation of the lattice sums ζ Λ (s, k, ω) and ζ (0) Λ (k, ω) we disclose in our supplemental material [6].
Insertion of (56) gives at once the field amplitudesẽ a (s, k, ω) in terms of a finite number of amplitudesẽ .., M } counting the positions s = η (j) of the polarizable atoms (ions, molecules) inside the Wigner-Seitz cell C Λ and a ∈ {x, y, z} denoting Cartesian components: Taking subsequently for j = 1, 2, ..., M the limit s → η (j) , then from (62) a 3M × 3M system of linear equations determining those amplitudesẽ ext,a (k, ω) of the external field is obtained:ẽ This result clearly brings out the advantage of the basis system {w (r; s, k)} k∈C Λ −1 ,s∈CΛ over the conventional basis system of plane waves e i(k+G)·r k∈C −1 Λ ,G∈Λ −1 , thus dispensing the consideration of large full up dielectric matrix kernels G + k| G •χ aa ′ |G ′ + k labeled by an infinite number of reciprocal lattice vectors G, G ′ ∈ Λ −1 . As an incidental remark see [25].
The determination of the expansion coefficientsẽ a (s, k, ω) representing the local electric field amplitude via (62) requires solving (63), a small sized 3M × 3M linear system of equations. Introducing the 3M × 3M matrices and the explicit solution of (63) reads Finally, substituting the expansion coefficients (66) into (62) we then obtain from (52) the following explicit representation for the microscopic local electric field: Incidentally, the general result (67) for the microscopic local electric field inside a crystal comprises also the case of multiple beams of incident signals, all with frequency ω but possibly with different wavevectors k, in this case represented by a multitude of expansion coefficientsẽ (j ′ ) ext,a ′ (k, ω) carrying different wave vectors k, for example see [14].
It should be emphasized that the well known dispersion relation ω = c |q| of light in vacuum expresses the solvability condition for the external field (69) solving the homogeneous Maxwell equations in free space. However, in reality external signals represent solutions of the inhomogeneous Maxwell equations with a source currentj ext,a (r, ω) = q ′jext,a (q ′ , ω) e iq ′ ·r flowing inside a source volume Ω S , the latter being usually positioned at a large (but finite) distance to the probe volume Ω P , see Eq. (57) in supplemental material [6]. So a typical external signalẼ ext,a (r, ω) at fixed (circular) frequency ω is composed of a bunch of mixed wave vectors q' forming that signal: This feature enables to regard ω and q from now on as independent variables. Making use of the superposition principle it suffices to specialize to a single beam. For convenience we consider now an external signal in the guise of a plane wave, with (circular) frequency ω and wavevector q so that Assuming then q ∈ C Λ −1 , certainly not a strong limitation for optical signals propagating in crystalline materials, it follows at once from the defining equation (54) that the expansion coefficients of the external signal are independent on the value of the mode index s ∈ C Λ . All the more the coefficientsẽ (j ′ ) ext,a ′ (k, ω) then being independent on η (j) there holdsẽ ext,a (s, k, ω) = |Λ P |δ k,qĒext,a (q, ω) ≡ẽ (j ′ ) ext,a ′ (k, ω) .
Exploiting next the lattice periodicity of ζ Λ (s, q, ω), see (59), there holds for s = 0 the Fourier series representation withζ Λ (G, q, ω) denoting the Fourier coefficients, see supplemental material [6]: Insertion of the Fourier series representation (72) into (71) leads finally to the following result for the microscopic local electric field: a,a ′′ (q + G, ω) K Λ (G, q, ω) a ′′ ,a ′′′Ēext,a ′′′ (q, ω) . Here we introduced the kernel a quantity being directly connected to the microscopic atom-individual polarizabilities, see (57). Eq. (74) constitutes the central result of this section. It explicitely determines the microscopic local electric field inside a crystal, that is the electric field that polarizes atoms (ions, molecules) at their positions R (j) = R + η (j) inside a unit cell of the lattice in reaction to an incident plane wave (69). While the external fieldẼ ext,a (r, ω) displays spatial variations on the length scale set by the wavelength λ of the incident optical signal in free space, the contributions of the reciprocal lattice vectors G = 0 to the Fourier series in (74) make it manifest that the microscopic fieldẼ a (r, ω) displays on the back of the slowly varying envelope e iq·r rapid spatial variations on the (atomic) scale set by the lattice constant a Λ ≪ λ, see Fig. 2.
The result (71) also reveals that the strength of the microscopic local field amplitude inside the crystal strongly depends on the choice of the propagation directionq = q |q| , a feature being directly connected to the photonic band structure implicitely encoded in the eigenvalues of the matrix Γ (q, ω). The huge size of the induced electric field strength, represented by the differenceẼ a (r, ω)−Ẽ ext,a (r, ω), can indeed be prorated to the predominant longitudinal character of the microscopic local electric field (71) inside a probe volume with a high density of polarizable atoms (ions, molecules). This is intuitively accessible in view of the quasi static electric field inside a material probe at a point r ∈ Ω P originating from nearby positioned induced atomic dipoles.
To elucidate the nature of the microscopic local electric field in (71) and (74) respectively, let us decompose it into longitudinal and transversal parts. With (74) we write now a,a ′ (q + G)Ḡ a ′ ,a ′′ (q + G, ω) K Λ (G, q, ω) a ′′ ,a ′′′Ēext,a ′′′ (q, ω) . Taking into account (30) we have The searched-for longitudinal and transversal parts of the microscopic local electric field are thus explicitely determined:Ẽ The spatial variations of the transversal amplitudeẼ  of the eigenvalue problem The dispersion relation of photons, i.e. the photonic bandstructure ω n (q), can now be readily determined for any number M of basis atoms inside the unit cell C Λ by first solving (numerically) the eigenvalue problem (82) for a given wave vector q ∈ C Λ −1 as a function of ω, thus obtaining a family of 3M eigenvalue curves γ n (q, ω) vs. ω, and then solving (numerically) for n = 1, 2, 3, ...3M the implicit equations (81) for the unknown ω so that Varying then the wavevector q along (widely) different symmetry lines inside the Brillouin zone C Λ −1 various pieces of the photonic bandstructure ω n (q) emerge, as is exemplarily displayed for the diamond lattice (M = 2) in Fig.  6(a). Like electrons moving in a periodic potential also electromagnetic waves propagating in a crystal are governed by a bandstructure, for example [16,17,26]. While the wave function for electrons (discarding spin-orbit forces and Zeeman splitting ) is a scalar, propagating electromagnetic waves are vectorfields. Incident lightsignals composed of frequencies within an omni-directional band gap will be reflected from such a crystal irrespective of the light source being polarized or unpolarized, which is interesting for technical applications, for instance dielectric mirrors, filters or antenna-substrates. Further examples of photonic band structures calculated from (82) and (83) with our method, exemplarily for sc-, fcc-and bcc-lattices (M = 1), are given in supplemental material [6]. Our results comply well with calculations carried out within the frame of the Fano-Hopfield model [27,28] in the case of weak coupling of the oscillators. Earlier calculations for super-lattices on the basis of the macroscopic Maxwell equations [16,17,26], carried out in a basis of plane waves e i(q+G)r , require for q fixed a large number N of reciprocal lattice vectors G to be taken into account, for accurate computations typically N ≥ 2000 [29], thus leading to a huge 3N × 3N matrix eigenvalue problem for the modes and mode freqencies. Unfortunately, these calculations suffer from an inconsistency, as a number N of spurious zero eigenvalue modes needs explicit elimination by an ad hoc transversality constraint. Because in our approach the extension of the scatterers is tiny compared to the lattice constant, our results cease agreement with theirs in certain details of the photonic bandstructure at high photon frequency, while for optical frequencies and below our results are in full agreement with theirs. Regarding the calculational cost of our approach: in the case of sc-, fcc-and bcc-lattices with one atom in the unit cell we solve (for each q-vector) a 3 × 3 eigenvalue problem, and in the case of diamond with two atoms in the unit cell then a 6 × 6 eigenvalue problem. Parameters of interest to propagation of light pulses, for instance group velocity and effective photon mass, can be determined conveniently using k · p-perturbation theory [29], a method often employed in solid state electronic band structure theory [24].
Finally, it should be considered that higher bands n > 1 of the photonic band structure ω n (q) in real crystalline materials are not credibly calculable. While the concept of a photonic band structure exhibiting many band branches ω n (q) certainly applies for (artificial) superlattice structures with large mesoscopic lattice constant a Λ ≫ a, it can be accepted only with reserve for a real dielectric material. This is because for a microscopic lattice constant a Λ ≃ a photon frequencies above ω ≥ ω Λ = c × π aΛ are far and beyond the ultra-violet. In this case the effects of radiation damping, as represented by the imaginary part of the lattice sums, can no longer be neglected. For a derivation of (84) see supplemental material [6]. In Fig.6(b) the effect of taking into account radiation damping is exemplarily shown for the diamond lattice (M = 2). The corresponding plots revealing the influence of radiation damping for the sc-, fcc-and bcc-lattices (M = 1) we also present in [6]. Of course, if the wavelength of the external electromagnetic field is ultra short then the point dipole ansatz (35) for the dielectric susceptibility should be extended to include also magnetic dipoles and electric quadrupoles on equal footing [30].

IV. THE DIELECTRIC TENSOR OF MACROSCOPIC ELECTRODYNAMICS
The macroscopic electric fieldẼ a (r, ω) inside the probe we conceive as a low-pass filter applied to the Fourier-transformationĒ a (q, ω) of the spatially rapidly varying microscopic local electric fieldẼ a (r, ω). Introducing a cut-off wavenumber q c so that q < q c implies q ∈ C Λ −1 , theñ E a (q, ω) = 1 |Ω P |ˆΩ P d 3 re −iq·rẼ a (r, ω) .
Let us abbreviate for G = 0: On the other hand there holds keeping in mind the restriction q ∈ C Λ −1 : Restricting to the long wavelength limit |q| < q c and thus identifyingP a (q, ω) =P a (q, ω) andĒ a (q, ω) =Ē a (q, ω), and then combining (95) (97) for CsI with microscopic polarizabilities with the parameters from Table II. The ionic polarizability with frequency dependence given by (38) also includes a small damping parameter γ > 0. ≡ εΛ ωT · 10 2 = 3.05 above and beyond all optical phonon frequencies: Insisting both lines should be identical for any external field amplitudeĒ ext,a (q, ω) immediately leads (with help of elementary matrix algebra) to the identification This is a central result. The macroscopic dielectric tensor [ε Λ (q, ω)] aa ′ is solely determined by the lattice sums ζ Λ (s, q, ω) of the Bravais lattice Λ under consideration and the individual polarizations α a ′′ ,a ′ η (j ′′ ) , η (j ′ ) , ω of the atoms (ions) inside the unit cell. As a test of the analytic structure of the dielectric function ε Λ (q = 0, ω) in the complex frequency domain we checked the Lyddane-Sachs-Teller relation, see Fig. 7. In agreement with general considerations under ω → −ω the real part of ε Λ (q, ω) is an even function of ω and the imaginary part is an odd one.

Macroscopic Electric Field
It follows from what has beeen said that the macroscopic electric field amplitude is determined directly from the Fourier series representation (74) discarding all contributions of reciprocal lattice vectors G = 0: In position space then the longitudinal (L) and transversal (T) macroscopic electric field amplitudes read Comparing now with (79) and (80) we see at once that Accordingly the microscopic local electric field and the macroscopic electric field differ by the contributions of the sums over all reciprocal lattice vectors G = 0: In Fig. 8 we compare the spatial variation of the transversal macroscopic electric field amplitude (100) with the spatial variation of the transversal microscopic local electric field (80) along a path as shown in Fig.3, assuming the external electric field was purely transversal, i.e.Ẽ

Macroscopic Magnetic Induction Field
The amplitude of the microscopic local magnetic induction fieldB (r, ω) is of course directly connected to the amplitude of the local microscopic electric fieldẼ (r, ω) via Faraday's law: Insertion of the representation (103) for the microscopic local electric field amplitude,Ẽ a (r, ω) =Ẽ a (r, ω)+δẼ a (r, ω), leads in this way immediately tõ Identifying now the macroscopic magnetic induction field amplitude viã where the correction term representing the difference to the microscopic magnetic induction field amplitude is Like in the electric field case, the macroscopic magnetic induction field amplitudeB c (r, ω) represents the low pass filtered microscopic local magnetic induction field amplitudeB c (r, ω). The plot of the residue δB y (r, ω) =B y (r, ω)− B y (r, ω) is displayed in Fig.9. Clearly,B y (r, ω) andB y (r, ω) essentially coincide. Note that in the electric field case the relative size of the residue δẼ (T ) z (r, ω) along the same path turned out to be even smaller, see Fig.8.

Deducing the Differential Equations of Macroscopic Electrodynamics
Restricting to long wavelengths so that |q| < q c let us rewrite (95) in the guisē

Multiplication on both sides with the inverse of the kernel (30) gives
In terms of the projection operators (8) then so that (110) leads to a,a ′ (q) Ē ext,a ′ (q, ω) . ext,a (q, ω) and the external charge distribution̺ ext (q, ω): So the right hand side in (112) reduces together with the Fourier transformed relation (21) to Consequently (112) assumes the guise If the dependence on wavevector q of the dielectric tensor in (118) can be ignored, we replace ε Λ (q, ω) → ε Λ (ω) and obtain then in position space the well known (so called) vector wave equation determining the macroscopic electric field: It should be noted, that hereẼ (r, ω) still may be decomposed into divergence-free (transversal) and curl-free (longitudinal) parts,Ẽ (r, ω) =Ẽ (T ) (r, ω) +Ẽ (L) (r, ω). Thus it is deceptive to interpret (119) as a wave equation determining electromagnetic radiation as propagating photons with speed determined by the eigenvalues of the dielectric tensor ε Λ (ω), unlessẼ (L) (r, ω) ≡ 0. To find the differential equations for the transversal and longitudinal parts of the macroscopic field amplitude let us first introduce block matrix notation specifying transversal and longitudinal projections of the dielectric tensor: A, B ∈ {L, T } Then the vector wave equation (118) separates into two coupled equations for the respective transversal and longitudinal Fourier amplitudes of the macroscopic field: ext,a (q, ω) Choosing Eq. (119) as a starting point for the transport theory of radiation (light intensity) inside a (possibly disordered) material appears according to what has been said questionable, as the fluctuation contributionẼ a (r, ω) − E a (r, ω) ≡ δẼ a (r, ω), see Eq.(103), is in this case not included, despite the product δẼ a (r, ω)δẼ b (r ′ , ω) apparently comprising a spatially slowly varying interference contribution.

Wave Equation with Renormalized Speed of Light
If the external field was purely transversal, i.e.Ē (L) ext,a (q, ω) ≡ 0, then the longitudinal component of the macroscopic field is readily eliminated in (121), provided the inverse of the longitudinal block ε (L,L) Λ (q, ω) exists: Insertion leads to ext,a (q, ω) .
If the dependence on wavevector q of the (transversal) dielectric function can be ignored, then ε ab (ω), and equation (125) corresponds in position space to a (scalar) wave equation determining the Cartesian components of the transversal macroscopic electric field amplitudeẼ (T ) (r, ω) propagating inside a dielectric crystal, in full agreement with the standard theory of the propagation of polarized light in transparent dielectric materials: With a choice of a coordinate frame such that the dielectric tensor is diagonal in that frame, ε (T ) aa ′ (ω) = δ aa ′ n 2 a (ω), one finds for light propagating along a high symmetry axis corresponding to an eigenvector of that dielectric tensor the usual reduction of the speed of light, c → c/n a (ω) characteristic for in general birefringent crystalline dielectrics. Note that because of chromatic dispersion of the index of refraction then those eigenvectors may undergo a corresponding chromatic dispersion of axes as well, yet this effect being observable only for monoclinic and triclinic crystalline symmetry [12].
If the dependence on wavevector q of the (transversal) dielectric function cannot be ignored, the Taylor expansion around q = 0 in (125) provides then insight into various optical phenomena connected to retardation and non-locality of the dielectric tensor, in full agreement with the phenomenological reasoning of Agranovich and Ginzburg [3]: the eigenvalues of the symmetric tensor ε ba (ω) describing chromatic dispersion of the index of refraction and birefringence, the antisymmetric first order term γ abc (ω) specifying rotary power (natural optical activity), the second order terms α abcd (ω) shaping the intrinsic effects of a spatial-dispersion-induced-birefringence. Of course, in a centrosymmetric crystal there exists no natural optical activity: γ abc (ω) ≡ 0. The tensor α abcd (ω) originates from the symmetry-breaking evoked by the finite q-vector of the light [3], it displays in general 3 × 3 × 3 × 3 = 81 components. But for crystals with cubic symmetry the number of independent components of that tensor reduces substantially: for symmetry group T and T h there exist four independent components, for symmetry group T d ,O h and O there exist three independent components and for isotropic systems that number reduces to two [3]. The (weak) effects of a dispersion-induced-birefringence indeed give as a matter of fact reason for concern regarding the image quality of dielectric lenses made from CaF 2 and BaF 2 , a topic of prime importance designing modern lithographic optical systems in the ultraviolet [19,20].
The dielectric tensor ε Λ (q, ω) being a functional of the microscopic polarizabilities α η (j) , ω of atoms (ions, molecules) located at position η (j) in the unit cell C Λ , for example see (2), undergoes variations in proportion to changes of those atom individual polarizabilities caused by external static fields. For instance, in the presence of a static external magnetic induction field B 0 , the polarizability of a single atom displays now a magnetic field induced anisotropy [32], even though the atom-individual polarizability in zero field α η (j) , ω was isotropic: Such changes then reflect in corresponding changes of the dielectric tensor. As a result the Taylor expansion of the transversal dielectric tensor of the crystal now reads to leading order in the small quantities q and B 0 : Note the symmetry Correspondingly, if the dependence of the transversal dielectric tensor on wavevector q can be ignored, then light propagation in the presence of a static field B 0 inside a dielectric crystal is governed by the wave equation (126), but with a dielectric tensor now dependent on that static magnetic field: The presence of the antisymmetric tensor β abc (ω) = c ′ ǫ abc ′ λ c ′ c (ω) in (131) and with that said in (126), now leads to left and right circularly polarized waves propagating at slightly different speeds, thus giving rise to a magnetic field induced circular birefringence. If B 0 ∧ q = 0, i.e. if B 0 is orientated parallel or anti-parallel to the directionq of light propagation, this is the well known Faraday rotation effect of a light waves linear polarization in a static magnetic field.

Electric-Field Screening
Conversely, if the external current source was purely longitudinal, i.e.j with the longitudinal dielectric tensor describing electric-field screening. From equation (132) we readily infer For optical frequencies and below it is (obviously) adequate to approximate ε (L) (q, ω) ≃ ε (L,L) (q, ω), provided det |q| 2 I − ω 2 c 2 ε (T,T ) (q, ω) = 0. Then In position space, and representing the longitudinal electric fieldẼ (L) (r, ω) = −∇φ (r, ω), this conforms in the long wavelength limit |q| → 0 to the usual Poisson type equation of electrostatics determining the scalar potentialφ (r, ω) of a given charge distribution̺ ext (r, ω): Now it is manifest that it is the dielectric 3 × 3 tensor that describes electric-field screening, like in electrostatics.

Results of Calculations for the Chromatic Dispersion, Natural Optical Activity and Spatial-Dispersion-Induced Birefringence in Various Dielectric Crystals
Conceiving the polarizability of atoms or ions not as given input from microscopic theory, but as a fit function of the type (36), the diagonal parts depending on two parameters ω for each pair (j, j ′ ) of ions j, j ′ ∈ {1, 2, ...M } in the unit cell C Λ , we find the well known Sellmeier fit [33] of the frequency dependence of the index of refraction n (ω) is nicely reproduced from the eigenvalues of the dielectric tensor ε (T ) ab (q = 0, ω) for a variety of ionic compounds, see exemplarily our results for CsI and RbCl in Fig.10 and CaF 2 and BaF 2 in Fig.11. The relative error of our calculations compared to a fit of experimental data for n (ω) with the Sellmeier formula for the mentioned crystals is less than 1%.
Let us emphasize our approach warrants notably fewer fit parameters compared to a Sellmeier fit. For example, to reproduce the experimentally observed chromatic dispersion of CsI over a wide frequency intervall, ranging from the ultraviolet to far-infrared, a satisfactory fit to the experimental data within our approach needs only two functions of the type (37) to model the induced electronic polarization of individual ions and a third fit function of the type (38) to model the ionic polarization effect. So our fit relies only on 6 parameters modelling the microscopic polarizabilities of atoms (ions, molecules, ion pairs) compared, for example, to the 17 parameters required by the Sellmeier fit [34], describing chromatic dispersion of the refractive index of CsI.
Furthermore, for (ultraviolet) light propagating along the diagonal of the x-y plane, i.e. q = |q| √ 2 e (x) + e (y) , the dielectric tensor ε (T ) ab (q, ω) reveals two transversal modes capable to propagate with slightly different speeds inside the crystals mentioned above, thus causing an intrinsic birefringence ∆n (ω) induced by spatial dispersion. Our calculation of ∆n (ω) for the afore mentioned ionic crystals and a comparison with experimental data can be found in Fig.10 and Fig.11. The applied fit parameters entering the calculations of n (ω) and ∆n (ω) are listed in Table II. Having thus determined the model polarizabilities (37) and (38) for each (different) atom species and each ion pair, the dependence of the dielectric function on wave vector q is in our approach already fixed by the crystalline structure of the material under consideration, i.e. the rotary power γ abc (ω) and the dispersion induced anisotropy α abcd (ω) are already implicitely encoded in the q−dependence of the transversal dielectric tensor ε (T ) ab (q, ω). To what large extend our calculations agree with published experimental data over a wide range of optical frequencies for a series of quite different crystalline materials we summarize in Table I, and in particular in Fig. 10 and Fig.11.
While the refractive indices are deduced from the square root of the (real) eigenvalues of the transversal dielectric function for q = 0, the rotary power is determined from the imaginary part of the off-diagonal elements of ε (T ) ab (q, ω) for wave propagation along the crystals' (optical) z-axis, i.e. q = |q| e (z) . The examples of crystal structures listed in Table I cover the cubic crystal system as well as all uniaxial crystal systems, where the number M of ions comprising the unit cell C Λ varies between M = 4 (for e.g. hexagonal BeO) and M = 66 (for e.g. cubic Bi 12 SiO 20 and Bi 12 T iO 20 ). It should be pointed out, that in contrast to the results shown in Fig.10 and 11, our calculations for the refractive index as well as for the rotary power, both presented in Table I, solely rest on published data of (anisotropic) electronic polarizabilities being reported in the particular cited references.
As a side remark let us point out, that the described principal effects of non locality, the optical activity γ abc (ω) and/or the dispersion induced anisotropy α abcd (ω), at first sight being small effects compared to n 2 a − 1 with n 2 a representing the eigenvalues of the tensor ε (T ) ab (q, ω) in ordinary crystalline materials, could well be comparable to n 2 a − 1 in artificial periodic structures choosing appropriately taylored super lattices, see [38]. I. Calculation of the index of refraction n (calc) and rotary power ρ (calc) for wave propagation along the (optical) zaxis for various crystals and comparison with experimental results. The applied data for crystal structures and electronic polarizabilities entering our calculations, as well as the experimental data for the refractive index n (exp) and rotary power ρ (exp) are taken from the publications cited in the column "references". In case of anisotropic polarizabilities, see e.g. the uniaxial crystals TiO2 and CaCO3, α ′ and α ⊥ ′ denote the polarizabilities parallel and perpendicular to the optical z-axis, respectively.   Table  II. The relative error compared to a fit of experimental data for n (ω) with the multi-parameter Sellmeier formula [34] (blue line) is less than 1%. The spatial dispersion induced birefringence ∆n as calculated from the q-dependence of the macroscopic dielectric function (97) is displayed (red) in (c) for CsI and (d) for RbCl. To compare with experimental data [35] (blue dots) the orientation of the q−vector was chosen along the diagonal of the x-y plane. The estimated error in reading from the plots of the experimental data in [35] is about ±0.2 · 10 −6 .

Static Limit of the Dielectric Function for Monoatomic Bravais Lattices
Assuming for simplicity M = 1, then there is no loss of generality setting η (1) = 0. Identifying then, see (57), one readily infers from (94) the explicit representation   Table II. The relative error compared to a fit of experimental data for n (ω) with the multi-parameter Sellmeier formula [36] (solid lines) is less than 1%. The spatial dispersion induced birefringence ∆n as calculated from the q-dependence of the macroscopic dielectric function (97) is displayed (red) in (b) for CaF2 and (c) for BaF2. To compare with experimental data [37] (blue dots) the orientation of the q−vector was chosen along the diagonal of the x-y plane.
Elementary matrix algebra leads then to the result In the supplemental material [6] it is shown that The lattice sum (141) is conveniently evaluated along the lines indicated by Ewald [7,10], splitting the sum into two absolutely converging sums, one converging rapidly in the Fourier domain and the other converging rapidly in the spatial domain. For details and a discussion of our (modified) splitting method, see supplemental material [6]. For simple cubic lattices the expression (141) can also be evaluated employing Jacobi theta functions [67,68].
In the static limit, first |q| → 0 and then ω → 0, we write A comprehensive analysis of the general properties of electromagnetic response functions, quite apart from a particular model description of a material and only based on general principles such as causality and thermodynamic stability, has been given by Kirzhnitz [21], who derived for the isotropic case [ε Λ ] a,a ′ = εI a,a ′ as a lower bound of allowed static values of the dielectric function: Introducing the (dimensionless) Lorentz factors solely dependent on the lattice geometry, and identifying 1 |CΛ| = |ΛP | |ΩP | = ν P with the density of polarizable atoms in the probe volume, there follows from (140) for mono-atomic Bravais lattices (M = 1) the following exact formula for the static dielectric tensor: While the polarizability α refers to individual atomic (ionic, molecular) properties, the dielectric tensor ε Λ also depends via the 3 × 3 matrix of Lorentz factors L on how the atoms are assembled to build the crystal Λ. Note that the particle density ν P of a dielectric crystalline probe volume is bounded from above by the a critical value ν (c) P , with L max denoting the largest positive eigenvalue of the matrix L. This condition indeed implies the matrix ε Λ − I being positive definit, thus generalizing the stability criterion (143) to the anisotropic case.
It should be noted that the trace of the matrix L as defined in (144) is normalized to unity: For a proof see supplemental material [6]. Numerical values for the Lorentz factors are readily calculated for all 14 monoatomic Bravais lattices Λ along the lines indicated in the supplemental material [6]. In the special case of a lattice with tetragonal or hexagonal symmetry the matrix L becomes diagonal, L a,a ′ = δ a,a ′ L a . In this case the trace identity (147) leads to L x = L y = 1 2 (1 − L z ), with L z being a universal function, solely dependent on the ratio az ax of the lattice constants, a z parallel and a x perpendicular to the crystalline z−axis. In Fig.  12 that function L z is plotted vs. the ratio of lattice constants az ax for simple tetragonal (st), body-centered tetragonal (bct) and also hexagonal (hex) lattice symmetry. Our results for the Lorentz-factors obtained with Eq. (144) coincide with previous works, for example [69][70][71].

Clausius-Mossotti Relation
For cubic symmetry there holds L a,a ′ = 1 3 δ a,a ′ and the static dielectric tensor (145) reduces to the well known Clausius-Mossotti relation for isotropic systems, see supplemental material [6]: Incidentally, the relation (148) applies for a wide class of (isotropic, non polar) materials, including dielectric liquids and gases. On a final note: our derivation of (148) completely avoids the usual trick of introducing the Lorentz sphere, where the medium outside of this sphere is considered as a continuum. For an in-depth explanation of that trick see for example [24].

V. CONCLUSIONS
The field-integral equation approach presented in this article differs from traditional presentations of crystal optics, for example [12,13,15]. Based on the Helmholtz-Hodge theorem the source term in the microscopic Maxwell equations representing the current density has been decomposed into longitudinal and transversal parts, thus establishing the formulation of equivalent (inhomogenous) field-integral equations with a kernel modelling the induced microscopic polarizationP (r, ω) inside a dielectric crystal as a convolution integral of the dielectric susceptibility χ aa ′ (r, r ′ , ω) (35) and the Fourier amplitudeẼ (r, ω) of the microscopic local electric field, see (32). Exploiting the lattice periodicity of the dielectric susceptibility tensor χ aa ′ (r + R, r ′ + R, ω) = χ aa ′ (r, r ′ , ω) it is then natural to expand the field amplitudeẼ (r, ω) into a complete and orthonormal basis of eigenfunctions of the operator T R generating translations by a lattice vector R ∈ Λ, see (41). But instead of expanding the solution to the field-integral equations in the well known basis of plane waves e i(q+G)r , constructed from eigenfunctions of the momentum operator, thus requiring to handle for each wave vector q in the Brillouin zone C Λ −1 of the lattice Λ then (infinite dimensional) matrices labelled by reciprocal lattice vectors G, G ′ ∈ Λ −1 , we use non-standard Bloch functions representing a complete and orthonormal system of eigenfunctions of the translation operator T R constructed from eigenfunctions of the position operator, see (50). Our choice of basis functions indeed enables to sidesteps the inversion (and truncation) of large matrices with regard to the reciprocal lattice vectors G, G ′ ∈ Λ −1 , thus easing significantly the numerical effort.
Considering the homogenous field-integral equations the electromagnetic modes and the photonic band structure ω n (q) of a dielectric crystal have been identified solving a small sized 3M × 3M matrix eigenvalue problem, with M denoting the number of polarizable atoms (ions, molecules) in the elementary cell C Λ of the lattice. A radiation damping term ∝ ω 3 , that originates not from the damping terms in the atom-individual polarizabilities (2), but from the retarded Helmholtz propagator evaluating the lattice sum without the self field term, see (84), has been taken into account in all our calculations. Exemplarily we then presented the photonic band structure of diamond (M = 2). An overview of our results for the photonic bandstructure of (primitive) sc-, fcc-and bcc lattices we present in the supplemental material [6]. In all cases we find quantitative agreement with previously published work on photonic band structures calculated with other methods.
Considering the inhomogenous field-integral equations, choosing a (circular) frequency ω and a wavevector q obeying to ω = ω n (q), the microscopic local electric field amplitudeẼ (r, ω) in the presence of a slowly varying time harmonic external field with amplitudeẼ ext (r, ω) is found to display in general a radiative transversal part (T) and also a longitudinal part (L), the size of the longitudinal amplitudeẼ (L) (r, ω), see (79), compared to the size of the transversal amplitudeẼ (T ) (r, ω), see (80), being strongly depend on the density of polarizable atoms (ions, molecules) in the crystal, see Fig.5. In a sufficiently dense packed dielectric crystal we find the longitudinal amplitudeẼ (L) (r, ω) being substantially larger compared to the transversal amplitudeẼ (T ) (r, ω), but in a dilute (artificial) superlattice with polarizable subunits positioned at the lattice sites the longitudinal field amplitude becomes small compared to the transversal one, see Fig.5. While the microscopic local electric field amplitudeẼ (r, ω) displays inside a densely packed lattice rapid spatial variations with large amplitude traversing a distance set by the microscopic inter-particle distance a, it turns out that the transversal (radiative) partẼ (T ) (r, ω) essentially coincides with the slowly varying macroscopic electric field amplitudeẼ (r, ω), that fieldẼ (r, ω) being conceived as the low pass filtered microscopic local field amplitude, see (85). If the external fieldẼ ext (r, ω) was purely transversal and slowly varying, then the amplitude of the residue δẼ (T ) (r, ω) =Ẽ (T ) (r, ω) −Ẽ (T ) (r, ω) is tiny, see Fig.8.
Regarding the propagation of light signals inside a dielectric crystal then the question was clarified, which parts of the dielectric tensor ε Λ (q, ω) describe the renormalization of the speed of light and possibly birefringence, chromatic dispersion, rotary power and spatial-dispersion induced birefringence, and which parts govern electric-field screening. Accordingly we derived directly from the field-integral equations, see (109), a set of (coupled) differential equations for the longitudinal and transversal components determining the macroscopic electric fieldẼ (r, ω) directly, without any prior knowledge of the microscopic local electric fieldẼ (r, ω). If the external field was purely transversal then an effective wave equation for the transversal (radiative) macroscopic field emerged, thus identifying the pieces of ε Λ (q, ω) comprising the transversal dielectric tensor ε (T ) aa ′ (q, ω) , see (124). Conversely, if the external field was purely longitudinal, then an effective Poisson type equation for a scalar potentialφ (r, ω) turned up, so that the longitudinal macroscopic field is represented byẼ (L) (r, ω) = −∇φ (r, ω), thus identifying the pieces of the dielectric tensor comprising the longitudinal dielectric tensor ε (L) aa ′ (q, ω) being liable for electric-field screening (like in electrostatics), see (133). In the static limit, q → 0 and ω → 0, an exact expression for the dielectric tensor ε Λ is derived that applies for all 14 mono-atomic Bravais lattices, see (145). Our result for ε Λ in particular conforms with general (thermodynamic) stability criteria [21], and for cubic symmetry the well known Clausius-Mossotti relation is recovered.
The Taylor expansion of the transversal dielectric tensor, ε (T ) ab (q, ω) = ε (T ) ab (ω)+i c γ abc (ω) q c + c,d α abcd (ω) q c q d + .. around q = 0, provides then insight into various optical phenomena connected to retardation and non locality of the dielectric response, in full agreement with the phenomenological reasoning of Agranovich and Ginzburg [3]: the eigenvalues of the tensor ε (T ) ab (ω) describing chromatic dispersion of the index of refraction and birefringence, the first order term γ abc (ω) specifying rotary power (natural optical activity) in crystals lacking inversion symmetry, the second order term α abcd (ω) shaping the (weak) effects of a spatial-dispersion induced birefringence, nowadays a critical parameter for the design of lenses made from CaF 2 and BaF 2 for optical lithograpy systems in the ultraviolet.
Considering various dielectric crystalline materials comprising atoms (molecules, ions) with known polarizabilities from the literature, in all cases the calculated indices of refractions, the rotary power and the dispersion induced birefringence have been shown to coincide well with the experimental data displayed in table I, thus illustrating the utility of the theory. For ionic crystals, exemplarily for CsI and RbCl, a satisfactory agreement between theory and the measured chromatic dispersion of the index of refraction is manifest over a wide frequency interval, ranging from ultraviolet to far infra-red, accomplishing this with an appreciably smaller number of adjustable parameters directly linked to microscopic atomindividual polarizabilities compared to the well known Sellmeier fit. Further we computed for the cubic fluorites CaF 2 and BaF 2 the frequency dependence of the intrinsic birefingence, as represented by three independent components of the tensor α abcd (ω), and found good agreement with published data [19,37].
Even though all our calculations are based on a phenomenological model of microscopic polarizabilities, see (37) and (38), due to the conformity with the fundamental frequency dependence of atom-individual polarizabilities as predicted by quantum mechanics, see (2), our results for the dielectric tensor ε Λ (q, ω) may well claim general validity in the range of optical frequencies and below. While the frequency dependence of ε Λ (q, ω) is deeply anchored in the retarded response of polarizable atoms (ions, molecules) comprising a dielectric crystal, the dependence of ε Λ (q, ω) on the wave vector q of the propagating light describes the effects attributed to the non-locality of that response, for instance optical activity and intrinsic birefringence, all theses effects being primarily dependent on details of the crystalline symmetry. Finally it should be noted, that our theory of the macroscopic electric field doesn't make use of the notion of a displacement field D, and therefore we avoided to mention it.

Outlook
We are confident the presented field-integral equation approach can be extended to calculations of the microscopic local electric field near to the surface of a dielectric crystal and also to thin films. In addition we consider it possible to extend our approach to disordered systems, for instance crystals subject to substitutional site disorder, thus enabling a theory of the dielectric tensor for disordered systems within the frame of the coherent potential approximation (CPA).

ACKNOWLEDGMENTS
The corresponding author thanks Oleg Dolgov for insightfull remarks and for bringing the references [1, 2, 21] to his attention, all of this a good while ago. He also thanks Mario Liu for sharing with him over the years his thoughts on macroscopic electrodynamics. Last not least we thank Oliver Eibl, Klaus-Peter Federsel, József Fortágh, Reinhold Kleiner and Claus Zimmermann for useful discussions. [5] We conform to the convention, that the Fourier transformation of a function F (t) of time t and its Fourier inverseF (ω) as a function of (circular) frequency ω are defined by