Rough 1D Photonic Crystals: a transfer matrix approach

We assume that the nominal system is time invariant and has translational symmetry along the $xy$ plane, so that we can consider fields with a well defined frequency $\omega$ and 2D wavevector $\bm{Q}$ along $xy$. For local isotropic media and for a given TE or TM polarization, the transfer matrix $\mathbf{M}(z_{2},z_{1})$ is a 2$\times$2 matrix that relates the components of the electric field $E_{\parallel}$ and the magnetic field $H_{\parallel}$ parallel to the $xy$ plane, and normal to the axis of the structure which we take as the $z$ axis, evaluated at two arbitrary heights $z_{1}$ and $z_{2}$,

$$\left(\begin{array}[]{c}E_{\|}\\ H_{\|}\end{array}\right)_{z_{2}}=\mathbf{M}(z_{2},z_{1})\left(\begin{array}[]{c}E_{\|}\\ H_{\|}\end{array}\right)_{z_{1}}.$$

(1)

Several equivalent formulations have been proposed to obtain $\mathbf{M}$ [39, 40, 41]. Within a homogeneous nonmagnetic layer $\alpha$ characterized by a dielectric function $\epsilon_{\alpha}$ and refractive index $n_{\alpha}=\sqrt{\epsilon_{\alpha}}$, the field is in general the sum of upward ($+$) and downward ($-$) going plane waves, with wavevector components $k_{\alpha}^{\pm}=\pm k_{\alpha}$ along $z$ given by the dispersion relation

$$Q^{2}+k_{\alpha}^{2}=\epsilon_{\alpha}\frac{\omega^{2}}{c^{2}}.$$

(2)

For economy in notation, we will use the nomenclature of plane waves even if $Q$ is large, or $\epsilon$ is negative or complex, in which case, $k_{\alpha}$ may acquire an imaginary part and the waves become evanescent. We define $k_{\alpha}$ as the solution of (2) with a positive imaginary part. The surface impedance is defined as the quotient

$$Z=E_{\parallel}/H_{\parallel}.$$

(3)

From Faraday’s and Ampere-Maxwell’s equations we obtain

$$Z^{\pm}_{\alpha}=\pm Z_{\alpha}=\pm\begin{cases}\frac{\omega}{k_{\alpha}c},&\text{(TE)}\\ \frac{k_{\alpha}c}{\omega\epsilon_{\alpha}},&\text{(TM)}\end{cases}$$

(4)

for upwards ($+$) and downwards ($-$) moving plane waves. It is convenient to write $E_{\parallel}(z)$ and $H_{\parallel}(z)$ in terms of the electric field amplitude of the $\pm$ waves for TE polarization and in terms of the magnetic field amplitudes for TM polarization. Denoting these amplitudes as $\gamma^{\pm}(z)$ in both cases, we write

$$\left(\begin{array}[]{c}E_{\parallel}\\ H_{\parallel}\\ \end{array}\right)_{z}=\left(\begin{array}[]{cc}1&1\\ Y_{\alpha}&-Y_{\alpha}\\ \end{array}\right)\left(\begin{array}[]{c}\gamma^{+}\\ \gamma^{-}\\ \end{array}\right)_{z},\quad\text{(TE)}$$

(5)

and

$$\left(\begin{array}[]{c}E_{\parallel}\\ H_{\parallel}\\ \end{array}\right)_{z}=\left(\begin{array}[]{cc}Z_{\alpha}&-Z_{\alpha}\\ 1&1\\ \end{array}\right)\left(\begin{array}[]{c}\gamma^{+}\\ \gamma^{-}\\ \end{array}\right)_{z}.\quad\text{(TM)}$$

(6)

where $Y_{\alpha}\equiv 1/Z_{\alpha}$ is the surface admittance. As $\gamma^{\pm}(z)$ propagate as plane waves for both polarizations, $\gamma^{\pm}(z)\propto\exp(\pm ik_{\alpha}z)$ and

$$\left(\begin{array}[]{c}\gamma^{+}\\ \gamma^{-}\\ \end{array}\right)_{z_{b}}=\left(\begin{array}[]{cc}\exp(\rmi k_{\alpha}(z_{b}-z_{a}))&0\\ 0&\exp(-\rmi k_{\alpha}(z_{b}-z_{a}))\\ \end{array}\right)\left(\begin{array}[]{c}\gamma^{+}\\ \gamma^{-}\\ \end{array}\right)_{z_{a}}$$

(7)

