Reflection and transmission in nonlocal susceptibility models with multiple resonances

TL;DR

This paper develops a comprehensive model for electromagnetic reflection and transmission in nonlocal dielectric materials with multiple resonances, incorporating boundary effects and exciton dispersion, and analyzes their spectral energy densities.

Contribution

It introduces exact formulas for reflection and transmission in multi-resonance nonlocal media with arbitrary boundary conditions, extending existing models to include exciton dispersion and multiple resonances.

Findings

Improved single-band approximation for exciton bands

Extended model to include linear $k$ terms in dispersion

Calculated spectral energy density for various multi-resonance models

Abstract

We consider a semi-infinite dielectric with multiple spatially dispersive resonances in the susceptibility. The effect of the boundary is described by an arbitrary reflection coefficient for polarization waves in the material at the surface, with specific values corresponding to various additional boundary conditions (ABCs) for Maxwell's equations. We derive exact expressions for the electromagnetic reflection and transmission coefficients and present the results for a variety of materials with multiple exciton bands. We find an improved single-band approximation for heavy/light exciton bands and extend our model to exciton dispersion relations with linear $k$ terms which occur in uniaxial crystals. Finally, we calculate the spectral energy density of thermal and zero-point radiation for a variety of multi-resonance models and ABCs.

Figures

Tables

Table 1. Table 1: List of ABC’s

	$U_x$	$U_y$	$U_z$
Agarwal et al.agarwal1971a ; agarwal1971b ; agarwal1972 ; birman1972 ; agarwal1973 ; maradudin1973 ; birman1974 ; mills1974 ; foley1975 ; bishop1976	0	0	0
Ting et al.ting1975	1	1	1
Fuchs-Kliewerting1975 ; kliewer1968 ; kliewer1971 ; kliewer1975 ; ruppin1981	1	1	-1
Rimbey-Mahanrimbey1974 ; rimbey1975 ; rimbey1976 ; rimbey1977 ; rimbey1978	-1	-1	1
Pekarpekar1958a ; pekar1958b ; pekar1958c ; pekar1959	-1	-1	-1

Table 2. Table 2: List of model parameters

	ZnOlagois1981			GaAssell1973
m	1	2	3	1	2
${\chi }_0$	5.2	5.2	5.2	11.6	11.6
$\hslash {\omega }_T$ (eV)	3.3758	3.3810	3.4198	1.514	1.514
$\hslash {\omega }_L$ (eV)	3.3776	3.3912	3.4317	1.515	1.515
$\hslash \gamma$ (meV)	0.7	0.7	0.7	0.05	0.05
$m_{\mathrm{ex}}$ ($m_{\mathrm{ex}}$)	0.87	0.87	0.87	0.183	0.805
$\hslash {\omega }_p$ (eV)	0.5334	0.6055	0.5983	0.138	0.138
$D$ (${10}^{11}{m}^2{s}^{-2}$)	6.82	6.84	6.91	14.55	3.31

Table 3. Table 3: List of χ lin ​ y subscript 𝜒 lin 𝑦 \chi_{{\rm lin}y} expressions for 𝒄 𝒄 {\bm{c}} orientations.

	${\chi }_{\mathrm{lin}y}$
$𝒄\parallel \widehat{𝒙}$	$\frac{2{\omega }_p^2\left({\omega }_T^2+D_{⟂}{q}^2+D_{\parallel }{K}^2-{\omega }^2-i\gamma \omega \right)}{{\left({\omega }_T^2+D_{⟂}{q}^2+D_{\parallel }{K}^2-{\omega }^2-i\gamma \omega \right)}^2-{\xi }^2{q}^2}$
$𝒄\parallel \widehat{𝒚}$	$\frac{2{\omega }_p^2}{\left({\omega }_T^2+D_{⟂}\left(K^2+q^2\right)\right)-{\omega }^2-i\gamma \omega }$
$𝒄\parallel \widehat{𝒛}$	$\frac{2{\omega }_p^2\left({\omega }_T^2+D_{⟂}{K}^2+D_{\parallel }{q}^2-{\omega }^2-i\gamma \omega \right)}{{\left({\omega }_T^2+D_{⟂}{K}^2+D_{\parallel }{q}^2-{\omega }^2-i\gamma \omega \right)}^2-{\xi }^2{K}^2}$

Table 4. Table 4: List of χ lin ​ x subscript 𝜒 lin 𝑥 \chi_{{\rm lin}x} and χ lin ​ z subscript 𝜒 lin 𝑧 \chi_{{\rm lin}z} expressions for 𝒄 𝒄 {\bm{c}} orientations.

	${\chi }_{\mathrm{lin}x}$	${\chi }_{\mathrm{lin}z}$
$𝒄\parallel \widehat{𝒙}$	$\frac{2{\omega }_p^2}{\left({\omega }_T^2+D_{⟂}{q}^2+D_{\parallel }{K}^2\right)-{\omega }^2-i\gamma \omega }$	$\frac{2{\omega }_p^2\left({\omega }_T^2+D_{⟂}{q}^2+D_{\parallel }{K}^2-{\omega }^2-i\gamma \omega \right)}{{\left({\omega }_T^2+D_{⟂}{q}^2+D_{\parallel }{K}^2-{\omega }^2-i\gamma \omega \right)}^2-{\xi }^2{q}^2}$
$𝒄\parallel \widehat{𝒚}$	$\frac{2{\omega }_p^2\left({\omega }_T^2+D_{⟂}\left(K^2+q^2\right)-{\omega }^2-i\gamma \omega \right)}{{\left({\omega }_T^2+D_{⟂}\left(K^2+q^2\right)-{\omega }^2-i\gamma \omega \right)}^2-{\xi }^2\left(K^2+q^2\right)}$	$\frac{2{\omega }_p^2\left({\omega }_T^2+D_{⟂}\left(K^2+q^2\right)-{\omega }^2-i\gamma \omega \right)}{{\left({\omega }_T^2+D_{⟂}\left(K^2+q^2\right)-{\omega }^2-i\gamma \omega \right)}^2-{\xi }^2\left(K^2+q^2\right)}$
$𝒄\parallel \widehat{𝒛}$	$\frac{2{\omega }_p^2\left({\omega }_T^2+D_{⟂}{K}^2+D_{\parallel }{q}^2-{\omega }^2-i\gamma \omega \right)}{{\left({\omega }_T^2+D_{⟂}{K}^2+D_{\parallel }{q}^2-{\omega }^2-i\gamma \omega \right)}^2-{\xi }^2{K}^2}$	$\frac{2{\omega }_p^2}{\left({\omega }_T^2+D_{⟂}{K}^2+D_{\parallel }{q}^2\right)-{\omega }^2-i\gamma \omega }$

