Nonlocal Effective Medium (NLEM) for Quantitative Modelling of Nanoroughness in Spectroscopic Reflectance

Abstract

Spectroscopic reflectance is a versatile optical methodology for the characterization of transparent and semi-transparent thin films in terms of thickness and refractive index. The Fresnel equations are used to interpret the measurements, but their accuracy is limited when surface roughness is present. Nanoroughness can be modelled through a discretized multi-slice and effective medium approach, but to date, this offered mostly qualitative and not quantitative accuracy. Here we introduce an adaptive and nonlocal effective medium approach, which considers the relative size and environment of each discretized slice. We develop our model using finite-difference time-domain simulation results and demonstrate its ability to predict nanoroughness size and shape with relative errors < 3% in a variety of test systems. The accuracy of the model is directly compared to the prediction capabilities of the Bruggeman and Maxwell–Garnett models, highlighting its superiority. Our model is fully parametrized and ready to use for exploring the effects of roughness on reflectance without the need for costly 3D simulations and to be integrated into the Fresnel simulator of spectroscopic reflectance tools.

1. Introduction

Spectroscopic Reflectance (SR) is a simple and well-established methodology for accessing the thickness of transparent and semi-transparent films. It relies on inverting the reflectance signal using the Fresnel equations [1,2] and the known index of refraction of film and substrate. By employing SR, the thickness of films in a very wide thickness range, from 1 nm to more than 1 mm can be measured accurately and in a very short time [3,4] by employing well-calibrated and high optical resolution spectrometers. In addition, spectroscopic reflectance has been successfully employed in the measurement of the optical properties of transparent and semi-transparent films, such as real and imaginary parts of the refractive index and band gap [5]. Besides thin-film metrology [3], SR has also been demonstrated in sensing applications [6,7], whereby nanometer-scale changes in thickness due to the presence of an analyte can be easily monitored and recorded.

In the standard SR setup, the light from a light source is coupled to the illumination branch (6 fibers on a circle) of a reflection probe and illuminates at 90° the sample under characterization; the specular reflectance is collected by the receiving fiber of the reflection probe. Successful application of SR requires the thin film and substrate to have flat surfaces. For non-uniform films, the spot size, i.e., the area from which the specular signal is collected, should be quite small to exclude the macroscale roughness effect and allow for discrimination of the interference fringes. In those cases, the optical set-up is built around a microscope and the spot size is mainly defined by the magnification of the objective lens. When this is not the case, e.g., in the case of nanoscale roughness, additional film layers are required in order to model the rough thin film [8]. In the simplest approximation, an extra “roughness” layer is assumed on top of the thin film. Its refractive index is usually modelled by Maxwell–Garnett (MG) [9] or Bruggeman (BR) [10] Effective Medium Theories (EMT) assuming a volume mixture of the thin film material and air. The thickness and refractive index of the roughness layer follow from the average roughness (top-to-bottom) height and the volume coverage of the thin film material, both of which now become free parameters for fitting the SR signal. In more elaborated schemes, the roughness layer is discretized into several thinner slices, each having a different volume coverage depending on its vertical position and the actual roughness shape [8]. A Fresnel layout considering the effective multilayer system can now be easily solved. Nevertheless, we will show that using a static (i.e., with the same fixed depolarization factor [11]) EMT model in all cases [12] does not yield an optimal description, and thus the post-processing of the SR spectra will not lead to a quantitative evaluation of the film surface quality.

However, it turns out that even an EMT that varies according to the volume coverage (i.e., non-static) is not enough to represent all the roughness details, such as different aspect ratios, shapes, etc. In fact, different roughness shapes set different geometrical overlaps between successive slices, and thus slightly different light interactions, pointing towards a nonlocal effective medium model. That is, the effective medium of a slice depends not only on the microstructure of the slice itself (local model) but also on the microstructure of the neighboring (top and bottom) slices. Here, we present such a nonlocal and non-static effective medium (NLEM) model, extracted with the aid of rigorous numerical simulations, and parametrized within the context of traditional EMTs. We show that the model results in a remarkably good description of roughness (RMS and shape) effects in the reflectance spectrum. For simplicity, we assume non-dispersive non-absorbing materials and create a unified NLEM model that covers a range of film indices and roughness sizes and shapes.