for any $z_{a}$ and $z_{b}$ within the homogeneous layer. Evaluating (5) and (6) at both nominal edges $z^{(0)}_{\alpha-1}$ and $z^{(0)}_{\alpha}=z^{(0)}_{\alpha-1}+d_{\alpha}$ of the layer $\alpha$ of width $d_{\alpha}$ (superscript $(0)$ denotes absence of roughness) and using (7) we obtain

$$\left(\begin{array}[]{c}E_{\parallel}\\ H_{\parallel}\\ \end{array}\right)_{z^{(0)}_{\alpha}}=\mathbf{M}_{\alpha}^{(0)}\left(\begin{array}[]{c}E_{\parallel}\\ H_{\parallel}\\ \end{array}\right)_{z^{(0)}_{\alpha-1}}$$

(8)

for both TE and TM polarizations, where

$$\mathbf{M}_{\alpha}^{(0)}=\mathbf{M}_{\alpha}(d_{\alpha})=\left(\begin{array}[]{cc}\cos(k_{\alpha}d_{\alpha})&iZ_{\alpha}\sin(k_{\alpha}d_{\alpha})\\ iY_{\alpha}\sin(k_{\alpha}d_{\alpha})&\cos(k_{\alpha}d_{\alpha})\\ \end{array}\right),$$

(9)

is the transfer matrix of a homogeneous layer $\alpha$ that transfers the fields across its nominal width $d_{\alpha}$. We remark that the determinant of this matrix is $1$, as required to comply with time reversal symmetry.

The continuity of $E_{\|}$ and $H_{\|}$ imply that the fields may be transferred from the ambient $\alpha=0$ towards the substrate $\alpha=N+1$ through the layers $\alpha=1\ldots N$ by the matrix

$$\mathbf{M}^{(0)}=\mathbf{M}^{0}_{N}\mathbf{M}^{0}_{N-1}\ldots\mathbf{M}^{0}_{2}\mathbf{M}^{0}_{1}.\quad\text{(nominal)}$$

(10)

We consider a system as above, but in which a single interface $\alpha$ separating medium $\alpha$ from $\alpha+1$ is displaced from its nominal position at $z^{(0)}_{\alpha}$ to a new position $z_{\alpha}=z^{(0)}_{\alpha}+\zeta_{\alpha}$. The transfer matrix of the modified system can then be obtained by replacing $\mathbf{M}_{\alpha}^{(0)}$ and $\mathbf{M}_{\alpha+1}^{(0)}$ in (10) by $\mathbf{M}_{\alpha}(d_{\alpha}+\zeta_{\alpha})=\mathbf{M}_{\alpha}(\zeta_{\alpha})\mathbf{M}_{\alpha}^{(0)}$ and $\mathbf{M}_{\alpha+1}(d_{\alpha+1}-\zeta_{\alpha})=\mathbf{M}_{\alpha+1}^{(0)}\mathbf{M}_{\alpha+1}(-\zeta_{\alpha})$ respectively. This is equivalent to the replacement of the product $\mathbf{M}^{(0)}_{\alpha+1}\mathbf{M}^{(0)}_{\alpha}$ by $\mathbf{M}^{(0)}_{\alpha+1}\mathbf{M}^{I}_{\alpha}\mathbf{M}^{(0)}_{\alpha}$, where we introduced an effective interface transfer matrix

$$\mathbf{M}^{I}_{\alpha}=\mathbf{M}_{\alpha+1}(-\zeta_{\alpha})\mathbf{M}_{\alpha}(\zeta_{\alpha}).$$

(11)

Assuming now that all interfaces $\alpha=0\ldots N$ are shifted by a corresponding displacement $\zeta_{\alpha}$, we can build the complete transfer matrix as the product

$$\mathbf{M}=\mathbf{M}^{I}_{N}\mathbf{M}^{(0)}_{N}\mathbf{M}^{I}_{N-1}\mathbf{M}^{(0)}_{N-1}\mathbf{M}^{I}_{N-2}\ldots\mathbf{M}^{I}_{2}\mathbf{M}^{(0)}_{2}\mathbf{M}^{I}_{1}\mathbf{M}^{(0)}_{1}\mathbf{M}^{I}_{0}.\quad\text{(displaced)}$$

(12)

Now we introduce the roughness through $x$ and $y$ dependent height profiles $\zeta_{\alpha}(x,y)$, $\alpha=0\ldots N$, as illustrated in figure 1.