Table 5. Table 5: List of model parameters

	CdShonerlage1984
${\chi }_0$	6.5
$\hslash {\omega }_T$ (eV)	2.5674
$\hslash {\omega }_L$ (eV)	2.5688
$\hslash \gamma$ (meV)	0.075
$m_{\mathrm{ex}}$ ($m_{\mathrm{ex}}$)	1.3 ($⟂𝒄$) , 1.02 ($\parallel 𝒄$)
$\zeta$ (eV m)	5.6$\times {10}^{-12}$
$\hslash {\omega }_p$ (eV)	0.164
$D$ (${10}^{11}{m}^2{s}^{-2}$)	3.47 ($⟂𝒄$), 4.43 ($\parallel 𝒄$)
$\xi$ (${\mathrm{ms}}^{-2}$)	6.637$\times {10}^{19}$

Equations

$$\chi (\mathbf{k},\omega )={\chi }_0+\sum _{m=1}^M\frac{{\omega }_{pm}^2}{{\omega }_{Tm}^2(k)-{\omega }^2-i{\gamma }_m\omega },$$

$$\hslash {\omega }_T(k)=\hslash {\omega }_T+\frac{{\hslash }^2{k}^2}{2m_{\mathrm{ex}}},{\omega }_T^2(k)\approx {\omega }_T^2+D{k}^2$$

$$\nabla \times \nabla \times \mathbf{E}(\mathbf{r},\omega )-\frac{{\omega }^2}{c^2}\mathbf{E}(\mathbf{r},\omega )=\frac{{\omega }^2}{c^2}\mathbf{P}(\mathbf{r},\omega ),$$

$$P_i(\mathbf{r},\omega )=\int d^3{\mathbf{r}}^{\prime }\chi (\mathbf{r}-{\mathbf{r}}^{\prime },\omega )E_i({\mathbf{r}}^{\prime },\omega ).$$

$$P_i(\mathbf{r},\omega )=\frac{1}{(2\pi {)}^3}\int d^3\mathbf{k}{P}_i(\mathbf{k},\omega )e^{i\mathbf{k}\cdot \mathbf{r}}$$

$$P_i(\mathbf{k},\omega )=\sum _j\chi (\mathbf{k},\omega )E_i(\mathbf{k},\omega ).$$

$$(\omega /c{)}^2\left[1+\chi (\mathbf{k},\omega )\right]=k^2,$$

$$1+\chi (\mathbf{k},\omega )=0,$$

$$\chi (K,0,q)=$$

$${\chi }_m(K,0,q)=\frac{{\omega }_{pm}^2/D_m}{q^2-{\Gamma }_m^2(K)}$$

$${\Gamma }_m^2(K)=\frac{{\omega }^2-{\omega }_{Tm}^2+i{\gamma }_m\omega }{D_m}-K^2.$$

$$P_i(K,0,z)={\int }_0^{\infty }d{z}^{\prime }{\tilde{\chi }}_i^{\prime }(K,0,z,z^{\prime })E_i(K,0,z^{\prime }).$$