The model parametrization and its relevant parameters are tabulated and ready to use, which may offer added value in roughness characterization to SR systems. The developed methodology can be used on-the-fly during numerical fitting of SR data without significant effect on the time-to-result calculation time. At the same time, however, by using the NLEM model, we find the existence of groups of roughness geometries that are optically similar (i.e., they produce the same reflectance spectrum), which precludes the simultaneous determination of both roughness height and roughness shape by a single reflectance measurement. As a result, prior knowledge of either the surface morphology or roughness height is needed to uniquely determine the roughness geometry. This is not a weakness of the NLEM model but an inherent property of nanoroughness in the reflection spectra.

2. Methodology and Results

In its simplest embodiment, thin film roughness can be modelled as a periodic square array of linear square-based pyramids, as depicted in Figure 1a. Note that diffuse scattering and/or diffraction are completely suppressed in our system because we assume the roughness as a periodic array of identical scatterers arranged at a subwavelength periodicity (i.e., below the diffraction limit $L<\Phi$ where $n$ is the maximum refractive index and Φ is the diffraction angle). While this is not realistic in the case of a rough film (where diffuse scattering is always present), it allows us to develop our model with well-controlled simulation data. The chosen pyramid period $L$ is fixed (without loss of generality) to a subwavelength value $L=32$ nm (a tenth of the smallest wavelength typically used in RS), with the half-angle $\theta$ defined by the roughness height $H$ as $\tan \theta =L/2H$.

Thus, $L$ is much smaller than the wavelength range considered (300–1200 nm) ensuring the validity of an effective medium approach. The bulk film underneath has thickness $D$, which in turn is on top of the substrate. In a SR set-up, light is incident normally from the top (i.e., from the air) and due to the absence of diffraction (subwavelength surface features), it is partially reflected vertically and partially transmitted through the substrate. For simplicity, we will assume the refractive index of the film and substrate to be constant (non-dispersive) and real (transparent). The theoretical reflection spectra, on top of which we will develop our NLEM model, are obtained using the finite-difference time-domain (FDTD) method [13] (using the Lumerical FDTD solver).

To create the NLEM model we split the pyramid into $N$ thin slices, as depicted in Figure 1b. The slice thickness $h=H/N$ (where $H$ is the pyramid thickness and N the number of slices) needs to be small enough to ensure a good description of fine pyramidal shapes and to avoid staircasing artifacts. We use $h=1 \mathrm{nm}$. The effective refractive index $\tilde{n}$ of each slice is calculated via the FDTD method using representative cells [8], allowing the construction of an effective 1D stratified model shown in Figure 1c, which can then be straightforwardly solved by the Transfer Matrix Method (TMM) [14]. For N material layers, the final transfer matrix $M$ is defined as the reverse order product:

$$M={\displaystyle \prod }_{j=N}^1 P_j{I}_{j,j-1}=P_N{I}_{N,N-1}{P}_{N-1}{I}_{N-1,N-2}\dots P_1{I}_{1,0}$$

(1)

where $I_{j,j-1}{}_{ }$ is the interface matrix between two material layers $j-1$ and $j$ with corresponding refractive indices $n_{j-1}$ and $n_j$ ($n_0=1$ for air), defined as

$$I_{j,,j-1 }=\frac{1}{n_j}\left(\begin{array}{cc}{n}_j+n_{j-1} & n_j-n_{j-1}\\ n_j-n_{j-1}& n_j+n_{j-1}\end{array}\right)$$

(2)

and $P_j$ is the propagation matrix for the material $j$

$$P_j=\left(\begin{array}{cc}{e}^{i{k}_j{d}_j}& 0\\ 0& e^{-i{k}_j{d}_j}\end{array}\right),$$

(3)

where $k_j=2\pi n_j/\lambda$ is the wavevector and $d_j$ is the thickness of the material. According to the Fresnel equations the reflectance is calculated as $R={\left|M_{21}/M_{22}\right|}^2$.

This approach was first shown in [8], where Aspnes et al. split the roughness layer into sublayers and developed the corresponding effective medium models. This method has also been used in analyzing spectroscopic ellipsometry measurements [15,16,17,18] of rough surfaces. However, in all cases, the EMT was static, without the flexibility to be adjusted to different shapes and slopes.

We use explicit FDTD simulations to extract an effective refractive index [19]. We distinguish two main representations of a slice cross-section, as shown in Figure 2: (a) the infinitesimal thickness case with a slice of thickness $h$ suspended in a vacuum, and (b) the semi-infinite thickness case, where only the top surface is probed, and multiple reflections are absent. The reason we need both representations is that case (a) alone is inadequate to properly describe columnar (high slope) roughness. We will show that the correct effective index is of the form