In the spirit of the well known Kirchhoff approximation [26, 31] in the limit of low angle roughness [42] it seems tempting to apply (12) for each $x,y$ and average it over the $xy$ plane, or, for a randomly rough system, over all realizations of an ensemble. For the case of mutually uncorrelated profiles, this would be equivalent to a replacement of $\mathbf{M}^{I}_{\alpha}$ by its average $\langle\mathbf{M}^{I}_{\alpha}\rangle$. Nevertheless, this would be wrong; the transfer matrix is not a quantity that may be meaningfully averaged, as the determinant of the average doesn’t agree in general with the average of the determinant, and the transfer matrix ought to be unimodular in order to be consistent with time inversion symmetry; even though $\det\mathbf{M}^{I}_{\alpha}=\det\mathbf{M}_{\alpha+1}(-\zeta_{\alpha})\det\mathbf{M}_{\alpha}(\zeta_{\alpha})=1$, $\det\langle\mathbf{M}^{I}_{\alpha}\rangle=\det\langle\mathbf{M}_{\alpha+1}(-\zeta_{\alpha})\mathbf{M}_{\alpha}(\zeta_{\alpha})\rangle\neq 1$. Nevertheless, the scattering matrix of each interface is a quantity that may be safely and meaningfully averaged.

We define the scattering $\mathbf{S}_{\alpha}$ matrix for the $\alpha$-th interface as

$$\mathbf{S}_{\alpha}=\left(\begin{array}[]{cc}r_{\alpha}^{+}&t_{\alpha}^{-}\\ t_{\alpha}^{+}&r_{\alpha}^{-}\end{array}\right),$$

(13)

as it produces the outgoing waves $\gamma^{-}_{\alpha}(z)=O^{-}_{\alpha}\rme^{-\rmi k_{\alpha}(z-z_{\alpha}^{(0)})}$ and $\gamma^{+}_{\alpha+1}(z)=O^{+}_{\alpha}\rme^{\rmi k_{\alpha+1}(z-z_{\alpha}^{(0)})}$ when applied to the incoming waves $\gamma^{+}_{\alpha}(z)=I^{+}_{\alpha}\rme^{\rmi k_{\alpha}(z-z_{\alpha}^{(0)})}$ and $\gamma^{-}_{\alpha+1}(z)=I^{-}_{\alpha}\rme^{-\rmi k_{\alpha+1}(z-z_{\alpha}^{(0)})}$, namely,

$$\left(\begin{array}[]{c}O^{-}_{\alpha}\\ O^{+}_{\alpha}\end{array}\right)=\mathbf{S}_{\alpha}\left(\begin{array}[]{c}I^{+}_{\alpha}\\ I^{-}_{\alpha}\end{array}\right),$$

(14)

where we defined the amplitudes $O^{\pm}_{\alpha}$ and $I^{\pm}_{\alpha}$ using the nominal plane $z=z_{\alpha}^{(0)}$ as a reference, and where $r^{\pm}_{\alpha}$ and $t^{\pm}_{\alpha}$ are the reflection and transmission amplitudes corresponding to waves that impinge on the interface $\alpha$ moving upwards ($+$) within layer $\alpha$ or moving downwards ($-$) within layer $\alpha+1$. We may characterize the fields above and below the interface in terms of the components of $\mathbf{S}_{\alpha}$ through the matrices

$$\mathbf{A}_{\alpha}=\left(\begin{array}[]{cc}t^{+}_{\alpha}&1+r^{-}_{\alpha}\\ Y_{\alpha+1}t^{+}_{\alpha}&-Y_{\alpha+1}(1-r^{-}_{\alpha})\end{array}\right)\quad\text{(TE)}$$

(15)

and

$$\mathbf{B}_{\alpha}=\left(\begin{array}[]{cc}1+r^{+}_{\alpha}&t^{-}_{\alpha}\\ Y_{\alpha}(1-r^{+}_{\alpha})&-Y_{\alpha}t^{-}_{\alpha}\end{array}\right)\quad\text{(TE)}$$

(16)

for TE polarization, and

$$\mathbf{A}_{\alpha}=\left(\begin{array}[]{cc}Z_{\alpha+1}t^{+}_{\alpha}&-Z_{\alpha+1}(1-r^{-}_{\alpha})\\ t^{+}_{\alpha}&1+r^{-}_{\alpha}\end{array}\right)\quad\text{(TM)}$$

(17)

and