$${\tilde{\chi }}_{mi}^{\prime }(z,z^{\prime })=\{\begin{array}{rlr}& {\tilde{\chi }}_m(z-z^{\prime })+U_{mi}{\tilde{\chi }}_m(z+z^{\prime })& \text{if }z,z^{\prime }>0,\\ & 0& \text{otherwise}\end{array}$$

$${\tilde{\chi }}_i^{\prime }(z,z^{\prime })=\{\begin{array}{rlr}& {\chi }_0\delta (z-z^{\prime })+\sum _m^M{\tilde{\chi }}_{mi}^{\prime }(z,z^{\prime })& \text{if }z,z^{\prime }>0,\\ & 0& \text{otherwise}.\end{array}$$

$$P_i(z)=$$

$$\displaystyle\bigg{[}e^{iq(z-z^{\prime})}+U_{mi}e^{iq(z+z^{\prime})}\bigg{]}\chi_{m}(q)E_{i}(z^{\prime}),\quad z>0.$$

$$E_i(z)=\sum _{n=1}^N{E}_i^{(n)}{e}^{i{q}_{n}z},$$

$$P_i(z)=$$

$$\displaystyle\bigg{[}\frac{1}{q_{n}-q}+\frac{U_{mi}}{q_{n}+q}\bigg{]}\chi_{m}(q)E_{i}^{(n)},\quad z>0.$$

$$P_i(z)$$

$$-\frac{{\omega }_{p1}^2}{2D_1{\Gamma }_1}\sum _{n=1}^N\left(\frac{1}{q_n-{\Gamma }_1}+\frac{U_{1i}}{q_n+{\Gamma }_1}\right)E_i^{(n)}{e}^{i{\Gamma }_{1}z}$$

$$-\frac{{\omega }_{p2}^2}{2D_2{\Gamma }_2}\sum _{n=1}^N\left(\frac{1}{q_n-{\Gamma }_2}+\frac{U_{2i}}{q_n+{\Gamma }_2}\right)E_i^{(n)}{e}^{i{\Gamma }_{2}z}$$

$$-\dots$$

$$-\frac{{\omega }_{pM}^2}{2D_M{\Gamma }_M}\sum _{n=1}^N\left(\frac{1}{q_n-{\Gamma }_M}+\frac{U_{Mi}}{q_n+{\Gamma }_M}\right)E_i^{(n)}{e}^{i{\Gamma }_{M}z}.$$

$$\sum _{n=1}^N{\varphi }_{mi}^{(n)}{E}_i^{(n)}=0,m=1,\dots ,M$$

$${\varphi }_{mi}^{(n)}=\left(\frac{1}{q_n-{\Gamma }_m}+\frac{U_{mi}}{q_n+{\Gamma }_m}\right).$$

$$Z_p=\frac{E_x(0^{+})}{H_y(0^{+})}.$$

$$B_y^{(n)}(z)$$

$$=\left[\frac{q_n-K{\eta }^{(n)}}{k_0}\right]E_x^{(n)}{e}^{i{q}_{n}z}$$

$$={\tau }^{(n)}{E}_x^{(n)}{e}^{i{q}_{n}z},$$

$$E_z^{(n)}={\eta }^{(n)}{E}_x^{(n)},$$

$${\tau }^{(n)}=\{\begin{array}{rlr}& \frac{q_n^2+K^2}{q_n{k}_0},& \text{transverse waves},\\ & 0,& \text{longitudinal waves}.\end{array}$$

$$Z_p=$$

$$=$$

$$\sum _{n=1}^N\left[{\varphi }_{mx}^{(n)}\right]E_x^{(n)}=0,\sum _{n=1}^N\left[{\eta }^{(n)}{\varphi }_{mz}^{(n)}\right]E_x^{(n)}=0.$$

Full text

I Introduction

The susceptibility describing the material response to an applied electromagnetic field contains both temporal and spatial dispersionLLcm ; rukhadze1961 . As a result of the latter, the induced polarization at any point in the medium depends on the applied field in the region surrounding that point. This behavior can be due to a range of excitations within the material, such as phonons or excitons, each with their own associated resonance in the susceptibility. While this nonlocal response is often negligible in comparison to the frequency dependence there are cases where it can play a significant role, such as metallic nanostructuresraz11 ; wie12 ; pen12 ; tos15 ; sch16 , radiative heat transfersin15a ; sin15b , spontaneous emissionpur46 ; dre68 ; bar98 ; Novotny ; hor14 , spectral energy densitychurchill2016b and Casimir self-forceshor14 . Spatial dispersion is also important in semiconductors, where the complex electronic band structure can lead to many excitationshonerlage1984 , each with their own nonlocal response. In this paper, we continue our previous work on nonlocal responsechurchill2016b by extending results on reflection and transmission at planar boundaries to the case where the medium has multiple spatially-dispersive resonances.

The susceptibility is typically expressed as a sum of resonances. Nonlocal response is included as a $k$ dependence in the model parameters, but is usually limited to a $k$-dependent resonant frequency. Hopfield and Thomashopfieldthomas proposed the following model, based on the properties of semiconductors, but it can also be derived from a simple classical modelchu16 :

[TABLE]

where $\omega_{T}(k)$ is the resonant frequency, $\gamma$ quantifies the absorption and $\omega_{p}$ is the oscillator strength. The term $\chi_{0}$ collects contributions from other resonances and acts as a background susceptibility. For the sake of simplicity, the parabolic dispersion

[TABLE]

was usedhopfieldthomas to describe the exciton bands in the medium, where $m_{\rm ex}$ is the exciton mass and $D=\hbar\omega_{T}/m_{\rm ex}$

The difficulties involved in the calculation of electromagnetic reflection and transmission coefficients for nonlocal media are well known. While there is only one transmitted wave for a local medium, the nonlocal medium has several transmitted waves due to the $k$ dependence in (1)rukhadze1961 . The Maxwell boundary conditions are no longer sufficient to solve for the unknown amplitudes of the multiple transmitted waves. Historically, this need for extra information was resolved with the introduction of Additional Boundary Conditions (ABCs) on the polarization ${\bm{P}}$ associated with the resonances in (1). Various authorsagarwal1971a ; agarwal1971b ; agarwal1972 ; birman1972 ; agarwal1973 ; maradudin1973 ; birman1974 ; mills1974 ; foley1975 ; bishop1976 ; ting1975 ; kliewer1968 ; kliewer1971 ; kliewer1975 ; ruppin1981 ; rimbey1974 ; rimbey1975 ; rimbey1976 ; rimbey1977 ; rimbey1978 ; pekar1958a ; pekar1958b ; pekar1958c ; pekar1959 have proposed different ABCs under certain assumptions that suit different types of material. The Pekar ABC, where ${\bm{P}}$ vanishes at the boundary, is the simplest and most commonly used.

The majority of work on the subject has focused on susceptibilities with a single, isolated resonance. In this case, Halevi and FuchsHalevi have derived reflection coefficients for a generalized ABC model, containing all the previously suggested ABCs. In a previous paperchurchill2016b we have adapted this model to the tensor case with different transverse and longitudinal susceptibilities.

In general, the susceptibility of real materials is far more complex than simple isolated resonanceshonerlage1984 . In exciton bands, for example, there can be multiple closely spaced bands, degenerate bandshonerlage1984 ; kane1975 ; altarelli1977 and more complex $k$ dependencekoteles1979 ; koteles1980 ; Fiorini1980 . While some authorslagois1977 ; lagois1981 ; venghaus1978 ; sell1973 ; sermage1981 ; mahan1964 ; halevi1985 ; halevi1987 have considered multi-resonance systems, they are typically limited to a maximum of two resonances and a specific ABC.

The first aim of this paper is to extend the Halevi and FuchsHalevi generalized ABC model to a multi-resonance susceptibility and derive expressions for the reflection and transmission coefficients. This derivation is first applied to a system with simple parabolic exciton bands and then to bands that are degenerate at $k=0$, where we find improved parameters for the single band approximation.

The second aim is to modify the derivation further to include alternate wave vector dependences, specifically the case where the dispersion in (2) contains a $\pm k$ term. This behavior is known as linear splitting and is typically found in uniaxial crystals such as Wurtzite. This case has been previously calculatedmahan1964 ; halevi1985 ; halevi1987 , but only for a specific ABC and orientation of the crystal axis. We will show that linear $k$ splitting can easily be incorporated into the multi-resonance model and that the orientation of the crystal axis has significant effects on the result.

Finally, we will use the derived electromagnetic reflection coefficients to calculate the spectral energy density of thermal and zero-point radiation outside the spatially dispersive medium. The results for various multi-resonance systems are calculated and compared to those in our previous paper on the isolated resonancechurchill2016b .

The assumptions made in the following derivation are discussed in greater detail in our previous paperchurchill2016b . In summary, we consider a smooth boundary that does not contain any features such as slits or other nontrivial structures and is sufficiently far from any other boundaries such that multiple reflections can be ignored. In addition we do not include any quantum mechanical effects not directly encoded in the macroscopic susceptibility

The paper is organized as follows. In Sec. II we present the spatially dispersive susceptibility model for a half-infinite dielectric with a multi-resonance permittivity and derive the field equations. In Sec. III and IV we derive the general expressions for the reflection and transmission coefficients for $p$ and $s$ polarized light and present the results for parabolic exciton bands in Sec. V. In Sec. VI we extend our derivation to a uniaxial crystal by including a linear splitting term and in Sec. VII we calculate the zero-point and thermal spectral energy density.

II Theory

II.1 Infinite Medium

We first consider an infinite, homogeneous, spatially-dispersive dielectric with the susceptibility (1). The electric field ${\bm{E}}$ and polarization field ${\bm{P}}$ satisfy the wave equation:

[TABLE]

where the polarization field is given by

[TABLE]

In general, the spatially dispersive susceptibility $\chi$ is a tensorLLcm , but here we consider a scalar. Using the Fourier transformation

[TABLE]

we have

[TABLE]

The wave equation (3) has solutions for ${\bm{E}}$ when the frequency and wave vector satisfy the dispersion relationrukhadze1961

[TABLE]

for transverse waves with ${\bm{E}}\cdot{\bm{k}}=0$ or

[TABLE]

for longitudinal waves with ${\bm{E}}\times{\bm{k}}=0$. With the field dependence $\textrm{exp}(ik_{z}z)$ we restrict ourselves to Im$[k_{z}]>0$, leading to $M+1$ transverse and $M$ longitudinal waves for the susceptibility in (1) and the parabolic dispersion in (2).

II.2 Half-Infinite Medium

We now consider the half-infinite dielectric occupying the $z>0$ region as shown in Fig. 1. The vacuum contains the incident wave (${\bm{E}_{0}}$) and reflected wave (${\bm{E}_{r}}$) with wave-vectors $\bm{k}_{0}$ and $\bm{k}_{r}$ ($k_{0}=k_{r}=\omega/c$). Inside the dielectric there are $N=2M+1$ transmitted waves (${\bm{E}_{n}}$) with the corresponding wave vectors ${\bm{k}_{n}}$. The coordinate system has been chosen such that the $xz$-plane coincides with the plane of incidence, with $k_{nx}=K$, $k_{ny}=0$ and $k_{nz}=q_{n}$.

The bulk susceptibility (1) in this co-ordinate system can be rewritten as:

[TABLE]

where:

[TABLE]

and:

[TABLE]

With the presence of the boundary at $z=0$, the polarization field now depends on a position-dependent susceptibility $\chi^{\prime}_{i}$ ($i\in\{x,y,z\}$)Halevi . After a Fourier transformation in the $xy$ plane:

[TABLE]

We subsequently omit $K$ dependence for notational simplicity. We assume that each resonance in the half-infinite medium can be expressed in terms of the bulk susceptibility in the same manner as Halevi and FuchsHalevi :

[TABLE]

and the overall susceptibility is given by:

[TABLE]

The first term in (13) is the position-independent nonlocal bulk response. The second describes a polarization wave propagating from $z^{\prime}$ to the surface before reflecting and continuing to $z$. The reflection amplitude coefficient $U_{i}$ is (in general) complex and frequency dependent, with $|U_{i}|=1$ implying elastic reflection. Halevi & FuchsHalevi demonstrated that specific values of $U_{i}$ correspond to certain ABC’s, shown in Table 1. Each ABC was developed for a particular type of medium or excitation. For example, Pekarpekar1958a ; pekar1958b ; pekar1958c ; pekar1959 and Rimbey-Mahanrimbey1974 ; rimbey1975 ; rimbey1976 ; rimbey1977 ; rimbey1978 were developed for Frenkel (tight-binding) excitons, Ting et al.ting1975 for Wannier-Mott (weak-binding) excitons, Fuchs-Kleiwerting1975 ; kliewer1968 ; kliewer1971 ; kliewer1975 ; ruppin1981 for metals and Agarwal et al.agarwal1971a ; agarwal1971b ; agarwal1972 ; birman1972 ; agarwal1973 ; maradudin1973 ; birman1974 ; mills1974 ; foley1975 ; bishop1976 for the general case where surface effects can be ignored.

Substituting (14) into (12) gives:

[TABLE]

At this point we introduce an ansatz for the ${\bm{E}}$ field inside the mediumHalevi — a linear combination of $N=2M+1$ plane waves from (7) and (8):

[TABLE]

where $n=1$ to $M+1$ are transverse waves and $n=M+2$ to $N$ are longitudinal waves.

After substitution of the ansatz (16) into (15) and evaluating the $z^{\prime}$ integral, we find

[TABLE]

The $q$ integral is evaluated by performing a contour integration in the upper half-plane. This encloses the poles at $q=q_{n}$ and $\Gamma_{m}$, giving

[TABLE]

For the wave equation (3) to be valid for all values of $z$, we require each of the right-hand side sums proportional to $\textrm{exp}(i\Gamma_{m}z)$ in (18) to equal zero:

[TABLE]

where

[TABLE]

This leads to a set of $M$ equations of the form (19) for each of the $E_{i}$ components.

III $p$-polarization

The field can be decomposed to components with $\bm{E}$ perpendicular to ($s$-polarized) or in the plane of incidence ($p$-polarized). For $p$-polarized light $E_{y}=0$, $E_{x}\neq 0$ and $E_{z}\neq 0$. Both the transverse and longitudinal waves appear in the medium.

III.1 Surface Impedance

The reflection coefficient is calculated using the surface impedance, which for $p$-polarized light is given by:

[TABLE]

(Here ${\bm{H}}=\mu_{0}{\bm{B}}$.) The magnetic field $B_{y}$ can be expressed in terms of the electric field using $k_{0}{\bm{B}}={\bm{k}}\times{\bm{E}}$ and (16):

[TABLE]

where we have substituted $E_{x}$ for $E_{z}$ using

[TABLE]

where $\eta^{(n)}=-K/q_{n}$ for transverse waves and $\eta^{(n)}=q_{n}/K$ for longitudinal waves. This leads to:

[TABLE]

The surface impedance can now be expressed solely in terms of $E_{x}$ field amplitude ratios:

[TABLE]

III.2 Field Amplitude Ratios

To proceed any further, we require the $E_{x}$ field amplitude ratios in (25). By using (23), we can rewrite the $E_{z}$ equations in (19) in terms of $E_{x}$:

[TABLE]

We now have $2M$ equations relating the $2M+1$ waves inside the medium and have sufficient information to solve for the reflection coefficient. After dividing by $E^{(1)}_{x}$ and rearranging we can express (26) in matrix form. As an example, we present the result for a two-resonance system:

[TABLE]

where $n=1,2,3$ are transverse waves and $n=4,5$ are longitudinal waves. By inverting the $2M\times 2M$ matrix, we can find the field amplitude ratios.

III.3 Reflection and Transmission Coefficients

The $p$-polarization reflection coefficient can be expressed in terms of surface impedancekliewer1968

[TABLE]

where $Z_{p}$ is given by (25) and $Z_{p}^{(0)}=\sqrt{k_{0}^{2}-K^{2}}/\mu_{0}k_{0}$ is the vacuum surface impedance.

We can find the transmission coefficients for the $N$ transmitted waves by imposing the continuity of the tangential ${\bm{E}}$ field across the boundary. Our choice of coordinate system means we simply equate the $E_{x}$ components on each side:

[TABLE]

This can be expressed in terms of the previously calculated field amplitude ratios using (28) and $\cos{\theta_{i}}=\sqrt{k_{0}^{2}-K^{2}}/k_{0}$:

[TABLE]

By rewriting:

[TABLE]

we can derive the transmission coefficient:

[TABLE]

for transverse waves:

[TABLE]

and longitudinal waves:

[TABLE]

We now have a single method for $r_{p}$ and $t_{p}$ in the presence of multiple resonances that can cover a wide frequency range. This derivation could be extended further to a tensor susceptibility using the method described in our previous paperchurchill2016b . In this case the polarization reflection coefficents are a tensor $U_{ij}$ with certain restrictions on the components.

IV $s$-polarization

We now consider the simpler case of $s$-polarized light, where $E_{y}\neq 0$, $E_{x}=0$ and $E_{z}=0$. As ${\bm{k}_{n}}$ all lie in the $xz$-plane, this leads to the absence of longitudinal waves in the medium, leaving the $M+1$ transverse waves.

IV.1 Surface Impedance

The surface impedance for $s$-polarized light is given by:

[TABLE]

As in (25) we express $H_{x}$ in terms of $E_{y}$ and field amplitude ratios:

[TABLE]

As we only have the $M+1$ transverse waves in the $s$-polarization, the set of $M$ equations from (19) can be rewritten as

[TABLE]

which is sufficient to solve for the amplitude ratios.

IV.2 Reflection and Transmission Coefficients

Using the vacuum surface impedance $Z_{s}^{(0)}=k_{0}/\mu_{0}\sqrt{k_{0}^{2}-K^{2}}$ and (36), we can construct the $s$-polarization reflection coefficient:

[TABLE]

As in the previous section, we impose the continuity of the tangential ${\bm{E}}$ field across the boundary. As we only have $E_{y}$ components, this leads to:

[TABLE]

Using (38) this can be expressed in terms of the field amplitude ratios previously found from (37):

[TABLE]

which leads to:

[TABLE]

V Reflection Coefficient Results

We now use the derivations in the previous sections for a number of materials with a variety of exciton band structures. In Fig. 2 we show some example exciton bands and the corresponding dispersion relations for transverse ${\bm{E}}$ waves in the absence of damping.

In this section we first consider the simple case of ZnOlagois1981 with three non-interacting exciton bands and GaAssell1973 with two bands that are degenerate at $k=0$. The model parameters are given in Table 2. The values of $\omega_{p}$ are calculated from the measured values of $\omega_{L}$, which are the solutions of the dispersion relation for transverse ${\bm{E}}$ waves at $k=0$ in the absence of damping. Similarly $D$ is found from the measured exciton mass $m_{\rm ex}$, which is given in units of the rest electron mass $m_{\rm e0}$.

V.1 Simple Resonances

We first consider a three-resonance model for ZnOlagois1981 , involving the A, B and C excitons which we label $m=1$, $2$ and $3$ respectively. The exciton bands are of the form in (2) and do not interact.

Figure 3 shows $r_{p}(\omega)$ and $r_{s}(\omega)$ for a fixed incident angle. The peak locations are determined by the $\omega_{T}$ and $\omega_{L}$ values, indicated by solid and dashed vertical lines respectively. The $r_{p}$ behavior is mostly determined by $U_{x}$, with $U_{x}=1$ giving the largest maxima and smallest minima. In contrast $U_{z}$ only affects the results at the reflection minima. The frequency region around $\omega_{T3}$ can be accurately described using the single-resonance model. This is because the overlap with the other resonances in the susceptibility is very small. The same is not true of the $m=1$ and $2$ resonances. Here the proximity of $\omega_{T1}$ and $\omega_{T2}$ lead to the one-resonance models failing, particularly in intermediate frequency region. The effect of $U_{i}$ remains the same as the isolated resonance. However, if there is significant overlap in the resonant peaks (e.g. if the $\omega_{T}$ values are separated by less that the full-width half-maximum), then $U_{x}=1$ gives not only the the largest peaks but also the largest value in the intermediate frequency region.

V.2 Heavy & Light Excitons

We now move on to consider an exciton band structure with degeneracy at $k=0$. Kanekane1975 showed that interactions in a medium could lead to the splitting of degenerate exciton bands. In the case of isotropic valence bands, this can lead to a “heavy” and “light” exciton band with parabolic dispersion relations

[TABLE]

substituted into the susceptibility.

We consider a two-resonance model for GaAssell1973 containing only the heavy and light exciton bands, using the parameters in Table 2. We have slightly simplified the model by assuming isotropic valence bands and using the $\left<100\right>$ exciton masses for all directions. The results for $r_{p}$ are shown in Fig. 4 and display the same basic features as a single resonance model. The same behavior is seen in the $r_{s}$ result, which we do not plot.

Previous worksell1973 has suggested that the two-resonance system of heavy and light exciton bands can be approximated by a single band with an effective $D_{\rm eff}$ term in the susceptibility, given by

[TABLE]

and multiplying $\omega_{p}^{2}$ by a factor of 2. However, simply taking the average of the nonlocal parameter $D$ as in (44) does not lead to the best approximation. Figure 5 compares the heavy/light exciton model to the effective single band result and shows that it underestimates the peak of the heavy/light system. This is true even when $D_{\rm h}$ and $D_{\rm l}$ are close. Instead, if we re-express the nonlocal term as

[TABLE]

we find that a better fit is given by taking the average value of the coefficient $\sigma$, which leads to:

[TABLE]

This new value provides an excellent fit when $D_{\rm l}$ and $D_{\rm h}$ have similar values. For larger differences between $D_{\rm l}$ and $D_{\rm h}$ (such as this model where $D_{\rm l}\approx 4D_{\rm h}$), both single-resonance approximations begin to fail as the the two resonance model has a larger peak just above $\omega_{T}$. This difference is greatest for $U_{x}=-1$ and smallest for $U_{x}=1$. Despite this, we find that (46) gives a better fit to $r_{p}$ and $r_{s}$ than (44) for all values of $D_{\rm l/h}$ and $U_{i}$.

VI Linear $k$ terms

So far we have only considered materials with isotropic, parabolic energy bands of the form (2). However the symmetry of the crystal structure can lead to the introduction of additional $k$ terms. For example, in 1964 Mahan and Hopfieldmahan1964 used a linear $k$ term to explain a shoulder in the reflection spectra of CdS, a uniaxial medium with wurtzite crystal symmetry. This behavior was only observed when ${\bm{E}}$ was perpendicular to the crystal axis ${\bm{c}}$, with the exciton dispersion relation:

[TABLE]

where $m_{\rm ex}$ is the exciton mass and $k_{\parallel}$ and $k_{\perp}$ are wave vector components parallel and perpendicular to ${\bm{c}}$ respectively. In the susceptibility, (47) is approximated tomahan1964 :

[TABLE]

where $D_{\perp/\parallel}=\hbar\omega_{T}/m_{\rm ex\perp/\parallel}$ and $\xi=2\omega_{T}\zeta/\hbar$. This leads to two resonances:

[TABLE]

that can be collected to a single fraction $\chi_{\rm lin}^{\perp}$, giving

[TABLE]

for ${\bm{E}}\perp{\bm{c}}$. When ${\bm{E}}\parallel{\bm{c}}$, the exciton bands were found to be degenerate with $\zeta=0$, leading to the resonance:

[TABLE]

The bulk susceptibility of the uniaxial crystal is no longer a scalaer and takes a vector form when ${\bm{c}}$ is aligned with one of the co-ordinate axes, such as:

[TABLE]

for ${\bm{c}}\parallel\hat{{\bm{z}}}$. In Tables 3 and 4 we present the $\chi_{{\rm lin}i}$ components relevant for $s$ and $p$ polarized light when ${\bm{c}}$ is aligned with each of the co-ordinate axes as defined in Fig. 1.

Previous workmahan1964 ; halevi1985 ; halevi1987 has only considered such a model for ${\bm{c}}\parallel\hat{{\bm{y}}}$ (perpendicular to the incident plane) using the Pekar ABC. We will now modify the derivation of the previous section to include resonances of the form (52) in a multi-resonance system with arbitrary $U$ values when ${\bm{c}}$ is aligned with each of the co-ordinate axes.

VI.1 Field Amplitude Ratios

Due to the fact that the bulk components in (52) are no longer equal, the dispersion relations take the form

[TABLE]

for $s$-polarized light and

[TABLE]

for $p$-polarized light, where we have omitted $K$ dependence for notational simplicity. Unlike the previous section, the results of (53) are not also solutions of (54). We substitute the susceptibility (50-51) and the ansatz (16) into (15), using the $N$ values of $q_{n}$ that satisfy (53) and (54). If the field is aligned with ${\bm{c}}$, there is no linear splitting and the derivation in the previous section is sufficient to find the amplitude ratios. If linear splitting is present then $\chi_{{\rm lin}i}$ has two poles with Im$[q]>0$, which we label $\Gamma^{(+)}_{i}$ and $\Gamma^{(-)}_{i}$. Evaluating the contour integral in (17) gives

[TABLE]

where $F_{i}^{(\pm)}$ is a simple prefactor and we have used a single value of $U_{{\rm lin}i}$ for both resonances. The additional terms not proportional to ${\rm exp}(iq_{n}z)$ lead to the same set of equations as (19), but with different values of $\Gamma_{m}$:

[TABLE]

In summary, the inclusion of linear $k$ terms to the exciton dispersion relation does not significantly affect the derivation of field amplitude ratios presented in the previous section. The only changes required are to use the appropriate $\Gamma$ values in the field amplitude ratio matrix and the $q$ values that satisfy the dispersion relation.

VI.2 ${\bm{c}}\parallel\hat{{\bm{y}}}$ Orientation

The simplest case is to align the crystal axis ${\bm{c}}$ with $\hat{{\bm{y}}}$. This particular orientation has been looked at previously only for the Pekar ABC. As there is no linear splitting in $\chi_{{\rm lin}y}$ for this orientation, the derivation in the previous section is sufficient to calculate $r_{s}$ and $t_{s}^{(n)}$ in the $s$-polarization. The $\chi_{{\rm lin}y}$ term leads to two transverse wave from (53), if it is the only resonance.

In comparison, $\chi_{{\rm lin}x}$ and $\chi_{{\rm lin}z}$ are equal and contain linear splitting. Equation (54) for the $p$-polarization can be simplified to

[TABLE]

Solutions of the first bracket give two longitudinal waves and the second give three transverse waves, for a total of five if $\chi_{\rm lin}$ is the only resonance. The transverse waves are no longer the same as those in the $s$-polarization.

As $\chi_{{\rm lin}x}$ and $\chi_{{\rm lin}z}$ contain linear splitting, the integral in (17) leads to two equations of the form (56) for both $E_{x}$ and $E_{z}$. The $E_{z}$ equations can be converted to $E_{x}$ using (23) with $\eta^{(n)}=-K/q_{n}$ for transverse waves and $\eta^{(n)}=q_{n}/K$ for longitudinal waves as in the previous section. This gives a total of four $E_{x}$ equations, the same as the number of transmitted waves added by the resonance, which is sufficient to solve for the reflection coefficient.

VI.3 ${\bm{c}}\parallel\hat{{\bm{x}}}$ or ${\bm{c}}\parallel\hat{{\bm{z}}}$ Orientation

The two other orientations present additional challenges. As splitting is present in $\chi_{{\rm lin}y}$, the $s$-polarization now has three transverse waves, and two equations of the form (56) for $E_{y}$. There is sufficient information to solve for $r_{s}$ and $t_{s}^{(n)}$.

The fact that the expressions for $\chi_{{\rm lin}x}$ and $\chi_{{\rm lin}z}$ in Table 4 are different means (54) cannot be simplified and the waves are no longer purely transverse or longitudinal. The $\chi_{\rm lin}$ resonances give a total of four waves in the absence of other resonances. The contour integration in (17) leads to a total of three equations of the form (56) - one from the un-split $\chi_{{\rm lin}i}$ and two from the split $\chi_{{\rm lin}i}$. Again, we have sufficient information to solve for the field amplitude ratios.

As the waves are no longer purely transverse or longitudinal, care must be taken when converting $E_{z}$ to $E_{x}$. From the wave equation, we find:

[TABLE]

Similarly, the relation between $B_{y}$ and $E_{x}$ used in (25) is modified from (24) to:

[TABLE]

These expressions and the calculated field amplitude ratios are substituted into (25) and (36) to find the surface impedances and the subsequent reflection coefficients.

VI.4 Results

As an example, we present the results for CdShonerlage1984 using the model parameters in Table 5, where $\xi$ in (47) has been calculated from the measured $\zeta$ in (48). Figure 6 shows $r_{p}(\omega)$ and $r_{s}(\omega)$ at an incident angle of $60^{\circ}$ for ${\bm{c}}$ aligned with $\hat{{\bm{x}}}$, $\hat{{\bm{y}}}$ and $\hat{{\bm{z}}}$.

As expected, the ${\bm{c}}\parallel\hat{{\bm{y}}}$ result for $r_{p}$ is significantly different to the others due to linear splitting in both $\chi_{{\rm lin}x}$ and $\chi_{{\rm lin}z}$ components used for $p$ polarized light. This leads to an additional peak in $r_{p}$ below $\omega_{T}$. While the behavior of the new peak is still mostly determined by the value of $U_{x}$, the order in which the ABC’s appear is different to that of the previous section where $U_{x}=1$ is the largest and $U_{x}=-1$ the smallest. The $r_{s}$ result for ${\bm{c}}\parallel\hat{{\bm{y}}}$ is identical to that of the previous section as there is no splitting in $\chi_{{\rm lin}y}$ for this orientation.

The ${\bm{c}}\parallel\hat{{\bm{x}}}$ and ${\bm{c}}\parallel\hat{{\bm{z}}}$ cases have splitting in only one component of $\chi_{{\rm lin}i}$ for $p$ polarized light, resulting in an $r_{p}$ that is closer to the $\xi=0$ result but also displays new features just below $\omega_{T}$ at a smaller scale. The ${\bm{c}}\parallel\hat{{\bm{z}}}$ results for $r_{p}$ and $r_{s}$ are the closest to the $\xi=0$ case as the linear term in Table 4 only contains $K$. As $\theta_{i}$ (and $K$) is decreased, the peak at $\omega_{T}$ becomes smaller, returning to the $\xi=0$ result for normal incidence. As in the previous section, the new peak is unaffected by $U_{z}$ and is larger for $U_{x}=1$. In contrast, the ${\bm{c}}\parallel\hat{{\bm{x}}}$ case has $q$ in the linear splitting term. The differences in the $r_{p}$ result are similar in magnitude to the ${\bm{c}}\parallel\hat{{\bm{z}}}$ case, but the new peak at $\omega_{T}$ is more pronounced for $U_{x}=-1$ and is now affected by the value of $U_{z}$. The effect of the splitting in $q$ is even more pronounced in $r_{s}$, which displays features similar to $r_{p}$ in the ${\bm{c}}\parallel\hat{{\bm{y}}}$ case. This is because the linear splitting in $q$ is present in every $\chi_{{\rm lin}i}$ term used in their derivation.

In all cases, the difference between $m_{\perp}$ and $m_{\parallel}$ has very little effect compared to the choice of ABC. This agrees with our previous work on the tensor susceptibilitychurchill2016b .

VII Spectral Energy Density

We now focus on the electromagnetic zero-point and thermal radiation at a perpendicular distance $|z|$ from the boundary of the nonlocal medium. In our previous paper on the one-resonance systemchurchill2016b we found that the inclusion of spatial dispersion removed the unphysical $1/|z|^{3}$ divergence present in the spectral energy density of the local modelcandelas1982 ; henkel2000 . We now investigate how the behavior of the materials considered in the previous sections changes due to the presence of multiple resonances.

The average energy density of electromagnetic zero-point and thermal radiation in the vacuum outside a medium is given byjoulain2005

[TABLE]

where $u_{\rm tot}\left(z,\omega\right)$ is the spectral energy density. Assuming that the nonlocal medium is in thermal equilibrium with its surroundings and the system is rotationally invariant around the $z$ axis, this can be written in terms of the previously calculated reflection coefficients:

[TABLE]

The first integral in (61) is the contribution of propagating waves while the second comes from evanescent waves. The term $u_{0}$ is the spectral energy density in the absence of the material, given by

[TABLE]

where the mean energy of a harmonic oscillator in thermal equilibrium is

[TABLE]

The first term of (63) gives rise to the electromagnetic zero-point energy.

In the $K\to\infty$ limit for the local medium $r_{s}\to 0$ and $r_{p}\to\chi(\omega)/(2+\chi(\omega))$, leading to the divergent result for the second integraljoulain2005 :

[TABLE]

In our previous paper, we showed that the inclusion of spatial dispersion lead to peaks in $\textrm{Im}[r_{p}]$ near the point where the $\textrm{Re}[\Gamma^{2}]$ changed sign from positive to negative, followed by a $1/K^{4}$ decay in the large $K$ limit. In Fig. 7 we find the same behavior for ZnO, GaAs and CdS ($\hat{{\bm{z}}}\parallel{\bm{c}}$), with a peak for every $\Gamma_{m}$ value.

Figure 8 shows the $z$ dependence of $u_{\rm tot}$ for ZnO at $\hbar\omega=3.44{\rm eV}$. The behavior is the same as found in our previous paper - the $1/|z|^{3}$ divergence is removed and $u_{\rm tot}$ saturates to a finite value. The ABC behavior is similar, with Ting et al. giving the largest result, Fuchs-Kleiwer, Rimbey-Mahan and Agarwal et al. have similar intermediate values and Pekar is the smallest. The other materials display the same behavior and so are omitted here.

The main differences to the one-resonance model are found in the $\omega$ dependence of $u_{\rm tot}$ at a fixed distance from the boundary. We first consider the simple case of ZnO, with multiple, non-intersecting parabolic exciton bands. Figure 9 shows $u_{\rm tot}(\omega)$ at a distance of 8nm from the boundary. The results are strongly dependent on the choice of ABC, with Ting et al. giving the largest peaks, followed by Fuchs-Kleiwer, Agarwal et al., Rimbey-Mahan and finally Pekar. This behavior agrees with our previous paperchurchill2016b . Each resonance has two associated peaks in $u_{\rm tot}(\omega)$, for a total of six. The peaks in $u_{\rm tot}$ near $\omega_{Tm}$ are due to the $s$-polarization contribution to the integral and the peaks at $\omega_{Lm}$ are due to the $p$-polarization contribution. In the one-resonance system, the $s$-polarization peaks in $u_{\rm tot}$ were typically much smaller than their $p$-polarization counterparts. However, the $u_{\rm tot}$ peak at $\omega_{L1}$ has been suppressed due to the proximity of the peak at $\omega_{L2}$ and is now comparable in size to the $s$-polarization peaks. This can be seen in Fig. 7, where the peak in $\textrm{Im}[r_{p}]$ associated with $\Gamma_{1}$ and the $m=1$ resonance at the largest $K$ value is very small due to the presence of the nearby $\Gamma_{2}$ peak at a smaller $K$ value.

We next consider the heavy/light exciton model of GaAs. Figure 10 shows $u_{\rm tot}(\omega)$ at a distance of 8nm from the boundary. At first glance the results appear similar to those of the single resonance model, but the comparison in Fig. 11 reveals that both the $D^{*}$ and $D_{\rm eff}$ single-exciton models both underestimate the peak values. This is due to the behavior of $r_{p}$. While both approximations provide a good fit for propagating waves, $u_{\rm tot}$ depends more on the evanescent wave contribution at small distances. Both approximations have a single peak in $\textrm{Im}[r_{p}]$ at large $K$ in contrast to the two peaks of the heavy/light exciton model in Fig. 7.

We finally consider the uniaxial crystal of CdS in the ${\bm{c}}\parallel\hat{{\bm{z}}}$ case, as the other orientations lack the rotational invariance about the $z$ axis required for Eq. (61). The underlying behaviour behind Fig. 12 is more complex than the previous cases with $k^{2}$ dispersion. There are now two peaks in each of the $s$-and $p$-polarization contributions to $u_{\rm tot}$. A large peak in the $p$-polarization term is still found near $\omega_{L}$, but a new smaller peak is also present just above $\omega_{T}$. However, the peak in the $s$-polarization contribution to $u_{\rm tot}$ previously found at $\omega_{T}$ is now a minimum, with a larger peak below this frequency and a smaller peak above that coincides with the position of the lower peak in the $p$-polarization contribution. This leads to an overall three-peak structure in $u_{\rm tot}$. This strongly contrasts with the results in Fig. 6, where the reflection coefficients for ${\bm{c}}\parallel\hat{{\bm{z}}}$ are nearly identical to the $\xi=0$ results. Such a contrast was also present in our previous paperchurchill2016b , where a difference in transverse and longitudinal nonlocal terms had little effect on $r_{p}$, but a significant effect on $u_{\rm tot}$. The overall effect of the ABC choice remains the same as in previous sections.

From this and the previous sections, it is clear that if multiple spatially-dispersive resonances are present in a medium then they cannot be considered separately. Each of these results displays behavior not present in the single-resonance case, such as the suppression of peaks in $u_{\rm tot}$ associated with the $p$-polarization in closely-spaced parabolic bands, the difference between the heavy/light exciton model to the single band approximations and finally the additional $u_{\rm tot}$ peak from linear splitting in exciton bands.

VIII Conclusions

We have extended the work of Halevi and FuchsHalevi to derive exact expressions for electromagnetic reflection and transmission coefficients at the boundary of a medium with multiple spatially dispersive resonances in the susceptibility. Surface effects are included by using phenomenological reflection coefficients $U_{mi}$ for the polarization waves at the boundary. We have compared the results for several multi-resonance media, using a variety of $U_{mi}$ values corresponding to ABC’s in the literature. In the case of heavy/light exciton bands, we have found an improved fit for the single band approximation with $\sqrt{D^{*}}=(\sqrt{D_{\rm h}}+\sqrt{D_{\rm l}})/2$.

The model has been extended to alternate exciton dispersion relations with the inclusion of a linear splitting term in $\omega(k)$ typical of uniaxial crystals. We have compared the results when the crystal axis ${\bm{c}}$ is aligned with each of the co-ordinate axes in our system The largest effects were seen with ${\bm{c}}$ perpendicular to the plane of incidence for $s$-polarization and in the plane of incidence, parallel to the surface, for $p$-polarization.

Finally we have used the calculated reflection coefficients to find the zero-point and thermal spectral energy density $u_{\rm tot}(z,\omega)$ outside the dielectric. Many features are the same as the single-resonance model, such as the effect of the ABC choice and the saturation of $u_{\rm tot}$ as $z\to 0$. However, there is new behavior that is only present when the multiple-resonances are considered together. We have found that close resonances can lead to significant suppression in the peaks of $u_{\rm tot}(\omega)$ and that single band approximations fail to capture the correct behavior of the heavy/light exciton band model. The linear splitting term led to significant changes for $u_{\rm tot}$ in the uniaxial crystal by splitting the peak at the resonant frequency to give an overall three-peak structure.

While the model presented here incorporates many more of the features found in real materials than Halevi and Fuchs, it could be extended further to include differences between the transverse and longitudinal susceptibilitieschurchill2016b or higher-order nonlocal terms. This could also be applied to other problems, such as identifying which ABC is most appropriate for a medium with a complex exciton band structure, or in the calculation of Casimir self-forceshor14 .