$$\tilde{n}\left(\lambda \right)=w{\tilde{n}}_0\left(\lambda \right)+\left(1-w\right){\tilde{n}}_{\infty }\left(\lambda \right)$$

(4)

where ${\tilde{n}}_0$, ${\tilde{n}}_{\infty }$ are the wavelength-dependent effective indices for the thin slice (case a) and semi-infinite (case b) respectively, and $w$ is a weight factor. $w$ will be shown to depend on the actual refractive index of the film material, on the relative volume fraction of the material in air, and on the slope of the pyramid walls.

In case (a), each slice $j$ (starting from the bottom) of the pyramid is treated as a square area of thickness $h$ and side length $l_j$ embedded in air, where $j=1, 2,\dots N$ is the slice index. In our linear pyramid, $l_j=L\left(N-j+0.5\right)/N$. The embedding air box (red dashed lines in Figure 2a) has the same thickness $h$ and side length L. To get the effective medium of this composite (slice in embedding air box) we perform FDTD simulation of a plane wave at normal incidence (Figure 2a) and record the complex reflection and transmition amplitudes $r\left(\lambda \right)$ and $t\left(\lambda \right)$ respectively, at the embedding box top and bottom surfaces. Then, the effective refractive index ${\tilde{n}}_0$(λ) is calculated from the relation [19]:

$$2\cos \left(2\pi {\tilde{n}}_{0}h/\lambda \right)=t-r^2{t}^{-1}+t^{-1}$$

(5)

In case (b), the complex reflection amplitude $r_{\infty }\left(\lambda \right)$ of the semi-infinite column (i.e., the bottom of the column enters into the PML boundaries [20,21]) is simulated by FDTD. For non-magnetic materials, the effective refractive index ${\tilde{n}}_{\infty }\left(\lambda \right)$ is extracted from the relation [19]:

$${\tilde{n}}_{\infty }=\left(1-r_{\infty }\right)/\left(1+r_{\infty }\right)$$

(6)

As the first test of our hypothesis of two different effective index extremes, we plot in Figure 2c the reflectance of a semi-infinite SiO₂ substrate with roughness modelled by pyramids of $H=32 \mathrm{nm}$, along with three effective refractive indices obtained by using ${\tilde{n}}_0$, ${\tilde{n}}_{\infty }$ and (${\tilde{n}}_0+{\tilde{n}}_{\infty }$)/2. A remarkable agreement is achieved with the average index, which validates our approach of combining the two extreme cases.

The above effective refractive indices were then calculated as a function of side length $l$ for four film refractive indices: $n_f=1.47$ (e.g., SiO₂), $n_f=2$ (e.g., Si₃N₄), $n_f=2.6$ (e.g., TiO₂), and $n_f=3.42$ (e.g., Si). The aim of the NLEM model is to be used on-the-fly during the numerical fitting of SR data. Thus, the FDTD-extracted effective indices need to be properly parametrized by an analytical formula. To this end we adopt the MG model [9]:

$$\frac{\tilde{\epsilon }-1}{\tilde{\epsilon }+p}=f\frac{{\epsilon }_f-1}{{\epsilon }_f+p}$$

(7)

where $\tilde{\epsilon }$ is the effective dielectric function $\left(\tilde{\epsilon }={\tilde{n}}^2\right)$, ${\epsilon }_f$ is the dielectric function of the film material, the host is assumed to be air, $f=l^2/L^2$ is the surface filling factor, and $p$ is the MG depolarization factor. The choice of the Maxwell–Garnett model over Bruggeman model is mainly motivated by its algebraic simplicity, i.e., it is first order in the effective dielectric constant, compared to second order for the latter. Here, instead of using one fixed value for $p$, we treat it as an adjustable parameter to ensure the best possible fit of Equation 7 to the FDTD effective index.

We use the filling factor $f$ as the independent variable, and find that the slice (case a) FDTD effective index is best reproduced by Equation (7) when the depolarization factor $p$ follows the parametrization:

$$p=p_0+\gamma f^{\prime }$$

(8)

where $f^{\prime }={(f^{\delta }-0.5)}^2$ and $p_0$, $\gamma$, and $\delta$ are adjustable parameters, which depend on the film refractive index $n_f$. To augment the NLEM model so that it describes all refractive indices, we repeat the above process for all modelled film indices and extract the index dependent parameters, as shown in Figure 3. These are parametrized as