$$\mathbf{B}_{\alpha}=\left(\begin{array}[]{cc}Z_{\alpha}(1-r^{+}_{\alpha})&-Z_{\alpha}t^{-}_{\alpha}\\ 1+r^{+}_{\alpha}&t^{-}_{\alpha}\end{array}\right)\quad\text{(TM)}$$

(18)

for TM polarization. These matrices yield the electromagnetic field extrapolated to each side, above and below the nominal surface when they act on the incoming amplitudes $(I_{\alpha}^{+},I_{\alpha}^{-})^{T}$. Thus, we may relate $\mathbf{A}_{\alpha}$ and $\mathbf{B}_{\alpha}$ through $\mathbf{M}^{I}_{\alpha}$ and write

$$\mathbf{M}_{\alpha}^{I}=\mathbf{A}_{\alpha}\mathbf{B}_{\alpha}^{-1}.$$

(19)

For a flat displaced interface we can obtain the elements of the scattering matrix either from Eqs. (11) and (19) or through a simple extrapolation of the fields from the actual towards the nominal interface,

$$\mathbf{S}_{\alpha}=\left(\begin{array}[]{cc}r^{(0)+}_{\alpha}\rme^{2\rmi k_{\alpha}\zeta_{\alpha}}&t^{(0)-}_{\alpha}\rme^{-\rmi\Delta k_{\alpha}\zeta_{\alpha}}\\ t^{(0)+}_{\alpha}\rme^{-\rmi\Delta k_{\alpha}\zeta_{\alpha}}&r^{(0)-}_{\alpha}\rme^{-2\rmi k_{\alpha+1}\zeta_{\alpha}}\end{array}\right)$$

(20)

where $r^{(0)\pm}_{\alpha}$ and $t^{(0)\pm}_{\alpha}$ are the Fresnel coefficients of the nominal interface [43] and $\Delta k_{\alpha}=k_{\alpha+1}-k_{\alpha}$. Thus, not unexpectedly, the effect of a rigid shift of an interface is to incorporate simple phase factors proportional to $\zeta_{\alpha}$ into the optical coefficients. The interfacial matrix (11) can be recovered in this case from Eqs. (15)-(20).

Going back to a rough interfacce, in the spirit of the Kirchhoff approximation for low angle surfaces [26, 27], we average (20),

$$\langle\mathbf{S}_{\alpha}\rangle\equiv\left(\begin{array}[]{cc}\Braket{r^{+}_{\alpha}}&\Braket{t^{-}_{\alpha}}\\ \Braket{t^{+}_{\alpha}}&\Braket{r^{-}_{\alpha}}\end{array}\right)=\left(\begin{array}[]{cc}r^{(0)+}_{\alpha}\Braket{\rme^{2\rmi k_{\alpha}\zeta_{\alpha}}}&t^{(0)-}_{\alpha}\Braket{\rme^{-\rmi\Delta k_{\alpha}\zeta_{\alpha}}}\\ t^{(0)+}_{\alpha}\Braket{\rme^{-\rmi\Delta k_{\alpha}\zeta_{\alpha}}}&r^{(0)-}_{\alpha}\Braket{\rme^{-2\rmi k_{\alpha+1}\zeta_{\alpha}}}\end{array}\right),$$

(21)

we replace the optical coefficients by their averages in Eqs. (15)-(18) to obtain $\langle\mathbf{B}_{\alpha}\rangle$ and $\langle\mathbf{A}_{\alpha}\rangle$ and transferring the averaged fields across the nominal interface we obtain, in analogy with (19), the macroscopic interface transfer matrix $\mathbf{M}^{MI}_{\alpha}$ as

$$\mathbf{M}_{\alpha}^{MI}=\langle\mathbf{A}_{\alpha}\rangle\langle\mathbf{B}_{\alpha}\rangle^{-1}.$$

(22)

From Eqs. (15), (16) and (21), we notice that for TE polarization

$$\det\langle\mathbf{A}_{\alpha}\rangle=-2Y_{\alpha+1}t_{\alpha}^{(0)+}\langle\rme^{-\rmi\Delta k_{\alpha}\zeta_{\alpha}}\rangle$$

(23)

and

$$\det\langle\mathbf{B}_{\alpha}\rangle=-2Y_{\alpha}t_{\alpha}^{(0)-}\langle\rme^{-\rmi\Delta k_{\alpha}\zeta_{\alpha}}\rangle.$$

(24)

Nevertheless, from Fresnel’s formulae we obtain

$$Y_{\alpha}t_{\alpha}^{(0)-}=Y_{\alpha+1}t_{\alpha}^{(0)+},$$

(25)

so that $\det\mathbf{A}_{\alpha}=\det\mathbf{B}_{\alpha}$ and from (22) we verify the required unimodularity condition

$$\det\mathbf{M}^{MI}_{\alpha}=1.$$

(26)

It is trivially verified that this condition also holds for TM polarization.

Finally, we get a macroscopic unimodular transfer matrix $\mathbf{M}^{M}$for the complete system by replacing all the interface transfer matrices in (12) by their macroscopic counterparts (22),

$$\mathbf{M}^{M}=\mathbf{M}^{MI}_{N}\mathbf{M}^{(0)}_{N}\mathbf{M}^{MI}_{N-1}\mathbf{M}^{(0)}_{N-1}\mathbf{M}^{MI}_{N-2}\ldots\mathbf{M}^{MI}_{2}\mathbf{M}^{(0)}_{2}\mathbf{M}^{MI}_{1}\mathbf{M}^{(0)}_{1}\mathbf{M}^{MI}_{0}.$$

(27)

We remark that although our formulation is close to that of [28], their $\mathbf{W}$ matrix is not unimodular.

From the transfer matrix of the system we may obtain its optical properties straightforwardly using standard procedures. For example, the reflection and transmission amplitudes, $r$ and $t$, when the system is illuminated from the ambient towards the substrate, may be obtained by solving the $2\times 2$ system of equations

$$\left(\begin{array}[]{c}t\\ Y_{N+1}t\\ \end{array}\right)=\mathbf{M}^{M}\left(\begin{array}[]{c}1+r\\ Y_{0}(1-r)\\ \end{array}\right),\quad(\mathrm{TE})$$

(28)

and

$$\left(\begin{array}[]{c}Z_{N+1}t\\ t\\ \end{array}\right)=\mathbf{M}^{M}\left(\begin{array}[]{c}Z_{0}(1-r)\\ 1+r\\ \end{array}\right).\quad(\mathrm{TM})$$

(29)

Naturally, we define the nominal positions $z^{(0)}_{\alpha}$ by demanding $\langle\zeta_{\alpha}\rangle=0$. Thus, in the limit of small height roughness we may approximate (21) to order $\zeta_{\alpha}^{2}$ by

$$\langle\mathbf{S}_{\alpha}\rangle=\left(\begin{array}[]{cc}r^{(0)+}_{\alpha}(1-2k_{\alpha}^{2}\tilde{\zeta}_{\alpha}^{2})&t^{(0)-}_{\alpha}(1-\frac{1}{2}(\Delta k_{\alpha})^{2}\tilde{\zeta}_{\alpha}^{2})\\ t^{(0)+}_{\alpha}(1-\frac{1}{2}(\Delta k_{\alpha})^{2}\tilde{\zeta}_{\alpha}^{2})&r^{(0)-}_{\alpha}(1-2k_{\alpha+1}^{2}\tilde{\zeta}_{\alpha}^{2})\end{array}\right),$$

(30)

where $\tilde{\zeta}_{\alpha}\equiv\langle\zeta_{\alpha}^{2}\rangle^{1/2}$ is the RMS height, so that all optical coefficients are reduced by a factor of order $\tilde{\zeta}_{\alpha}^{2}/\lambda^{2}$, with $\lambda$ the free space wavelength. This wavelength-dependent reduction is due to the energy that is lost through scattering out of the specular direction and plays a role analogous but not identical to dissipation [44]. On the other hand, if $\zeta_{\alpha}$ is not necessarily small but obeys a Gaussian distribution with zero mean, we may evaluate (21) to obtain

$$\langle\mathbf{S}_{\alpha}\rangle=\left(\begin{array}[]{cc}r^{(0)+}_{\alpha}\rme^{-2k_{\alpha}^{2}\tilde{\zeta}_{\alpha}^{2}}&t^{(0)-}_{\alpha}\rme^{-(\Delta k_{\alpha})^{2}\tilde{\zeta}_{\alpha}^{2}/2}\\ t^{(0)+}_{\alpha}\rme^{-(\Delta k_{\alpha})^{2}\tilde{\zeta}_{\alpha}^{2}/2}&r^{(0)-}_{\alpha}\rme^{-2k_{\alpha+1}^{2}\tilde{\zeta}_{\alpha}^{2}}\end{array}\right).$$

(31)