$$\phi =c_0+c_1{n}_f+c_2{n}_f^2$$

(9)

where $\phi$ stands for any of the three adjustable parameters $p_0$, $\gamma$, and $\delta$. The corresponding values of the expansion coefficients are shown in

Table 1. The NLEM model parametrization.

	${\mathit{c}}_0$	${\mathit{c}}_1$	${\mathit{c}}_2$
	Depolarization factor parametrization
$p_0$	5.50	1.822	−0.2326
$\gamma$	−18.5	−1.062	-
$\delta$	0.57	0.039	−0.0053
$p_{\infty }$	1.25	-	-
	Weighting factor parametrization
a	−0.4762	0.091	-
b	0.0166	0.00041	-

. For the semi-infinite column (case b), on the other hand, we find that the best value of the depolarization factor $p$ is always approximately equal to $p_{\infty }=1.25$, irrespective of $f$ and $n_f$.

Having parametrized ${\tilde{n}}_0$ and ${\tilde{n}}_{\infty }$ through their MG representation, the next step is to find the proper weighting factors $w$. These will depend on the pyramid shape, i.e., if it is closer to a rod shape (small half-angle $\theta$) or not. To this end, we fit the NLEM model to FDTD simulations of the full pyramids of various half-angles $\theta$, using $w$ as a free parameter. This procedure works well if there are multiple reflections involved (i.e., spectrally dense interference fringes), so we adopt a simulation system involving a thin film on top of a substrate, as depicted in Figure 1a. We will simulate the layouts SiO₂/Si, SiN/Si, TiO₂/Si, and Si/SiO₂.

The error function used to fit the simulated reflectance is

$$ϵ\left(w\right)=\sum {\left[R_{FDTD}-R_{NLEM}\left(w\right)\right]}^2$$

(10)

where the reflectance $R={\left|r\right|}^2$. The sum runs over all modelled wavelengths (300–1200 nm nm spectral range with wavelength step of 1 nm), $R_{FDTD}$ is the numerically obtained reflectance, and $R_{NLEM}\left(w\right)$ the NLEM reflectance using $N=H/h$ slices in a TMM to describe the pyramids. The simulated $H$ values ranged from $H=4 \mathrm{nm}$ to $H=90 \mathrm{nm}$, while the periodicity is kept fixed at $L=32 \mathrm{nm}$. The bulk film thickness is in all cases fixed at $D=2 \mathsf{\mu }\mathrm{m}$. The weight $w$ is optimized to minimize $ϵ$ for each different $H$ and each different film index $n_f$.

The result of this optimization for the SiO₂/Si system and three values of $H$ is shown in Figure 4a. A different value of $w$ produces the minimum error between the FDTD and NLEM reflectance for different $H$ values. Specifically, $w\approx 0$ (maximal weight for ${\tilde{n}}_{\infty }$) is needed for the larger $H$ (tall columnar-like roughness) and $w\to 1$ (maximal weight for ${\tilde{n}}_0$) is needed for the smaller $H$ (short wide-angle roughness). This finding aligns well with our intuition regarding shape-dependent behavior. Thus, it calls for a nonlocal approach, which we will develop by parametrization, using the basis of the two distinct topological representations (${\tilde{n}}_0$ and ${\tilde{n}}_{\infty }$) for the effective medium.

Specifically, we proceed with repeating this exercise for all four material systems and present in Figure 4b the scaling of the weight $w$ with the pyramid half-angle $\theta$ (see Figure 4b inset). Remarkably, a simplelinear scaling law is found, which is nicely parametrized by

$$w\left(\theta \right)=a+b\theta$$

(11)

where the constants $a$, $b$ are linear functions of film’s refractive index $a,b=c_0+c_1{n}_f$, with the coefficients $c_0$, $c_1$ given in Table 1. It is interesting that we may find values of $w>1$ for the large angles (flat pyramids) and $w<0$ for the very small angles (spike-like pyramids). This points to the fact that our ${\tilde{n}}_0$ and ${\tilde{n}}_{\infty }$ do not represent the two extreme limits of the effective medium. This is particularly true for ${\tilde{n}}_0$ (large angle limit, where w > 1) and could be expected: ${\tilde{n}}_0$ is calculated for a thin slice suspended in air, while in reality the film is sandwiched between two other films. For simplicity, rather than exploring a redefinition of what the ${\tilde{n}}_0$ extreme should stand for and how it should be calculated, we acknowledge that our extrapolation with the extended $w$ values does allow NLEM to provide an excellent representation of the nonlocal effective medium, and thus continue with this formulation.