A simple substitution of any of these into (22) yields the corresponding interface matrix $\mathbf{M}^{MI}_{\alpha}$. Thus, in the small roughness case (30) we obtain

$$\mathbf{M}_{\alpha}^{MI}=\left(\begin{array}[]{cc}1+\frac{(Z_{\alpha+1}-Z_{\alpha})(k_{\alpha}+k_{\alpha+1})^{2}}{2(Z_{\alpha}+Z_{\alpha+1})}\tilde{\zeta}_{\alpha}^{2}&-\frac{(Z_{\alpha}k_{\alpha}-Z_{\alpha+1}k_{\alpha+1})^{2}}{(Z_{\alpha}+Z_{\alpha+1})}\tilde{\zeta}_{\alpha}^{2}\\ -\frac{(Z_{\alpha+1}k_{\alpha}-Z_{\alpha}k_{\alpha+1})^{2}}{(Z_{\alpha}+Z_{\alpha+1})Z_{\alpha}Z_{\alpha+1}}\tilde{\zeta}_{\alpha}^{2}&1-\frac{(Z_{\alpha+1}-Z_{\alpha})(k_{\alpha}+k_{\alpha+1})^{2}}{2(Z_{\alpha}+Z_{\alpha+1})}\tilde{\zeta}_{\alpha}^{2}\\ \end{array}\right),$$

(32)

and in the general case, we obtain

	$\displaystyle\mathbf{M}_{\alpha}^{MI}=$	$\displaystyle\frac{\braket{t_{\alpha}^{+}}}{2}\left(\begin{array}[]{cc}1&Z_{\alpha}\\ Y_{\alpha+1}&Z_{\alpha}Y_{\alpha+1}\\ \end{array}\right)$		(35)
		$\displaystyle+\frac{1}{2\braket{t_{\alpha}^{-}}}\left(\begin{array}[]{cc}(1+\braket{r_{\alpha}^{-}})(1-\braket{r_{\alpha}^{+}})&-(1+\braket{r_{\alpha}^{-}})(1+\braket{r_{\alpha}^{+}})Z_{\alpha}\\ -(1-\braket{r_{\alpha}^{-}})(1-\braket{r_{\alpha}^{+}})Y_{\alpha+1}&(1-\braket{r_{\alpha}^{-}})(1+\braket{r_{\alpha}^{+}})Z_{\alpha}Y_{\alpha+1}\\ \end{array}\right).$		(38)

For a Gaussian roughness we can further substitute (31).

In summary, we have developed a formalism that allows us to calculate the macroscopic effective transfer matrix of a multilayered system with rough interfaces using the Kirchhoff small-angle approximation. The transfer matrix is the product of the usual transfer matrices for each layer, corresponding to their nominal width, alternating with interface transfer matrices which can be obtained from the average of the corresponding scattering matrices. We obtained explicit expressions in terms of the variance of the height of the interface for the case of small height roughness and for the case of Gaussian roughness. In the following sections, we will illustrate the use of these matrices to calculate the optical properties of some rough stratified systems and we will compare some of our results to experiment. In A we extend the above theory to account for the possibility of several mutually correlated interfaces.

In figure 4 we show the dispersion relation $\omega$ vs. $K$ of the Bloch waves of an infinite 1D-PhC, calculated as discussed above for parallel wavevector $Q=0$.

The system consists of two alternating films with thicknesses $d_{1}$=100nm and $d_{2}$=80nm. To elucidate the role of roughness and of dissipation, we calculated the modes for three cases: dispersionless, dissipationless media with flat surfaces ($\epsilon_{1}=12$, $\epsilon_{2}=1$, $\tilde{\zeta}_{1}=\tilde{\zeta}_{2}=0$), the same system with some dissipation artificially added ($\epsilon_{1}=12+1.2i$, $\epsilon_{2}=1+0.1i$, $\tilde{\zeta}_{1}=\tilde{\zeta}_{2}=0$), and the same system without dissipation but with some Gaussian roughness ($\epsilon_{1}=12$, $\epsilon_{2}=1$, $\tilde{\zeta}_{1}=8$nm, $\tilde{\zeta}_{2}=10$nm). For the dissipative case we added an imaginary part to $\epsilon_{\alpha}$ equal to 10% of its real part. For the rough case, we chose the roughness amplitude as 10% of the width of the corresponding layer. Notice that this roughness is not small for the highest frequencies in the figure. Thus, (32) is not applicable, but we use (35). For the flat system without dissipation (figure 4, left) we found, as expected, bands for which the modes may propagate, with $\text{Im}\,K=0$ and with alternating positive and negative group velocity, and gaps for which propagation is forbidden and the modes decay exponentially, with $\text{Im}\,K\neq 0$ and $\text{Re}\,K=0$ or $\pm\pi/L$, at the center or at the edges of the Brillouin zone. The center of the gaps are at frequencies $\omega_{m}$ whose corresponding free-space wavelength $\lambda_{m}=2\pi c/\omega_{m}$ is a sub-multiple of twice the optical width of a period, yielding constructive interference of the waves back-scattered by successive periods,

$$n_{1}d_{1}+n_{2}d_{2}=m\lambda_{m}/2,$$

(45)

and where $\text{Im}\,K$ peaks. In the case with dissipation (figure 4, center), $\text{Im}\,K$ also displays maxima corresponding to the gaps, but $\text{Re}K$ for those frequencies is no longer a constant at the center or at the edges of the Brillouin zone, corresponding to both decay and propagation. There are also regions where $\text{Im}K$ becomes relatively small, corresponding to the propagating bands, but it is not zero, so there is some decay due to extinction through absorption. For high enough frequencies the minima of $\text{Im}K$ is almost $1/L$, so that the decay distance becomes of the order of one period. Of course, these are not actual modes that may be excited within an infinite system, but they may be excited in truncated photonic crystals or in crystals with defects. In the case with roughness (figure 4, right), the results are very similar to those obtained in the case with dissipation, as roughness also yields extinction, though its origin is scattering out of the specular direction instead of absorption. We remark that we didn’t calculate explicitly the non-specular scattered fields, but the optical theorem implies that the energy they carry away manifests itself through destructive interference in the specular term [29].

In order to explore roughness effects for realistic systems, in figure 5 we show the band structure for a nanoporous anodic silicon (NAS) 1D-PhC made of two alternating layers with porosities $p_{1}=73$% and $p_{2}=54$% and widths $d_{1}=111$nm and $d_{2}=74$ nm, as obtained from [52]. We obtained $\epsilon_{\alpha}$ from the porosity $p_{\alpha}$ and the dielectric function of silicon, taken from Palik’s handbook [48], using the Bruggeman 2D model [47]. In figure 5 we show results corresponding to a flat structure and to a structure with rough surfaces considering two cases, a RMS height of 10% of the corresponding layer’s width, i.e., $\tilde{\zeta}_{1}=11.1$ nm and $\tilde{\zeta}_{2}=7.4$nm, and a RMS height of 20% of the layer’s width ($\tilde{\zeta}_{1}=22.2$ nm and $\tilde{\zeta}_{2}=14.8$nm). As expected, the middle of the first band gap at $\hbar\omega_{1}=2.24$eV complies with (45), with $m=1$. For this system, the effects of dissipation dominate those of roughness. For the 10% case the roughness corrections are barely visible, though for the 20% they are relatively small but clearly discernible. We expect roughness to be more important for materials with a smaller dissipation, such as NAA.

In figure 6 we show the band diagram of an NAA 1D-PhC made up two alternating layers of widths $d_{1}=64$ nm and $d_{2}=130$ nm and porosities $p_{1}=7\%$ and $p_{2}=35\%$. We used Bruggeman’s 2D model in order to calculate the dielectric functions $\epsilon_{\alpha}$ ($\alpha=1,2$) of the porous layers from that of Al [51, 48]. The calculation was done both for the case of flat interfaces but including dissipation within the alumina due to impurities of $\text{Fe}_{2}\text{O}_{3}$ and CuO, as in Sec. 3, and for that of rough interfaces characterized by $\tilde{\zeta}_{1}=6.4$ nm and $\tilde{\zeta}_{2}=13$nm (so that $\tilde{\zeta}_{\alpha}/d_{\alpha}=0.1$) and without dissipation within the alumina. We plot the imaginary part of the Bloch vector $K$. Two band gaps are clearly seen in the band structure, consistent with (45), in which the imaginary part of Bloch’s vector $K$ has large peaks (see inset of figure 6) with or without roughness of the order of $0.1/L$, corresponding to small penetration depths of $\approx 10L$. Besides the gaps, there are propagation bands for which $\text{Im}\,K$ would be null in the absence of roughness or dissipation. Nevertheless, $\text{Im}K\neq 0$ for the rough system and takes values as high as $K\approx 0.005/L$ corresponding to a penetration distance of at most $\approx 200$ periods. For the parameters chosen, the effects of roughness are qualitatively similar to those of dissipation.