The quantitative accuracy of NLEM is shown in Figure 5a, where we plot $R_{FDTD}$ and $R_{NLEM}$ (using $h=1 \mathrm{nm}$ discretization for the film roughness) in the SiO₂/Si case for $H=32 \mathrm{nm}$. The $R_{NLEM}$ displays a remarkably accurate reproduction of $R_{FDTD}$. This NLEM accuracy also ensures a strong predictive capability. We demonstrate this in Figure 5b, where the $R_{FDTD}$ spectra from a series of FDTD run with different $H$ values are used as the “measured” signal, and optimization runs of various models (NLEM vs. other EMT models, all at $h=1 \mathrm{nm}$ discretization) try to fit $R_{FDTD}$ by using $H$ as the free parameter. The relative error in the predicted $H$ value is shown in Figure 5b as a function of the real $H$ value used in the FDTD simulation. It is impressive that the NLEM error is close to zero for all $H$ values (i.e., accurate prediction of $H$), while it smoothly varies between overestimation and underestimation for the other EMTs, crossing at some H value the zero-error line. The latter means that there is always one specific height H for each EMT where it will make a correct prediction, but as evident from Figure 5b, this is accidental. This is clear proof that the effective medium of each slice does not depend on its geometrical features only, but on the overall slope it creates with its neighbors, i.e., it is nonlocal.

We now put the NLEM model to the test by simulating different roughness shapes modelled by square-based pyramids. We parametrize the pyramid sidewalls by the exponent $s$:

$$z=H\left[1-{\left(l/L\right)}^s\right]$$

(12)

where $z$ is the side-wall height at lateral position $x=l/2$ (see inset of Figure 6), where $l$ is the local width. Equation (12) satisfies $z=H$ at $l=0$ and $z=0$ at $l=L$, while for $s=1$ it falls back to the linear pyramid studied above. By inverting Equation (12), we solve for the width l (and thus the surface coverage $f=l^2/L^2$) of a slice at vertical position z. The surface coverage for slice $j$ (starting from the bottom) is then $f={\left[\left(N-j+0.5\right)/N\right]}^{2/s}$.

We assume a fixed pyramid height of $H=32 \mathrm{nm}$ and vary the shape factor s between 0.2 and 3. The FDTD simulated reflectance is then fitted by the NLEM model with s being the free parameter. The percentage error in s prediction is shown in Figure 6 as a function of the actual s value. For comparison, the prediction errors of the other EMTs are also shown. Similarly, with the linear pyramids studied above, the NLEM error is close to zero for all $s$ values (i.e., accurate prediction of pyramid shape), while other EMTs smoothly vary between overestimation and underestimation.

The final test of our model is on its predictive capability when both roughness height $H$ and shape $s$ are treated simultaneously as free parameters. For this, we use the $R_{FDTD}$ for $H=32 \mathrm{nm}$, $s=0.5$ (see inset of Figure 7a) and color-code the NLEM error function (Equation (10)) as a function of $H$ and $s$ in Figure 7a. A unique solution (i.e., with an error significantly different than that of any other H, s combination), is not found. Instead, a continuous thin band emerges of possible $H$, $s$ combinations (of practically the same minimum error), which include the correct $H$, $s$ values. All these combinations yield the same error to $R_{FDTD}$ and thus are characterized by the same exact reflectance. That is, they correspond to different roughness geometries that produce the same exact optical effect and are thus indistinguishable based on their optical response, i.e., they correspond to optically similar roughness. To confirm this, we pick two points from the two extremes, $H=12 \mathrm{nm}$, $s=2.27$ and $H=54 \mathrm{nm}$, $s=0.28$ and plot in Figure 7b,c respectively, the color-coded error in fitting the corresponding $R_{FDTD}$. The roughness shapes are shown in the figure insets. For clarity, we use a different color-code scheme. For each case, a continuous thin band is predicted, which includes the simulated $H$, $s$.

In Figure 7d, the 3 error functions are overlayed, showing that they are on top of each other. All three $R_{FDTD}$ spectra are practically identical. This band, thus, traces between different height–shape combinations (i.e., different geometries, as shown in the insets of Figure 7) that nevertheless produce the same reflectance spectrum, i.e., it traces between optically similar roughness geometries. This is a real, inherent, property of rough surfaces, predicted by the NLEM model and validated by the FDTD calculations. While this precludes the simultaneous and unique determination of both $H$ and $s$ by a reflection spectrum, it does highlight our NLEM as a unique analytical tool in exploring subtle effects coming from variations in the roughness characteristics. No other EMT, to our knowledge, can reliably explore and identify optical similarity in roughness. NLEM eliminates the need for heavy numerical simulations for roughness. It may thus evolve into an important supplement for spectroscopic reflectance, with its usefulness being maximal when roughness variations across the wafer are small and sources of noise and errors in the measurements are low.

3. Discussion

To further show that our model correctly predicts the morphology of the thin film roughness, we applied it to two test cases using material indices and layer thicknesses other than the ones used to parametrize our model. In the first case, the substrate of the structure is Si and films of thickness 500 nm with four different values of the refractive index are placed on top of it. The relative error in the predicted $H$ and $s$ value as a function of the real $H$ and $s$ value used in the FDTD simulation are shown in Figure 8, along with the results of different effective medium models. We find that overall, the prediction errors of the NLEM are quite small and always smaller than those of other EMTs, however, at increasing index, the NLEM error also increases. This could be explained by the increased difficulty in obtaining accurate effective media when the contrast between materials (i.e., air and dielectric) increases. This however does not invalidate the usefulness of our model as it still outperforms all other traditional EMTs.

The second simulated test is a more complex structure consisting of three material layers. The top (film) and bottom (substrate) layer have an intermediate refractive index n = 3, while the middle one (spacer) has n=1.7. The thickness of the top film is 250 nm and of the spacer 500 nm. The morphology of the rough surface is assumed to be square periodic pyramids with period $L=16 \mathrm{nm}$, $H=48$ nm, and shape $s=1$ (linear pyramid). An excellent agreement between the FDTD method and the NLEM model is achieved (Figure 9), with relative errors for the predicted size $H$ and shape $s$ close to −0.5%.

These tests on independently generated simulation data ensure that the NLEM model is applicable to any rough morphology (modelled by subwavelength square-based pyramids) and non-dispersive materials with refractive indices within the range [1.47–3.47]. Note that the model and analysis presented here are based on systems simulated by FDTD, an accurate numerical method that has no sources of noise or errors. This will not be the case when NLEM is used against experimental measurements, where the roughness is not going to be uniform (in size and/or shape) across the film surface and there will be, in addition, many sources of noise and experimental errors. It is understood that NLEM’s accuracy will be compromised as these errors increase. However, so will the accuracy of any other EMT method. It is generally expected that when noise and errors become significant, a simpler EMT model with the minimum possible number of adjustable parameters would be more appropriate. However, such a systematic analysis of the sensitivity of different EMTs as a function of different sources of errors is beyond the scope here and will be the subject of follow-up work.

4. Conclusions

A nonlocal effective medium model was presented for the study of nanoroughness in thin-film metrology. It assumes roughness as a periodic array of subwavelength square-based pyramids. It decomposes each roughness pyramid into several slices and assumes an effective medium description of each slice that depends on the film index, the filling ratio, and the filling ratios of the adjacent slices, i.e., it depends on the sidewall slope. This results in an adaptive and nonlocal effective medium, which goes well beyond the state of the art which uses a fixed effective medium model for any material index, filling ratio, and environment. The excellent quantitative agreement with exact FDTD simulations for any square-based subwavelength nanoroughness size and side-wall shape, proves the usefulness of our model in reliably exploring roughness effects with a simple 1D transfer matrix method. In addition to exploring optically similar structures, the accuracy of the NLEM model also allows it to be used in an engineering design mode, where the roughness size and shape are optimized under certain application restrictions so as to generate a specific target reflectance spectrum.

The NLEM model is fully parametrized and ready to use for transparent films. The outlook is for creating a complete model which considers full material dispersion and absorption effects, and directly compares with real experimental reflectance measurements, where the NLEM findings can be validated by detailed measurements of the surface morphology (e.g., by atomic force microscopy). Such measurements will also guide the extension of NLEM to cover a wider range of pyramid base shapes and sizes beyond the square-based subwavelength sizes studied here. This is part of ongoing efforts.