Now, we consider a finite stratified system of width $NL$ made from $N$ repetitions of a unit cell composed of two layers of nominal thicknesses $d_{1}$ and $d_{2}=L-d_{1}$, refractive indices $n_{1}$ and $n_{2}$, and with mutually uncorrelated rough surfaces characterized by the RMS heights $\tilde{\zeta}_{1}$ and $\tilde{\zeta}_{2}$, respectively. The refractive index of ambient is $n_{0}$ while for the substrate is $n_{2N+1}=n_{s}$. We can obtain the reflectance and transmittance of the system from Eqs. (27)-(29).

In figure 7 we show the normal incidence reflectance $R$ and transmittance $T$ spectra of a free standing finite 1D-PhC in air ($n_{0}=n_{s}=1$) consisting of $N=10$ periods of the same photonic crystals as in figure 4, with and without dissipation and roughness. We also present results for a system with flat surfaces and no roughness nor dissipation, but whose thicknesses have mutually uncorrelated fluctuations, and for a system without dissipation nor width fluctuations but with a roughness whose amplitudes grow linearly from the surface towards the substrate. As shown in the top row of figure 7, in the case without dissipation, roughness nor fluctuations, there are relatively wide regions of very high reflectance $R\approx 1$ and very small transmittance $T\approx 0$ corresponding to the band gaps shown in figure 4. Outside of these regions, the reflectance and transmittance show strong oscillations due to the interference of the waves reflected from the front and back surfaces of the system. These oscillations take $R$ to values as low as 0. Correspondingly, $T$ takes values as high as 1.

In the second row of figure 7, we show the reflectance and transmittance averaged over an ensemble of 1D-PhC’s as those corresponding to the first row, but with fluctuations in the widths of each film obeying a Gaussian distribution with standard deviations $\sigma_{1}=3$nm and $\sigma_{2}=2.4$nm corresponding to $3\%$ of their nominal widths. We used one thousand ensemble members for our calculation. The purpose of this calculation is to allow for a finite transverse coherence length $\xi$ [29] for the incident light, so that the contributions to the reflectance and transmittance from illuminated regions farther apart than $\xi$ add incoherently. This situation will prove relevant for the system studied in the next section.

The results show regions where $R\approx 1$ and $T\approx 0$ as in the top row, but the oscillating structure due to the interference is damped by the averaging process, as the frequencies corresponding to each maxima and minima differ for each member of the ensemble. This damping is larger at higher frequencies, for which the oscillations displayed by the first row become narrower and their frequency spacing becomes smaller. The third row of figure 7 shows the effects of roughness characterized by the amplitudes $\tilde{\zeta}_{1}=\tilde{\zeta}_{2}=10$nm. In this case, $R$ shows maxima corresponding to the band gaps of figure 4, but they don’t attain the value 1, as in the top and second rows. The reason is that in this case, energy is lost through scattering, reducing the specular reflectance. In our treatment of roughness, we averaged scattering matrices, thus incorporating phase information that is not present in our treatment of thickness fluctuations corresponding to the second row of figure 7, for which we averaged the reflectance and transmittance. Hence the qualitative differences. Corresponding to the maxima in $R$ for the case of a rough surface there are minima in $T$. Nevertheless, most of the energy is reflected or scattered away so that beyond $\omega\approx 0.3\times 2\pi c/L$ the transmittance becomes negligible. The third row of figure 7 also shows results for a flat absorptive system with $\epsilon_{1}=12+1.2i$ and $\epsilon_{2}=1+0.1i$. The result of dissipation is similar to that of roughness without dissipation.

In the bottom row of figure 7 we show the effect of a roughness that increases linearly in height from the front ($\tilde{\zeta}=0$) towards the back ($\tilde{\zeta}=20$nm) so that its average coincides with that considered in the third row. In this case, the reflectance does attain the value $1$ corresponding to the band gaps, as in the first row, since within the gaps most of the electromagnetic energy is reflected before reaching the deep layers and thus it doesn’t sense the larger roughness. Nevertheless, the interference oscillations of the first row corresponding to the propagating bands are smoothed out as in the second and third rows. On the other hand, the transmittance is very similar to that of the third row for all frequencies as the transmitted energy is scattered away by the rough deep layers. If we invert the roughness progression, the results would be similar to those in the third row.

In the next section we will explore a real stratified system for which the analysis of roughness that increases with depth plays a relevant role.