Optical Characterization of Non-Stoichiometric Silicon Nitride Films Exhibiting Combined Defects

Abstract

The study was devoted to optical characterization of non-stoichiometric silicon nitride films prepared by reactive magnetron sputtering in argon-nitrogen atmosphere onto cold (unheated) substrates. It was found that these films exhibit the combination of three defects: optical inhomogeneity (refractive index profile across the films), uniaxial anisotropy with the optical axis perpendicular to the boundaries and random roughness of the upper boundaries. The influence of the uniaxial anisotropy was included into the corresponding formulae of the optical quantities using the matrix formalism and the approximation of the inhomogeneous layer by a multilayer system consisting of large number thin homogeneous layers. The random roughness was described using the scalar diffraction theory. The processing of the experimental data was performed using the multi-sample modification of the least-squares method, in which experimental data of several samples differing in thickness were processed simultaneously. The dielectric response of the silicon nitride films was modeled using the modification of the universal dispersion model, which takes into account absorption processes corresponding to valence-to-conduction band electron transitions, excitonic effects and Urbach tail. The spectroscopic reflectometric and ellipsometric measurements were supplemented by measuring the uniformity of the samples using imaging spectroscopic reflectometry.

1. Introduction

Silicon nitride thin films are intensively utilized in practice. They are employed as antireflection coatings in silicon solar cells because of their suitable optical properties owing to the silicon optical properties (see, e.g., [1,2,3,4]). Silicon nitride films also exhibit appropriate properties of passivation on silicon surfaces (see, e.g., [5,6]). They are utilized as passivation layers and masks for local oxidation of silicon surfaces and optical waveguides. Non-stoichiometric silicon nitride (SiN${}_x$) films have been employed as inhomogeneous films for creating optical filters [7]. The SiN${}_x$ films, which are silicon-rich, have utilization in many other applications. For example, these films are utilized in memory devices [8], photonic devices [9] and membranes for micromechanical devices [10]. For many practical applications, it is important to know the optical properties of stoichiometric (Si${}_3$N${}_4$) and non-stoichiometric SiN${}_x$ films. Therefore, considerable attention has been devoted to the optical characterization of these films so far.

From the papers published, one can see that the optical constants of the SiN${}_x$ are strongly dependent on technological conditions of their preparation. This statement is true for the SiN${}_x$ rich in silicon in particular. Thus, individual methods or approaches must be used to characterize thin films of SiN${}_x$. Correct results of the optical characterization of these films can only be achieved if true dispersion and structural models are utilized. Below, the typical results of the optical characterization of the SiN${}_x$ films presented in the literature are introduced.

In [11], spectroscopic ellipsometry was employed for the optical characterization of the SiN${}_x$ films deposited onto silicon single substrates in spectral range 1.5–4.0 eV. An evolution of spectral dependencies of their refractive index was investigated in dependency on stoichiometric ratio of Si and Si${}_3$N${}_4$. It was stated that the refractive index is a precise indicator of the composition of the studied SiN${}_x$ films. For determining the refractive index spectral dependencies, the Farouhi–Bloomer dispersion model was used. This model is problematic from the physical point of view since it lacks the time reversal symmetry and the integral in the sum rule is divergent. Moreover, this model is not Kramers–Kronig consistent if the zero absorption below the band gap is considered. In the structural model of the SiN${}_x$ films investigated in [11], no defects were taken into account. Note that the films were prepared by magnetron sputtering.

In [12], the spectral dependencies of transmittance of the SiN${}_x$ films prepared onto quartz substrates were measured within the spectral range 170–800 nm. Experimental data of transmittance were not processed by any dispersion model. Only the values of the band gap were determined using the Tauc relation (for details, see [12]). The dependency of the band gap values on a negative substrate bias was determined (films were prepared using radio frequency magnetron sputtering). The structural model of the SiN${}_x$ films was not mentioned either. By means of an atomic force microscope, it was found that the upper boundaries of the films exhibited very fine random roughness in values of several Angstroms.

In [13], the optical characterization of the SiN${}_x$ prepared by plasma-enhanced chemical vapor deposition (PECVD) was performed by the processing experimental data obtained by spectroscopic ellipsometry applied at incidence angle of 65${}^{\circ }$ within the spectral range 245–800 nm. Experimental data achieved using variable angle of incidence ellipsometry applied for several wavelengths from the visible range together with reflectance measured for wavelength of 248 nm at normal incidence were also used for the optical characterization of the SiN${}_x$ films. The processing of the experimental data was based on expression of the optical constants of these films by means of effective medium approximation (EMA). By means of the EMA, an isotropic mixture consisting of three components, i.e., amorphous Si${}_3$N${}_4$, amorphous Si and voids, was created. This mixture served as an effective material of the SiN${}_x$ films studied. In this paper, the spectral dependencies of the optical constants of these components are considered to be known. This means that no dispersion model was assumed. Then, the values of thickness, optical constants and relative volume of amorphous silicon in SiN${}_x$ films were determined for different technological conditions. A structural model utilized at the optical characterization of these films was not presented either.

In [14], the realistic Tauc–Lorentz dispersion model was used to fit data of spectroscopic ellipsometry measured for SiN${}_x$ films prepared by the PECVD within the spectral range 1–5.5 eV. Random roughness of the upper boundaries of these films was modeled by using the Bruggeman formula of the EMA. It was assumed that 50% voids and 50% SiN${}_x$ form the effective layer representing this roughness. The thickness of this effective layer together with spectral dependencies of the optical constants and thickness of the SiN${}_x$ films were determined. It should be emphasized that the description of random roughness by the Bruggeman formula is not too suitable since the effective layer do not substitute this roughness in realistic way. The use of some theories describing an interaction of light and random roughness more exactly than the EMA such as the Rayleigh–Rice theory or scalar diffraction theory is much better from the physical point of view (see, e.g., [15,16,17,18,19,20,21,22,23]). Other defects were not considered in [14].

In this study, the optical characterization of the samples of the amorphous SiN${}_x$ films prepared by reactive magnetron sputtering onto silicon single substrates was performed. These films exhibit three defects, i.e., inhomogeneity represented by profiles of the optical constants across the films, artificial optical anisotropy in the optical constants corresponding to uniaxial anisotropy with the optical axis perpendicular to the boundaries and random roughness of the upper boundaries. A new, very complicated and efficient method of the optical characterization of these films had to be developed. The application of this method for practical purposes is illustrated by means of characterizing three samples of these films. On the basis of the results achieved for these samples, the practical usefulness of the method is presented.

2. Theory

2.1. Structural Model

By means of X-ray diffraction measurements, it was found that the non-stoichiometric silicon nitride films (SiN${}_x$) studied here were amorphous. The simultaneous processing of experimental data obtained by variable angle spectroscopic ellipsometry and spectroscopic reflectometry showed that formulae for the optical quantities corresponding to the thin film model without any defects could not be utilized. It was evident because fits achieved by this simultaneous processing of ellipsometric and reflectometric data were very poor. Differences between the experimental and theoretical data were too large. However, if the ellipsometric and reflectometric data were processed separately, better fits were achieved for each of these data. This indicates that the characterized films exhibited some defects. To our experience, inhomogeneity of the optical constants of the SiN${}_x$ films often occurs. Therefore, the inhomogeneity in the optical constants across the films was assumed. In other words, it was assumed that the refractive index and extinction coefficient changed in the direction perpendicular to the boundaries. The reality that the fits of the individual data processed separately were evidently better than the fits obtained by the simultaneous processing of the ellipsometric and reflectometric data indicated that the SiN${}_x$ films could exhibit an optical anisotropy. It is known that this optical anisotropy originates if mechanical stresses occur inside the films. This anisotropy mostly corresponds to the uniaxial anisotropy with the optical axis perpendicular to the film boundaries. Therefore, such the kind of the optical anisotropy was included in our structural model of the SiN${}_x$ films studied here. By using atomic force microscopy (AFM), random roughness of the upper boundaries of these films was revealed. The AFM image of the upper boundary of one selected sample proving this roughness is presented in Figure 8. After including all the three defects into the structural model of the studied SiN${}_x$ films, good fits of the ellipsometric and reflectometric data were achieved when their simultaneous processing was performed (see below).

Thus, the structural model of the SiN${}_x$ films utilized in this work is specified as follows:

• The ambient is formed by vacuum with the refractive index equal to unity.
• The substrate is silicon single crystal wafer, i.e., the substrate is optically homogeneous and isotropic.
• The films are uniform in thickness inside the areas used for spectrophotometric and ellipsometric measurements.
• The non-stoichiometric silicon nitride thin films are inhomogeneous with profiles of the optical constants across these films.
• These films exhibit uniaxial anisotropy with the optical axis perpendicular to the boundaries.
• Random roughness was assumed on the upper boundaries of the films. It was assumed that this roughness corresponds to the Gaussian stochastic process.

2.2. Optical Quantities of Inhomogeneous Films

The method that was used to calculate the Fresnel coefficients of the inhomogeneous film with uniaxial anisotropy having the optical axis perpendicular to the boundaries is described in this section. The dielectric response at depth z inside the inhomogeneous film is described by the dielectric tensor

$$\widehat{\mathit{\epsilon }}\left(z\right)=\left(\begin{array}{ccc}{\widehat{\epsilon }}_{\mathrm{o}}\left(z\right)& 0& 0\\ 0& {\widehat{\epsilon }}_{\mathrm{o}}\left(z\right)& 0\\ 0& 0& {\widehat{\epsilon }}_{\mathrm{e}}\left(z\right)\end{array}\right),$$

(1)

where ${\widehat{\epsilon }}_{\mathrm{o}}\left(z\right)$ and ${\widehat{\epsilon }}_{\mathrm{e}}\left(z\right)$ are the ordinary and extraordinary complex dielectric functions.

The Fresnel coefficients were calculated using a method in which the inhomogeneous film was approximated by a multilayer system consisting of large number of thin homogeneous layers. The inhomogeneous film was divided by N boundaries into $N+1$ layers in a manner depicted in Figure 1. The thickness of the inner layers is $h=d/N$, where d is the total thickness of the inhomogeneous film, while the first and the last layer have half the thickness of the inner layers. The ordinary and extraordinary dielectric functions of the individual layers were calculated as

$${\widehat{\epsilon }}_{\mathrm{o},k}={\widehat{\epsilon }}_{\mathrm{o}}\left(z_k\right),{\widehat{\epsilon }}_{\mathrm{e},k}={\widehat{\epsilon }}_{\mathrm{e}}\left(z_k\right),\mathrm{where}{z}_k=k\frac{d}{N},$$

(2)

and $k=0,1,2,\dots ,N$.

The Fresnel coefficients of the layered systems can be calculated very efficiently using the matrix formalism. We present the method based on the Yeh matrix formalism [24], which was introduced to deal with layered systems containing anisotropic media. However, it should be noted that there are also other matrix formalisms usable for the special type of anisotropy considered in this paper, e.g., the transfer matrix formalism [25,26].

The junction conditions imply that the components of the wavevectors parallel to the boundaries are the same in all the layers. If the coordinate system is chosen such that z-axis is perpendicular to the boundaries and x–z is the plane of incidence, then the y-components of the wavevectors vanish and the x-components are given as $k_x=k_{0}sin\theta$, where $k_0=2\pi /\lambda$ is the size of the wavevector in vacuum and $\theta$ is the incidence angle in vacuum. The z-components of the wavevectors of the ordinary and extraordinary waves inside kth medium are equal to $\pm {\widehat{k}}_{\mathrm{o}z,k}$ and $\pm {\widehat{k}}_{\mathrm{e}z,k}$, where

$${\widehat{k}}_{\mathrm{o}z,k}=k_0\sqrt{{\widehat{\epsilon }}_{\mathrm{o},k}-{sin}^2\theta },{\widehat{k}}_{\mathrm{e}z,k}=k_0\sqrt{\frac{{\widehat{\epsilon }}_{\mathrm{o},k}}{{\widehat{\epsilon }}_{\mathrm{e},k}}\left({\widehat{\epsilon }}_{\mathrm{e},k}-{sin}^2\theta \right)}.$$

(3)

In the Yeh matrix formalism, the layered system is described by means of the boundary matrices and phase matrices. The boundary matrices take into account reflection and refraction of the waves at boundaries between adjacent media while the phase matrices describe how the phases of the waves change when they propagate through layers. The boundary separating the $(k-1)$th and kth media is described using the boundary matrix

$${\widehat{\mathbf{B}}}_{k-1,k}=\left(\begin{array}{cccc}{\widehat{B}}_{k-1,k}^{\left(1\right)}& {\widehat{B}}_{k-1,k}^{\left(2\right)}& 0& 0\\ {\widehat{B}}_{k-1,k}^{\left(2\right)}& {\widehat{B}}_{k-1,k}^{\left(1\right)}& 0& 0\\ 0& 0& {\widehat{B}}_{k-1,k}^{\left(3\right)}& {\widehat{B}}_{k-1,k}^{\left(4\right)}\\ 0& 0& {\widehat{B}}_{k-1,k}^{\left(4\right)}& {\widehat{B}}_{k-1,k}^{\left(3\right)}\end{array}\right),$$

(4)

with the elements of the matrix given as

$${\widehat{B}}_{k-1,k}^{\left(1\right)}=\frac{1}{2}\frac{{\widehat{\alpha }}_k}{{\widehat{\alpha }}_{k-1}}\left(\frac{{\widehat{Y}}_{\mathrm{e},k}}{{\widehat{Y}}_{\mathrm{e},k-1}}+1\right),{\widehat{B}}_{k-1,k}^{\left(2\right)}=-\frac{1}{2}\frac{{\widehat{\alpha }}_k}{{\widehat{\alpha }}_{k-1}}\left(\frac{{\widehat{Y}}_{\mathrm{e},k}}{{\widehat{Y}}_{\mathrm{e},k-1}}-1\right),$$

(5)

$${\widehat{B}}_{k-1,k}^{\left(3\right)}=\frac{1}{2}\left(\frac{{\widehat{Y}}_{\mathrm{o},k}}{{\widehat{Y}}_{\mathrm{o},k-1}}+1\right),{\widehat{B}}_{k-1,k}^{\left(4\right)}=\frac{1}{2}\left(\frac{{\widehat{Y}}_{\mathrm{o},k}}{{\widehat{Y}}_{\mathrm{o},k-1}}-1\right).$$

(6)

The optical admittances for the ordinary and extraordinary waves and the quantity ${\widehat{\alpha }}_k$ are given as

$${\widehat{Y}}_{\mathrm{o},k}=\frac{{\widehat{k}}_{\mathrm{o}z,k}}{k_0},{\widehat{Y}}_{\mathrm{e},k}=\frac{k_0{\widehat{\epsilon }}_{\mathrm{o},k}}{{\widehat{k}}_{\mathrm{e}z,k}},{\widehat{\alpha }}_k=\frac{{\widehat{k}}_{\mathrm{e}z,k}}{\sqrt{k_0^2{\widehat{\epsilon }}_{\mathrm{o},k}^2/{\widehat{\epsilon }}_{\mathrm{e},k}-\left({\widehat{\epsilon }}_{\mathrm{o},k}/{\widehat{\epsilon }}_{\mathrm{e},k}-1\right){\widehat{k}}_{\mathrm{e}z,k}^2}}.$$

(7)

For isotropic media, ${\widehat{\epsilon }}_{\mathrm{o},k}={\widehat{\epsilon }}_{\mathrm{e},k}={\widehat{\epsilon }}_k$ holds and the expressions for the z-components of the wavevectors (Equation (3)), optical admittances and quantity ${\widehat{\alpha }}_k$ (Equation (7)) can be written as

$$\begin{array}{l}{\widehat{k}}_{\mathrm{o}z,k}={\widehat{k}}_{\mathrm{e}z,k}={\widehat{k}}_{z,k}=k_0\sqrt{{\widehat{\epsilon }}_k-{sin}^2\theta }=k_0{\widehat{n}}_{k}cos{\widehat{\theta }}_k,{\widehat{\alpha }}_k=\frac{{\widehat{k}}_{z,k}}{k_0\sqrt{{\widehat{\epsilon }}_k}}=cos{\widehat{\theta }}_k,\\ {\widehat{Y}}_{\mathrm{o},k}={\widehat{Y}}_{\mathrm{s},k}=\frac{{\widehat{k}}_{z,k}}{k_0}={\widehat{n}}_{k}cos{\widehat{\theta }}_k,{\widehat{Y}}_{\mathrm{e},k}={\widehat{Y}}_{\mathrm{p},k}=\frac{k_0{\widehat{\epsilon }}_k}{{\widehat{k}}_{z,k}}=\frac{{\widehat{n}}_k}{cos{\widehat{\theta }}_k},\end{array}$$

(8)

where ${\widehat{n}}_k=\sqrt{{\widehat{\epsilon }}_k}$ is the complex refractive index and ${\widehat{\theta }}_k$ is the refraction angle.

The propagation of the waves inside the layer with thickness h formed by kth medium is described using the phase matrix

$${\widehat{\mathbf{T}}}_k\left(h\right)=\left(\begin{array}{cccc}{\mathrm{e}}^{-\mathrm{i}h{\widehat{k}}_{\mathrm{e}z,k}}& 0& 0& 0\\ 0& {\mathrm{e}}^{\mathrm{i}h{\widehat{k}}_{\mathrm{e}z,k}}& 0& 0\\ 0& 0& {\mathrm{e}}^{-\mathrm{i}h{\widehat{k}}_{\mathrm{o}z,k}}& 0\\ 0& 0& 0& {\mathrm{e}}^{\mathrm{i}h{\widehat{k}}_{\mathrm{o}z,k}}\end{array}\right).$$

(9)

The transfer matrix for multilayer system depicted in Figure 1 approximating the inhomogeneous film is calculated as a product of the boundary and phase matrices

$${\widehat{\mathbf{S}}}^{\mathrm{I}}={\widehat{\mathbf{T}}}_0(h/2){\widehat{\mathbf{B}}}_{0,1}{\widehat{\mathbf{T}}}_1\left(h\right){\widehat{\mathbf{B}}}_{1,2}{\widehat{\mathbf{T}}}_2\left(h\right)\cdots {\widehat{\mathbf{B}}}_{N-2,N-1}{\widehat{\mathbf{T}}}_{N-1}\left(h\right){\widehat{\mathbf{B}}}_{N-1,N}{\widehat{\mathbf{T}}}_N(h/2){\widehat{\mathbf{B}}}_{N,\mathrm{S}}.$$

(10)

The above matrix does not take into account the reflection and refraction of light on the boundary between the ambient and the upper side of the film (0th film). This is because this boundary exhibits roughness and, therefore, it cannot be represented by a boundary matrix (Equation (4)) which describes smooth boundaries. The formulae which include also reflection and refraction on the boundary between the ambient and the film are presented in the next section. The matrix ${\widehat{\mathbf{B}}}_{N,\mathrm{S}}$ represents the boundary between the lower side of the film (Nth film) and the substrate. Since the substrate is formed by isotropic medium, Equation (8) can be used to calculate the optical admittances in this case. The Fresnel reflection and transmission coefficients for the p-polarization (extraordinary waves) are given as

$${\widehat{r}}_{\mathrm{p}}^{\mathrm{I}}=\frac{{\widehat{S}}_{21}^{\mathrm{I}}}{{\widehat{S}}_{11}^{\mathrm{I}}},{\widehat{t}}_{\mathrm{p}}^{\mathrm{I}}=\frac{1}{{\widehat{S}}_{11}^{\mathrm{I}}},{\widehat{r}}_{\mathrm{p}}^{\prime \mathrm{I}}=-\frac{{\widehat{S}}_{12}^{\mathrm{I}}}{{\widehat{S}}_{11}^{\mathrm{I}}},{\widehat{t}}_{\mathrm{p}}^{\prime \mathrm{I}}=\frac{{\widehat{S}}_{11}^{\mathrm{I}}{\widehat{S}}_{22}^{\mathrm{I}}-{\widehat{S}}_{12}^{\mathrm{I}}{\widehat{S}}_{21}^{\mathrm{I}}}{{\widehat{S}}_{11}^{\mathrm{I}}},$$

(11)

where $S_{ij}^{\mathrm{I}}$ denote the components of the matrix in Equation (10) and for the s-polarization (ordinary waves) they are given as

$${\widehat{r}}_{\mathrm{s}}^{\mathrm{I}}=\frac{{\widehat{S}}_{43}^{\mathrm{I}}}{{\widehat{S}}_{33}^{\mathrm{I}}},{\widehat{t}}_{\mathrm{s}}^{\mathrm{I}}=\frac{1}{{\widehat{S}}_{33}^{\mathrm{I}}},{\widehat{r}}_{\mathrm{s}}^{\prime \mathrm{I}}=-\frac{{\widehat{S}}_{34}^{\mathrm{I}}}{{\widehat{S}}_{33}^{\mathrm{I}}},{\widehat{t}}_{\mathrm{s}}^{\prime \mathrm{I}}=\frac{{\widehat{S}}_{33}^{\mathrm{I}}{\widehat{S}}_{44}^{\mathrm{I}}-{\widehat{S}}_{43}^{\mathrm{I}}{\widehat{S}}_{43}^{\mathrm{I}}}{{\widehat{S}}_{33}^{\mathrm{I}}}.$$

(12)

The Fresnel coefficients without prime correspond to the incident wave falling onto the layered system from the top (ambient side) while the Fresnel coefficients with prime correspond to the incident wave falling from the bottom (substrate side). The formulae presented above are written in the convention that the reflection coefficients satisfy ${\widehat{r}}_{\mathrm{p}}=-{\widehat{r}}_{\mathrm{s}}$ for the normal incidence.

To obtain exact results corresponding to the inhomogeneous film, we should perform the limit in which the number of approximating homogeneous layers goes to infinity. In [27], it was shown that the convergence to the exact results can be improved using the Richardson extrapolation. The paper dealt with isotropic inhomogeneous layers but the results can be easily generalized to anisotropic inhomogeneous layers. If two steps of the Richardson extrapolation are used, the final formulae for the reflection coefficients can be written as

$${\widehat{r}}_q^{\mathrm{I}}=\frac{16}{15}\left(\frac{4}{3}{\widehat{r}}_q\left(N\right)-\frac{1}{3}{\widehat{r}}_q(N/2)\right)-\frac{1}{15}\left(\frac{4}{3}{\widehat{r}}_q(N/2)-\frac{1}{3}{\widehat{r}}_q(N/4)\right),$$

(13)

where the index $q=\mathrm{p},\mathrm{s}$ is used to distinguish the polarizations and ${\widehat{r}}_q\left(N\right)$, ${\widehat{r}}_q(N/2)$ and ${\widehat{r}}_q(N/4)$ are the reflection coefficients calculated using Equations (10)–(12) with the number of dividing boundaries set to N, $N/2$ and $N/4$, respectively. The formulae for the other Fresnel coefficients ${\widehat{t}}_q^{\mathrm{I}}$, ${\widehat{r}}_q^{\prime \mathrm{I}}$, and ${\widehat{t}}_q^{\prime \mathrm{I}}$ are analogous. The number N (which must be divisible by four) determines the precision of the calculated Fresnel coefficients of the inhomogeneous film and should be fixed in sufficiently large value.

2.3. Roughness of the Upper Boundaries

The influence of random roughness on the optical quantities of the SiN${}_x$ films is described using the scalar diffraction theory (SDT). This theory was utilized for deriving the formulae for the optical quantities of homogeneous and isotropic thin films with randomly rough boundaries (see, e.g., [18,19,20,21,22]). In this paper, a generalization of this approach of the SDT is performed for inhomogeneous thin films exhibiting random roughness of the upper boundaries and uniaxial anisotropy with the optical axis perpendicular to the mean planes of the boundaries. From formulae expressing the optical quantities of the homogeneous isotropic films with randomly rough boundaries presented in the cited papers [18,19,20,21,22], it can be implied that the mean values of the reflection and transmission coefficients of the studied SiN${}_x$ films are given as follows:

$$\begin{array}{l}{\widehat{r}}_q^{}={\int }_{-\infty }^{\infty }w\left(z\right){\widehat{r}}_q^{\mathrm{Q}}\left(z\right){\mathrm{e}}^{-2\mathrm{i}z{\widehat{k}}_{z,\mathrm{A}}}\mathrm{d}z,{\widehat{t}}_q^{}={\int }_{-\infty }^{\infty }w\left(z\right){\widehat{t}}_q^{\mathrm{Q}}\left(z\right){\mathrm{e}}^{-\mathrm{i}z{\widehat{k}}_{z,\mathrm{A}}}\mathrm{d}z,\\ {\widehat{r}}_q^{\prime }={\int }_{-\infty }^{\infty }w\left(z\right){\widehat{r}}_q^{\prime \mathrm{Q}}\left(z\right)\mathrm{d}z,{\widehat{t}}_q^{\prime }={\int }_{-\infty }^{\infty }w\left(z\right){\widehat{t}}_q^{\prime \mathrm{Q}}\left(z\right){\mathrm{e}}^{\mathrm{i}z{\widehat{k}}_{z,\mathrm{A}}}\mathrm{d}z,\end{array}$$

(14)

where $z=\eta (x,y)$ denotes a random function describing heights of roughness related to the mean plane of the upper boundary, $w\left(z\right)$ is the one-dimensional probability density function for the distribution of the heights and ${\widehat{k}}_{z,\mathrm{A}}$ is the z-component of the wavevector in the ambient. The symbols ${\widehat{r}}_q^{\mathrm{Q}}\left(z\right)$, ${\widehat{t}}_q^{\mathrm{Q}}\left(z\right)$, ${\widehat{r}}_q^{\prime \mathrm{Q}}\left(z\right)$, and ${\widehat{t}}_q^{\prime \mathrm{Q}}\left(z\right)$ denote the local reflection and transmission coefficients on the boundaries of the SiN${}_x$ films, which are functions of the heights. The coefficients ${\widehat{r}}_q^{}$ and ${\widehat{t}}_q^{}$ correspond to the incidence of light from the top and coefficients ${\widehat{r}}_q^{\prime }$ and ${\widehat{t}}_q^{\prime }$ correspond to the incidence of light from the bottom. It should be noted that the formulae in Equation (14) are valid if the slopes of random roughness of the upper boundaries are sufficiently small so that they can be neglected. The local reflection and transmission coefficients are expressed in this way:

$$\begin{array}{l}{\widehat{r}}_q^{\mathrm{Q}}\left(z\right)={\widehat{r}}_q^{\mathrm{U}}+{\widehat{t}}_q^{\mathrm{U}}{\widehat{r}}_q^{\mathrm{I}}{\widehat{t}}_q^{\prime \mathrm{U}}{\mathrm{e}}^{2\mathrm{i}z{\widehat{k}}_{qz,0}}{\displaystyle \sum _{m=0}^{\infty }}{\left({\widehat{r}}_q^{\mathrm{I}}{\widehat{r}}_q^{\prime \mathrm{U}}{\mathrm{e}}^{2\mathrm{i}z{\widehat{k}}_{qz,0}}\right)}^m,\\ {\widehat{t}}_q^{\mathrm{Q}}\left(z\right)={\widehat{t}}_q^{\mathrm{U}}{\widehat{t}}_q^{\mathrm{I}}{\mathrm{e}}^{\mathrm{i}z{\widehat{k}}_{qz,0}}{\displaystyle \sum _{m=0}^{\infty }}{\left({\widehat{r}}_q^{\mathrm{I}}{\widehat{r}}_q^{\prime \mathrm{U}}{\mathrm{e}}^{2\mathrm{i}z{\widehat{k}}_{qz,0}}\right)}^m,\\ {\widehat{r}}_q^{\prime \mathrm{Q}}\left(z\right)={\widehat{r}}_q^{\prime \mathrm{I}}+{\widehat{t}}_q^{\prime \mathrm{I}}{\widehat{r}}_q^{\prime \mathrm{U}}{\widehat{t}}_q^{\prime \mathrm{U}}{\mathrm{e}}^{\mathrm{i}z{\widehat{k}}_{qz,0}}{\displaystyle \sum _{m=0}^{\infty }}{\left({\widehat{r}}_q^{\mathrm{I}}{\widehat{r}}_q^{\prime \mathrm{U}}{\mathrm{e}}^{2\mathrm{i}z{\widehat{k}}_{qz,0}}\right)}^m,\\ {\widehat{t}}_q^{\prime \mathrm{Q}}\left(z\right)={\widehat{t}}_q^{\prime \mathrm{I}}{\widehat{t}}_q^{\prime \mathrm{U}}{\mathrm{e}}^{\mathrm{i}z{\widehat{k}}_{qz,0}}{\displaystyle \sum _{m=0}^{\infty }}{\left({\widehat{r}}_q^{\mathrm{I}}{\widehat{r}}_q^{\prime \mathrm{U}}{\mathrm{e}}^{2\mathrm{i}z{\widehat{k}}_{qz,0}}\right)}^m,\end{array}$$

(15)

where ${\widehat{r}}_q^{\mathrm{U}}$, ${\widehat{t}}_q^{\mathrm{U}}$, ${\widehat{r}}_q^{\prime \mathrm{U}}$ and ${\widehat{t}}_q^{\prime \mathrm{U}}$ are the Fresnel reflection and transmission coefficients of the upper boundary. They can be expressed as

$$\begin{array}{l}{\widehat{r}}_{\mathrm{p}}^{\mathrm{U}}=-\frac{{\widehat{Y}}_{\mathrm{e},0}-{\widehat{Y}}_{\mathrm{p},A}}{{\widehat{Y}}_{\mathrm{e},0}+{\widehat{Y}}_{\mathrm{p},A}},{\widehat{t}}_{\mathrm{p}}^{\mathrm{U}}=\frac{{\widehat{\alpha }}_{\mathrm{A}}}{{\widehat{\alpha }}_0}\frac{2{\widehat{Y}}_{\mathrm{p},A}}{{\widehat{Y}}_{\mathrm{e},0}+{\widehat{Y}}_{\mathrm{p},A}},{\widehat{r}}_{\mathrm{p}}^{\prime \mathrm{U}}=-\frac{{\widehat{Y}}_{\mathrm{p},A}-{\widehat{Y}}_{\mathrm{e},0}}{{\widehat{Y}}_{\mathrm{p},A}+{\widehat{Y}}_{\mathrm{e},0}},{\widehat{t}}_{\mathrm{p}}^{\prime \mathrm{U}}=\frac{{\widehat{\alpha }}_0}{{\widehat{\alpha }}_{\mathrm{A}}}\frac{2{\widehat{Y}}_{\mathrm{e},0}}{{\widehat{Y}}_{\mathrm{p},A}+{\widehat{Y}}_{\mathrm{e},0}},\\ {\widehat{r}}_{\mathrm{s}}^{\mathrm{U}}=\frac{{\widehat{Y}}_{\mathrm{o},0}-{\widehat{Y}}_{\mathrm{s},A}}{{\widehat{Y}}_{\mathrm{o},0}+{\widehat{Y}}_{\mathrm{s},A}},{\widehat{t}}_{\mathrm{s}}^{\mathrm{U}}=\frac{{\widehat{\alpha }}_{\mathrm{A}}}{{\widehat{\alpha }}_0}\frac{2{\widehat{Y}}_{\mathrm{s},A}}{{\widehat{Y}}_{\mathrm{o},0}+{\widehat{Y}}_{\mathrm{s},A}},{\widehat{r}}_{\mathrm{s}}^{\prime \mathrm{U}}=\frac{{\widehat{Y}}_{\mathrm{s},A}-{\widehat{Y}}_{\mathrm{o},0}}{{\widehat{Y}}_{\mathrm{s},A}+{\widehat{Y}}_{\mathrm{o},0}},{\widehat{t}}_{\mathrm{s}}^{\prime \mathrm{U}}=\frac{{\widehat{\alpha }}_0}{{\widehat{\alpha }}_{\mathrm{A}}}\frac{2{\widehat{Y}}_{\mathrm{o},0}}{{\widehat{Y}}_{\mathrm{s},A}+{\widehat{Y}}_{\mathrm{o},0}},\end{array}$$

(16)

where the optical admittances ${\widehat{Y}}_{\mathrm{e},0}$, ${\widehat{Y}}_{\mathrm{o},0}$ and the quantity ${\widehat{\alpha }}_0$ at the upper side of the inhomogeneous film are calculated using Equation (7) and the optical admittances ${\widehat{Y}}_{\mathrm{p},A}$, ${\widehat{Y}}_{\mathrm{s},A}$ and the quantity ${\widehat{\alpha }}_{\mathrm{A}}$ are calculated using Equation (8) for the ambient. The symbols ${\widehat{r}}_q^{\mathrm{I}}$, ${\widehat{t}}_q^{\mathrm{I}}$ and ${\widehat{r}}_q^{\prime \mathrm{I}}$, ${\widehat{t}}_q^{\prime \mathrm{I}}$ denote the effective reflection and transmission coefficients representing the inhomogeneous thin film calculated using Equations (11) and (12) (direction of incidence light is distinguished in the same way as for the coefficients of the upper boundary). It is assumed that the SiN${}_x$ films exhibit the constant refractive index and extinction coefficient within the range of random roughness of the upper boundaries. After inserting Equation (15) into Equation (14), one obtains the following expressions for the reflection and transmission coefficients of the characterized SiN${}_x$:

$$\begin{array}{l}{\widehat{r}}_q^{}={\widehat{r}}_q^{\mathrm{U}}\widehat{\chi }(-2{\widehat{k}}_z^{\mathrm{A}})+{\widehat{t}}_q^{\mathrm{U}}{\widehat{t}}_q^{\prime \mathrm{U}}{\displaystyle \sum _{m=0}^{\infty }}{\left({\widehat{r}}_q^{\mathrm{I}}\right)}^{m+1}{\left({\widehat{r}}_q^{\prime \mathrm{U}}\right)}^m\widehat{\chi }\left({\widehat{B}}_m\right),{\widehat{t}}_q^{}={\widehat{t}}_q^{\mathrm{U}}{\widehat{t}}_q^{\mathrm{I}}{\displaystyle \sum _{m=0}^{\infty }}{\left({\widehat{r}}_q^{\mathrm{I}}\right)}^m{\left({\widehat{r}}_q^{\prime \mathrm{U}}\right)}^m\widehat{\chi }\left({\widehat{M}}_m\right),\\ {\widehat{r}}_q^{\prime }={\widehat{r}}_q^{\prime \mathrm{I}}+{\widehat{t}}_q^{\prime \mathrm{I}}{\widehat{t}}_q^{\mathrm{I}}{\displaystyle \sum _{m=0}^{\infty }}{\left({\widehat{r}}_q^{\mathrm{I}}\right)}^m{\left({\widehat{r}}_q^{\prime \mathrm{U}}\right)}^{m+1}\widehat{\chi }\left({\widehat{C}}_m\right),{\widehat{t}}_q^{\prime }={\widehat{t}}_q^{\prime \mathrm{I}}{\widehat{t}}_q^{\prime \mathrm{U}}{\displaystyle \sum _{m=0}^{\infty }}{\left({\widehat{r}}_q^{\mathrm{I}}\right)}^m{\left({\widehat{r}}_q^{\prime \mathrm{U}}\right)}^m\widehat{\chi }\left({\widehat{M}}_m\right).\end{array}$$

(17)

where symbol $\widehat{\chi }\left(v\right)$ denotes the characteristic function. For the Gaussian distribution, it is given as

$$\widehat{\chi }\left(v\right)=exp\left(-\frac{1}{2}{\sigma }^2{v}^2\right),$$

(18)

where $\sigma$ is the rms value of the heights of the roughness and v represents one of the following quantities

$${\widehat{B}}_m=(2m+2){\widehat{k}}_{q,z}^{\mathrm{U}}-2{\widehat{k}}_z^{\mathrm{A}},{\widehat{M}}_m=(2m+1){\widehat{k}}_{q,z}^{\mathrm{U}}-{\widehat{k}}_z^{\mathrm{A}},{\widehat{C}}_m=(2m+2){\widehat{k}}_{q,z}^{\mathrm{U}}.$$

(19)

2.4. Reflection on the Back Sides of the Substrates

The silicon substrates become transparent in the infrared region, therefore the depolarization caused by the reflections on the back sides of these substrates must be taken into account. The incoherent formalism working with the Mueller matrices was used to include this effect into the formulae for the measured optical quantities.

The system can be separated into three parts, the layered system on the front side of the substrate (i.e., the silicon nitride film), the silicon substrate, and the layered system on the back side (i.e., the native oxide layer). For the layered system on the front side of the substrate, the Mueller matrix representing reflection can be calculated as

$${\mathbf{R}}_{\mathrm{t}}=\frac{1}{2}\left(\begin{array}{cccc}|{\widehat{r}}_{\mathrm{p}}{|}^2+{\left|{\widehat{r}}_{\mathrm{s}}\right|}^2& |{\widehat{r}}_{\mathrm{p}}{|}^2-{\left|{\widehat{r}}_{\mathrm{s}}\right|}^2& 0& 0\\ |{\widehat{r}}_{\mathrm{p}}{|}^2-{\left|{\widehat{r}}_{\mathrm{s}}\right|}^2& |{\widehat{r}}_{\mathrm{p}}{|}^2+{\left|{\widehat{r}}_{\mathrm{s}}\right|}^2& 0& 0\\ 0& 0& {\widehat{r}}_{\mathrm{p}}^{*}{\widehat{r}}_{\mathrm{s}}+{\widehat{r}}_{\mathrm{p}}{\widehat{r}}_{\mathrm{s}}^{*}& \mathrm{i}{\widehat{r}}_{\mathrm{p}}^{*}{\widehat{r}}_{\mathrm{s}}-\mathrm{i}{\widehat{r}}_{\mathrm{p}}{\widehat{r}}_{\mathrm{s}}^{*}\\ 0& 0& -\mathrm{i}{\widehat{r}}_{\mathrm{p}}^{*}{\widehat{r}}_{\mathrm{s}}+\mathrm{i}{\widehat{r}}_{\mathrm{p}}{\widehat{r}}_{\mathrm{s}}^{*}& {\widehat{r}}_{\mathrm{p}}^{*}{\widehat{r}}_{\mathrm{s}}+{\widehat{r}}_{\mathrm{p}}{\widehat{r}}_{\mathrm{s}}^{*}\end{array}\right).$$

(20)

This matrix corresponds to the situation in which the wave is incident from the top. Analogous formulae can be used to calculate the Mueller matrix ${\mathbf{T}}_{\mathrm{t}}$ corresponding to transmission of light and Mueller matrices ${\mathbf{R}}_{\mathrm{t}}^{\prime }$, ${\mathbf{T}}_{\mathrm{t}}^{\prime }$ corresponding to reflection and transmission of the wave incident from the bottom. The propagation of the waves through the isotropic substrate is described using the Mueller matrices

$$\mathbf{U}={\mathbf{U}}^{\prime }=\mathbf{1}exp(-2d_{\mathrm{S}}\Im \left({\widehat{k}}_z^{\mathrm{S}}\right)),$$

where $\Im \left({\widehat{k}}_z^{\mathrm{S}}\right)$ is the imaginary part of the wavevector in the substrate and $d_{\mathrm{S}}$ is the thickness of the substrate. The Mueller matrices ${\mathbf{R}}_{\mathrm{b}}$, ${\mathbf{T}}_{\mathrm{b}}$, ${\mathbf{R}}_{\mathrm{b}}^{\prime }$, and ${\mathbf{T}}_{\mathrm{b}}^{\prime }$ corresponding to the reflection and transmission of waves on the layered system on the back side of the substrate are calculated in the same way as the Mueller matrices corresponding to the front side.

The Mueller matrix for the reflection of the sample is then calculated as

$$\begin{array}{l}\mathbf{R}={\mathbf{R}}_{\mathrm{t}}+c_1{\mathbf{T}}_{\mathrm{t}}^{\prime }{\mathbf{U}}^{\prime }{\mathbf{R}}_{\mathrm{b}}\mathbf{U}{\mathbf{T}}_{\mathrm{t}}+c_2{\mathbf{T}}_{\mathrm{t}}^{\prime }{\mathbf{U}}^{\prime }{\mathbf{R}}_{\mathrm{b}}\mathbf{U}{\mathbf{R}}_{\mathrm{t}}^{\prime }{\mathbf{U}}^{\prime }{\mathbf{R}}_{\mathrm{b}}\mathbf{U}{\mathbf{T}}_{\mathrm{t}}+c_3{\mathbf{T}}_{\mathrm{t}}^{\prime }{\mathbf{U}}^{\prime }{\mathbf{R}}_{\mathrm{b}}\mathbf{U}{\mathbf{R}}_{\mathrm{t}}^{\prime }{\mathbf{U}}^{\prime }{\mathbf{R}}_{\mathrm{b}}\mathbf{U}{\mathbf{R}}_{\mathrm{t}}^{\prime }{\mathbf{U}}^{\prime }{\mathbf{R}}_{\mathrm{b}}\mathbf{U}{\mathbf{T}}_{\mathrm{t}}+\cdots \\ ={\mathbf{R}}_{\mathrm{t}}+{\displaystyle \sum _{n=0}^{\infty }}{c}_{n+1}{\mathbf{T}}_{\mathrm{t}}^{\prime }{\mathbf{U}}^{\prime }{\mathbf{R}}_{\mathrm{b}}\mathbf{U}{\left({\mathbf{R}}_{\mathrm{t}}^{\prime }{\mathbf{U}}^{\prime }{\mathbf{R}}_{\mathrm{b}}\mathbf{U}\right)}^n{\mathbf{T}}_{\mathrm{t}},\end{array}$$

(21)

where the first term corresponds to the reflection on the front side of the substrate, the second term corresponds to the path of the beam reflected once from the back side and subsequent terms correspond to the paths involving multiple reflections on the back side. Since the substrate has finite thickness, the outgoing beams corresponding to paths involving reflections from the back side are shifted with respect to the beam reflected on the front side. This means that only part of these beams falls onto the detector. The factors $c_n$ are introduced to account for this effect. If the beam has circular cross section, then the measured areas on the sample surface are ellipses. The factor $c_n$ representing the fraction of the beam falling onto the detector is then proportional to the area in which the ellipse representing the area measured by the detector intersects with the ellipse representing the reflected beam.

The Mueller matrix can be written using the normalized Mueller matrix as

$$\mathbf{R}=R\left(\begin{array}{cccc}1& -I_{\mathrm{n}}& 0& 0\\ -I_{\mathrm{n}}& 1& 0& 0\\ 0& 0& I_{\mathrm{c}}& I_{\mathrm{s}}\\ 0& 0& -I_{\mathrm{s}}& I_{\mathrm{c}}\end{array}\right)$$

(22)

where $R=M_{00}$ is reflectance and $I_{\mathrm{s}}$, $I_{\mathrm{c}}$, and $I_{\mathrm{n}}$ are the quantities measured in ellipsometry. If there is no depolarization, then they correspond to the associated ellipsometric parameters.

2.5. Dispersion Model

The optical constants of crystalline float-zone silicon determined in [28] were used for the substrate. The optical constants determined for amphorous SiO${}_2$ in [29] were used for the native oxide layers on the back sides of the substrates.

The optical constants of the inhomogeneous films were described using the universal dispersion model [30,31]. The dielectric tensor obtained by this model satisfies three fundamental conditions that must be fulfilled by any physically correct model. The first condition is the time-reversal symmetry

$${\widehat{\mathit{\epsilon }}}^{*}\left(E\right)=\widehat{\mathit{\epsilon }}(-E),$$

(23)

where E is the photon energy and ∗ denotes the complex conjugation. The second condition is the Kramers–Kronig consistency

$${\mathit{\epsilon }}_{\mathrm{r}}\left(E\right)=\mathbf{1}+\frac{1}{\pi }-{\int }_{-\infty }^{\infty }\frac{{\mathit{\epsilon }}_{\mathrm{i}}\left(E^{\prime }\right)}{E^{\prime }-E}\mathrm{d}{E}^{\prime }$$

(24)

where ${\mathit{\epsilon }}_{\mathrm{r}}\left(E\right)$ and ${\mathit{\epsilon }}_{\mathrm{i}}\left(E\right)$ denote the real and the imaginary parts of the dielectric tensor and $\mathbf{1}$ is the unit matrix. The dash over the integral sign is used to emphasize that it must be evaluated as the Cauchy principal value. The last condition is the Thomas–Reiche–Kuhn sum rule

$${\int }_0^{\infty }E{\mathit{\epsilon }}_{\mathrm{i}}\left(E\right)\mathrm{d}E=N\mathbf{1},$$

(25)

where N is the total transition strength related to the density of atoms ${\mathcal{N}}_{\mathrm{a}}$ and the mean number of electrons per atom Z as [32]

$$N=\frac{{\left(eh\right)}^2}{8\pi {\epsilon }_0{m}_{\mathrm{e}}}Z{\mathcal{N}}_{\mathrm{a}}\mathcal{U},$$

(26)

where e is the electron charge, h is the Planck constant, ${\epsilon }_0$ is vacuum permittivity, $m_{\mathrm{e}}$ is electron mass and $\mathcal{U}\approx 1$ is the constant specific to individual materials.

The conditions presented above are expressed using the dielectric tensor. If the form of the dielectric tensor (Equation (1)) describing the uniaxial anisotropy is taken into account, then these conditions can be written as

$${\widehat{\epsilon }}_q^{*}\left(E\right)={\widehat{\epsilon }}_q(-E),{\epsilon }_{\mathrm{r},q}\left(E\right)=1+\frac{1}{\pi }-{\int }_{-\infty }^{\infty }\frac{{\epsilon }_{\mathrm{i},q}\left(E^{\prime }\right)}{E^{\prime }-E}\mathrm{d}{E}^{\prime },{\int }_0^{\infty }E{\epsilon }_{\mathrm{i},q}\left(E\right)\mathrm{d}E=N,$$

(27)

where the index q is either $\mathrm{o}$ (ordinary) or $\mathrm{e}$ (extraordinary). These conditions show that the ordinary and extraordinary dielectric functions could be modeled independently on each other, the only restriction is that they should give the same value of the sum rule integral.

In the universal dispersion model, the dielectric function is expressed as a sum of contributions, which model various elementary excitations contributing to the dielectric response, e.g., direct and indirect excitations of electrons, phonon absorption, excitations involving localized states, core electron excitations, etc. The number and types of contributions must be selected such that they reflect the type of the material and the spectral range for which the model is used.

Since the same model is used for the ordinary and extraordinary dielectric functions and repeating the index $q=\mathrm{o},\mathrm{e}$ in all places would make the formulae unwieldy, we omit this index. This means that the formulae are presented in the form that would be usable if the medium were isotropic. The index q is reintroduced when the need to distinguish between the ordinary and extraordinary dielectric functions arises.

A basic formula used to express the dielectric function in the universal dispersion model is

$$\widehat{\epsilon }\left(E\right)=1+\sum _p{N}_p{\widehat{\epsilon }}_p^0\left(E\right),$$

(28)

where the index p distinguishes individual contributions, the parameters $N_p$, called the transition strengths, determine the strengths of the contributions and the functions ${\widehat{\epsilon }}_p^0\left(E\right)$, which are in the nomenclature of the universal dispersion model called the normalized contributions to the dielectric function, determine the spectral distribution of the individual contributions. These functions are normalized to unity with respect to the sum rule integral, i.e. their imaginary parts fulfill the condition

$${\int }_0^{\infty }E{\epsilon }_{\mathrm{i},p}^0\left(E\right)\mathrm{d}E=1.$$

(29)

The imaginary parts of the normalized contributions are modeled on the basis of theoretical and empirical knowledge, the real parts are then calculated using the Kramers–Kronig relation as

$${\epsilon }_{\mathrm{r},p}^0\left(E\right)=\frac{2}{\pi }-{\int }_0^{\infty }\frac{E^{\prime }{\epsilon }_{\mathrm{i},p}^0\left(E^{\prime }\right)}{E^{\prime 2}-E^2}\mathrm{d}{E}^{\prime }.$$

(30)

In the case of the model for the dielectric function of the inhomogeneous film, two contributions are considered. The dielectric function is expressed as

$$\widehat{\epsilon }\left(E\right)=1+N_{\mathrm{vc}}{\widehat{\epsilon }}_{\mathrm{dt}}^0\left(E\right)+N_{\mathrm{ut}}{\widehat{\epsilon }}_{\mathrm{ut}}^0\left(E\right),$$

(31)

where the first contribution represents the valence to conduction band transitions of electrons and the second contribution represents the Urbach tail, which describes weak absorption below the band gap energy.

The imaginary part of the normalized contribution to the dielectric function for the valence to conduction band transitions of electrons is expressed as

$${\epsilon }_{\mathrm{i},\mathrm{dt}}^0\left(E\right)=\frac{1}{1+A_{\mathrm{ex}}}{\epsilon }_{\mathrm{i},\mathrm{vc}}^0\left(E\right)+\frac{A_{\mathrm{ex}}}{1+A_{\mathrm{ex}}}{\epsilon }_{\mathrm{i},\mathrm{ex}}^0\left(E\right)$$

(32)

where the first term with [30]

$${\epsilon }_{\mathrm{i},\mathrm{vc}}^0\left(E\right)=\frac{{(E-E_{\mathrm{g}})}^2{(E_{\mathrm{h}}-E)}^2}{C_{\mathrm{vc}}{E}^2}\mathsf{\Theta }(E;E_{\mathrm{g}},E_{\mathrm{h}})$$

(33)

models the basic absorption band with parameters $E_{\mathrm{g}}$ and $E_{\mathrm{h}}$ being the minimal and the maximal energy of transitions. The function $\mathsf{\Theta }(E;E_{\mathrm{g}},E_{\mathrm{h}})$ is defined as unity for $E_{\mathrm{g}}\le E\le E_{\mathrm{h}}$ and it is equal to zero for other photon energies. The second term with [33]

$${\epsilon }_{\mathrm{i},\mathrm{ex}}^0\left(E\right)=\frac{{(E-E_{\mathrm{g}})}^2{(E_{\mathrm{h}}-E)}^2}{C_{\mathrm{ex}}{E}^2[{(E-E_{\mathrm{ex}})}^2+B_{\mathrm{ex}}^2]}\mathsf{\Theta }(E;E_{\mathrm{g}},E_{\mathrm{h}})$$

(34)

introduces modification of the shape of the basic absorption band. It corresponds to the absorption band containing peak with central energy and width determined by parameters $E_{\mathrm{ex}}$ and $B_{\mathrm{ex}}$. The subscript $\mathrm{ex}$ is derived from the fact that it could be used to model the excitonic effects but in the most general context it should be interpreted just as a mechanism to modify the shape of the absorption band. The parameter $A_{\mathrm{ex}}$ is used to control the weights associated with ${\epsilon }_{\mathrm{i},\mathrm{vc}}^0\left(E\right)$ and ${\epsilon }_{\mathrm{i},\mathrm{ex}}^0\left(E\right)$.

The imaginary part of the normalized contribution to the dielectric function representing the Urbach tail is given as [33]

$$\begin{array}{l}{\epsilon }_{\mathrm{i},\mathrm{ut}}^0\left(E\right)=\frac{1}{C_{\mathrm{ut}}{E}^2}\left[\left(1+\frac{(E-E_{\mathrm{g}})(E_{\mathrm{m}}-E)}{E_{\mathrm{u}}(E_{\mathrm{m}}-E_{\mathrm{g}})}\right)\mathsf{\Theta }(E;E_{\mathrm{g}},E_{\mathrm{m}})\right.\\ \left.+{\mathrm{e}}^{\frac{E-E_{\mathrm{g}}}{E_{\mathrm{u}}}}\mathsf{\Theta }(E;E_{\mathrm{g}}/2,E_{\mathrm{g}})+{\mathrm{e}}^{\frac{E_{\mathrm{m}}-E}{E_{\mathrm{u}}}}\mathsf{\Theta }(E;E_{\mathrm{m}},E_{\mathrm{m}}+E_{\mathrm{g}}/2)-{\mathrm{e}}^{-\frac{E_{\mathrm{g}}}{2E_{\mathrm{u}}}}\right]\end{array}$$

(35)

where $E_{\mathrm{m}}=(E_{\mathrm{g}}+E_{\mathrm{h}})/2$ and the parameter $E_{\mathrm{u}}$, called the Urbach energy, determines how fast is the exponential tail decaying.

The constants $C_{\mathrm{vc}}$, $C_{\mathrm{ex}}$ and $C_{\mathrm{ut}}$ in Equations (33), (34) and (35) must be calculated such that the sum rule normalization condition in Equation (29) is satisfied.

The dielectric model described above is applicable to homogeneous isotropic films. The dielectric function given by Equation (31) is a function of the photon energy E and the dispersion parameters $N_{\mathrm{vc}},E_{\mathrm{g}},E_{\mathrm{h}},A_{\mathrm{ex}},E_{\mathrm{ex}},B_{\mathrm{ex}},N_{\mathrm{ut}},E_{\mathrm{u}}$, which for the sake of brevity are represented by a vector $p_{\alpha }$ with the index $\alpha$ distinguishing the individual dispersion parameters. Thus, if we want to indicate dependence on all the parameters, then the dielectric function should be written as $\widehat{\epsilon }(E,p_{\alpha })$.

To describe the dielectric response of media with uniaxial anisotropy, it is necessary to specify the ordinary and extraordinary dielectric functions. Both the dielectric functions are calculated using the model given in Equation (31) but each uses a separate set of dispersion parameters: the parameters $p_{\alpha ,\mathrm{o}}$ are used for the ordinary dielectric function and the parameters $p_{\alpha ,\mathrm{e}}$ are used for the extraordinary dielectric function. Therefore, the ordinary and extraordinary dielectric functions can be written as

$${\widehat{\epsilon }}_{\mathrm{o}}\left(E\right)=\widehat{\epsilon }(E,p_{\alpha ,\mathrm{o}}),{\widehat{\epsilon }}_{\mathrm{e}}\left(E\right)=\widehat{\epsilon }(E,p_{\alpha ,\mathrm{e}}).$$

(36)

We should emphasize that the transition strengths parameterizing the ordinary and extraordinary dielectric functions are not independent; they must be chosen such that the total transition strength (value of the sum rule integral) is the same in both cases. This means that they must fulfill the condition

$$N_{\mathrm{vc},\mathrm{e}}+N_{\mathrm{ut},\mathrm{e}}=N_{\mathrm{vc},\mathrm{o}}+N_{\mathrm{ut},\mathrm{o}}.$$

(37)

To describe the inhomogeneous films, it is also necessary to introduce the dependence on the coordinate z (depth inside film). The inhomogeneity of the studied films is modeled by introducing the dependence on z for the dispersion parameters, therefore the dielectric functions can be calculated as

$${\widehat{\epsilon }}_q(E,z)=\widehat{\epsilon }(E,p_{\alpha ,q}\left(z\right)),$$

(38)

where $q=\mathrm{e},\mathrm{o}$ and the functions $p_{\alpha ,q}\left(z\right)$ represent the profiles of the dispersion parameters. The inhomogeneity of the studied films was described using the profiles of the dispersion parameters $p_{\alpha ,q}\left(z\right)$ expressed as

$$p_{\alpha ,q}\left(z\right)={\mathrm{e}}^{-z/w}{p}_{\alpha ,q}^{\mathrm{U}}+(1-{\mathrm{e}}^{-z/w})p_{\alpha ,q}^{\infty },$$

(39)

where d is the film thickness, w is the parameter controlling the shape of the profile, $p_{\alpha ,q}^{\mathrm{U}}$ are the values at the upper boundary of the inhomogeneous film and $p_{\alpha ,q}^{\infty }$ are the values of the parameters for $z\to \infty$.

2.6. Imaging Spectroscopic Reflectometry

Imaging spectroscopic reflectometry (ISR) [34,35,36,37,38] is a technique in which the sample is illuminated by a monochromatic light and its image in reflected light is taken using the CCD camera. The signals registered by individual pixels of the CCD camera correspond to light reflected from small areas on the sample. By changing the wavelength of the incident light it is possible to obtain the spectral dependencies of the reflectance corresponding to small areas on the sample.

The imaging spectroscopic reflectometry is well suited for studying thickness non-uniformity of the samples. If the film is non-uniform in thickness, then different mean thicknesses are observed by individual pixels. However, since the areas on the sample corresponding to individual pixels are very small, the variations of the thickness within these small areas can be neglected. This means that the reflectance corresponding to the individual pixel can be calculated assuming the film is uniform in thickness within the corresponding region. By processing the experimental data obtained by the imaging spectroscopic reflectometry it is possible to reconstruct the maps of the local thicknesses. In this work, a simple single-pixel method is used to process the experimental data. In the single-pixel method, the spectral dependencies of reflectance are processed independently for each pixel.

3. Experiment

3.1. Sample Preparation

The thin films of non-stoichiometric silicon nitride were deposited onto one side of double-side polished silicon wafers. The deposition was performed by the reactive magnetron sputtering [39] of silicon target in argon-nitrogen atmosphere. The films were deposited using Vinci Technologies PVD 50S thin film deposition system. The pressure during the deposition was kept at constant value ${10}^{-2}$ mbar, and the flow rates of argon and nitrogen were $Q_{\mathrm{Ar}}=80$ sccm and $Q_{{\mathrm{N}}_2}=5$ sccm. The power supply for the magnetron was set at 100 W. The films were deposited onto cold (unheated) substrates. Three samples with the same deposition conditions but different deposition times 30, 45 and 90 min were prepared.

3.2. Experimental Data

The spectral dependencies of the ellipsometric parameters were measured using the Horiba Jobin Yvon UVISEL phase modulated ellipsometer (HORIBA Jobin Yvon, Irvine, CA, USA) for five incidence angles within the interval 55–75${}^{\circ }$ and within the spectral range 0.6–6.3 eV (197–2066 nm). The spectral dependencies of the reflectance were measured by the Perkin Elmer Lambda 1050 spectrophotometer (Perkin Elmer, Waltham, MA, USA) at near normal incidence angle 6${}^{\circ }$ within the spectral range 0.7–6.5 eV (190–1800 nm). The ISR data were measured at normal incidence within the spectral range 1.7–4.5 eV (275–700 nm), and the pixel size was $40\times 40\mu$m. The Bruker Dimension Icon atomic force microscope (Bruker, Billerica, MA, USA) was used to investigate topography of surfaces of the samples. The X-ray diffraction curves have been measured on the PaNalytical X’Pert MRD diffractometer (Malvern PaNalytical, Malvern, UK), using standard X-ray tube (Cu, long-fine focus, 45 kV, 30 mA).

3.3. Data Processing

To interpret the experimental data, it was necessary to assume that the non-stoichiometric silicon nitride films are inhomogeneous and that they also exhibit uniaxial anisotropy. The inhomogeneity of the films was modeled by considering the profiles of the dispersion parameters given in Equation (39). Two values must be provided for each dispersion parameter, one specifying its value at the upper boundary of the film (i.e., for $z=0$) and one specifying its value for $z\to \infty$. Moreover, the number of dispersion parameters was further doubled due to the uniaxial anisotropy because a separate set of dispersion parameters must be used for the ordinary and extraordinary dielectric functions. The dielectric function given in Equation (31) is a function of eight dispersion parameters, therefore $4\times 8$ parameters were used to describe the profile of the dielectric tensor. If all of these parameters were considered to be independent, the model would be overparameterized. Therefore, it was necessary to reduce the number of parameters sought in processing of the experimental data. It is reasonable to expect that the inhomogeneity of the film and the uniaxial anisotropy are relatively weak effects. Thus, it is possible to assume that some of the dispersion parameters do not exhibit profile (i.e., they have the same value at the upper boundary and for $z\to \infty$) and that some dispersion parameters have the same value for the ordinary and extraordinary dielectric functions. Even though the number of the parameters sought was reduced in this way, there were still enough degrees of freedom left to successfully interpret the measured experimental data.

The dispersion parameters used to describe the studied films could be separated into three groups. The four dispersion parameters $E_{\mathrm{h}}$, $B_{\mathrm{ex}}$, $N_{\mathrm{ut}}$ and $E_{\mathrm{u}}$ constituting the first group were assumed to be without a profile and, moreover, the same values were used for the ordinary and extraordinary dielectric functions. The second group consisting of $E_{\mathrm{g}}$, $A_{\mathrm{ex}}$ and $E_{\mathrm{ex}}$ represented dispersion parameters sought independently for the ordinary and the extraordinary dielectric functions and, furthermore, it was assumed that these parameters exhibit profiles. This means that four versions were considered for each of these parameters, two of them for the values of the ordinary and extraordinary dielectric functions at the upper boundary of the film and two for the values corresponding to the ordinary and extraordinary dielectric functions for $z\to \infty$. The last group consisted of only one parameter, $N_{\mathrm{vc}}$, which specifies the total strength of valence to conduction band transitions. It was assumed that this parameter exhibits profile, thus its value was sought independently for the upper boundary of the film and for $z\to \infty$. However, the same values were used for the ordinary and extraordinary dielectric functions. This is necessary to ensure the fulfillment of the condition in Equation (37), which ensures that the sum rule integral gives the same value for the ordinary and extraordinary dielectric functions (last condition in Equation (27)). The summary of the dispersion parameters describing the inhomogeneous film with uniaxial anisotropy is presented in

Table 1. Summary of dispersion parameters.

Parameters	Exhibit Profile	Different for ${\widehat{\mathit{\epsilon }}}_{\mathbf{o}}$ and ${\widehat{\mathit{\epsilon }}}_{\mathbf{e}}$	Different for Each Sample
$N_{\mathrm{vc}}$	yes	no	yes
$E_{\mathrm{g}}$, $A_{\mathrm{ex}}$, $E_{\mathrm{ex}}$	yes	yes	no
$E_{\mathrm{h}}$, $N_{\mathrm{ut}}$, $E_{\mathrm{u}}$, $B_{\mathrm{ex}}$	no	no	no

.

The multi-sample method [40,41,42] was used to process the experimental data. In the multi-sample method, the experimental data measured on several samples by various techniques (reflectometry and ellipsometry) are processed simultaneously. In our case, the reflectometric and ellipsometric data measured on three samples were processed simultaneously. The samples were prepared under the same deposition conditions but with different deposition times, thus the thickness of the silicon nitride films was different for each film. All dispersion parameters with the exception of the parameters $N_{\mathrm{vc}}$ were assumed to be the same for the silicon nitride films on all three samples. The values of the dispersion parameters $N_{\mathrm{vc}}$ were sought independently for each sample in order to account for small variations in the optical constants of the individual films. This is indicated in the column “different for each sample” in Table 1. Since it was not possible to ensure that the reflectometric and ellipsometric data were measured on the same spots on the samples and the films may exhibit slight non-uniformity in thickness, two film thicknesses were sought for each sample. The first one was used in the processing of the reflectometric data and the other was used in the processing of the ellipsometric data. The values of the profile parameter w used in Equation (39) and the rms values of the heights used to describe the random roughness of the surfaces of the films were sought independently for each sample. By simultaneous processing of experimental data measured on samples with considerably different thicknesses but otherwise very similar properties, it was possible to decrease the correlation among the sought parameters.

The experimental data were processed using the newAD2 open-source software by the least squares method. The objective function was given as

$$S=\sum _{i=1}^3{f}_i^{\mathrm{R}-\mathrm{NIR}}{S}_i^{\mathrm{R}-\mathrm{NIR}}+\sum _{i=1}^3{f}_i^{\mathrm{R}-\mathrm{UVV}}{S}_i^{\mathrm{R}-\mathrm{UVV}}+\sum _{i=1}^3{f}_i^{\mathrm{E}}{S}_i^{\mathrm{E}},$$

(40)

where the index i distinguishes individual samples and $S_i^{\mathrm{R}-\mathrm{NIR}}$, $S_i^{\mathrm{R}-\mathrm{UVV}}$ and $S_i^{\mathrm{E}}$ are the weighted sums of squares corresponding to the reflectometry in the near infrared region, reflectometry in the visible and ultraviolet region and ellipsometry measured on a given sample. The factors $f_i^{\mathrm{R}-\mathrm{NIR}}$, $f_i^{\mathrm{R}-\mathrm{UVV}}$ and $f_i^{\mathrm{E}}$ are the weights for sums of squares of individual measurements. Their values were determined automatically by the algorithm that ensures the contribution to the objective function was approximately the same for each measurement. If these factors were not used, then measurements with large number of points or incorrectly estimated errors could overweigh contributions from the other measurements. This would create a situation in which the influence of some measurements on the sought parameters would be strong while the influence of the other measurements would be negligible.

The experimental data obtained using the imaging spectroscopic reflectometry were processed using the same structural and dispersion models as the reflectometric and ellipsometric data. The thickness of the films was sought independently for each pixel using the least squares method. The values of all the other structural and dispersion parameters were fixed in the values found in the processing of the reflectometric and ellipsometric data.

4. Results and Discussion

The results determined by the simultaneous processing of experimental data of all three samples are presented in

Table 2. The values of the structural parameters.

Parameter			Sample 1	Sample 2	Sample 3
deposition time	t	[min]	30	45	90
thickness ellipsometry	$d_{\mathrm{e}}$	[nm]	$111.21\pm 0.05$	$167.10\pm 0.08$	$323.9\pm 0.1$
thickness reflectance	$d_{\mathrm{r}}$	[nm]	$112.12\pm 0.06$	$167.91\pm 0.07$	$320.80\pm 0.09$
roughness (rms)	$\sigma$	[nm]	$3.0\pm 0.2$	$4.6\pm 0.2$	$6.3\pm 0.2$
profile parameter	w	[nm]	$39.3\pm 1.8$	$65.6\pm 2.8$	$57.6\pm 2.6$

. Thickness values corresponding to the ellipsometry and reflectometry for the individual samples are mutually very close, which means that the thickness non-uniformity of the samples is very small. The thickness values of individual samples increase approximately linearly with the deposition times. The rms values of the heights $\sigma$ of the random roughness are also increasing with the deposition time.

The values of the dispersion parameters are presented in

Table 3. The values of the dispersion parameters.

Parameter
$E_{\mathrm{h}}$	[eV]		$70.0\pm 1.6$
$B_{\mathrm{ex}}$	[eV]		$3.90\pm 0.03$
$N_{\mathrm{ut}}$	[eV${}^2$]		$37.8\pm 2.6$
$E_{\mathrm{u}}$	[eV]		$0.405\pm 0.002$
Parameters without profile which are identical for the ordinary and extraordinary dielectric functions
Parameter			Upper Boundary		$\mathbf{z}\to \infty$
$N_{\mathrm{vc}}$	[eV${}^2$]	Sample 1	$401.6\pm 3.2$		$453.7\pm 3.4$
		Sample 2	$398.0\pm 3.5$		$437.1\pm 3.3$
		Sample 3	$381.6\pm 2.9$		$430.8\pm 3.2$
Parameters with profile which are identical for the ordinary and extraordinary dielectric functions
			Upper Boundary		$\mathbf{z}\to \infty$
Parameter			Ordinary	Extraordinary	Ordinary	Extraordinary
$E_{\mathrm{g}}$	[eV]		$2.65\pm 0.03$	$2.43\pm 0.0$	$2.77\pm 0.03$	$2.55\pm 0.02$
$A_{\mathrm{ex}}$			$0.21\pm 0.03$	$1.1\pm 0.1$	$11.8\pm 1.1$	$2.4\pm 0.1$
$E_{\mathrm{ex}}$	[eV]		$10.31\pm 0.05$	$9.15\pm 0.07$	$11.20\pm 0.02$	$8.65\pm 0.05$
Parameters with profile which are different for the ordinary and extraordinary dielectric functions

. This table is divided into three parts. In each part, the values of the dispersion parameters corresponding to groups introduced in Table 1 are listed. The spectral dependencies of the optical constants of the inhomogeneous films for the ordinary and the extraordinary waves at the upper and lower boundaries are presented in Figure 2 and Figure 3. Note that the complex refractive indices are defined as ${\widehat{n}}_{\mathrm{o}}=\sqrt{{\widehat{\epsilon }}_{\mathrm{o}}}$ for the ordinary wave and as ${\widehat{n}}_{\mathrm{e}}=\sqrt{{\widehat{\epsilon }}_{\mathrm{e}}}$ for the extraordinary wave. They could be written as ${\widehat{n}}_{\mathrm{o}}=n_{\mathrm{o}}+\mathrm{i}{k}_{\mathrm{o}}$ and ${\widehat{n}}_{\mathrm{e}}=n_{\mathrm{e}}+\mathrm{i}{k}_{\mathrm{e}}$, where $n_{\mathrm{o}}$ and $n_{\mathrm{e}}$ are the real refractive indices and $k_{\mathrm{o}}$ and $k_{\mathrm{e}}$ are the extinction coefficients.

The profiles of the optical constants corresponding to the ordinary and extraordinary waves are displayed in Figure 4. In this figure, it is evident that the profiles of the inhomogeneous films on Samples 2 and 3 are very similar while the profile exhibited by Sample 1 is somewhat different from the other two samples. This is also visible in Figure 3 on the spectral dependencies of the optical constants at the lower boundary, which are almost the same for Samples 2 and 3 while the refractive indices of Sample 1 are little higher. Moreover, the values of the profile parameters w are almost the same for Samples 2 and 3 but the value determined for Sample 1 is lower than the values determined for the other two samples. However, we should note that the spectral dependencies of the optical constants at the upper boundary exhibit different behavior. The refractive indices corresponding to Samples 1 and 2 are almost the same while the refractive index corresponding to Sample 3 is little lower. Nevertheless, the differences between the optical constants at the upper and lower boundary are relatively small. This means that the inhomogeneity in the optical constants of all three samples exhibits practically the same behavior. The optical constants change in the vicinity of the upper boundary while the lower parts of the films are almost homogeneous. In Figure 2, Figure 3 and Figure 4, it is evident that the refractive indices are smaller at the upper boundaries and higher at the lower boundaries. This is probably connected with the technological process of creating the films.

We also tried to use a similar model with exponential profile in which the optical constants changed near the lower boundary and the upper parts of the films were almost homogeneous. However, this model could not correctly describe the experimental data.

In Section 2.3, the formulae derived within the SDT use the approximation that the optical constants in the range of roughness are constant. Figure 4 shows that the optical constants exhibit non-negligible changes within this range. Thus, the approach that we used to include the influence of the roughness on the optical quantities should be considered to be only a rough approximation. Since the roughness is slight (the rms values of the heights are small), the use of this rough approximation is permissible. This is also supported by the very good fits of the experimental data plotted in Figure 5 and Figure 6.

The thickness of the silicon single crystal substrate was determined using the micrometer screw gauge as 0.585 mm. The thicknesses of the native oxide layers on the back sides of the substrates were determined on the basis of the ellipsometric data measured on the back sides at angle of incidence 65${}^{\circ }$. The determined values were 3.09 nm for Sample 1, 2.74 nm for Sample 2 and 3.38 nm for Sample 3.

The maps of local thicknesses obtained by ISR are presented in Figure 7 and Figure 8. From these maps, it is clear that the certain non-uniformity in thickness occurs near the edges of the samples. In the central parts of the samples, the inhomogeneous films are practically uniform in thickness. Since the measurements were performed in these uniform parts, it was not necessary to include the thickness non-uniformity into our structural model.

The topography of the surface of Sample 3 obtained by AFM is introduced in Figure 8. The rms value of the heights determined from the AFM image is $7.50\pm 0.04$ nm. This is in a relatively good agreement with the rms value of the heights determined by means of the optical method $6.3\pm 0.2$ nm.

5. Conclusions

In this study, the optical characterization of the non-stoichiometric silicon nitride films prepared by the reactive magnetron sputtering was performed. To interpret the experimental data, it was necessary to assume that the characterized films exhibit several defects. The original model combining these defects had to be developed. These defects are the inhomogeneity in the optical constants, uniaxial anisotropy with optical axis perpendicular to the boundaries and the roughness of the upper boundaries of the films. The inhomogeneity is represented by profiles of the optical constants across the films.

The simultaneous inclusion of inhomogeneity and anisotropy was achieved using the approximation in which the inhomogeneous layer was approximated by a multi-layer system consisting of large number of thin homogeneous layers. Since the inhomogeneous film was anisotropic, the individual layers in this multi-layer system were also assumed to be anisotropic. The optical quantities of the multi-layer system were calculated using the Yeh matrix formalism. The random roughness of the upper boundaries of the studied films was included using the SDT under the assumption that the optical constants exhibit only small changes in the range of this roughness.

The optical characterization was performed on the basis of the ellipsometric and reflectometric data measured within the spectral range 0.6–6.3 eV. The inclusion of all defects required a model depending on a large number of parameters. To reduce the correlation among the parameters sought within the data processing, the multi-sample method was utilized. The experimental data measured on three samples prepared with different deposition times (i.e., differing in thickness) were processed simultaneously.

It was found that the films exhibited profiles of the optical constants, which can be described using the exponential function so that the pronounced changes of the optical constants occur in the region near the upper boundary. Within the substantial part of the film, the optical constants were practically constant. The evident uniaxial anisotropy with the optical axis perpendicular to the boundaries corresponded to different values of the optical constants for the ordinary and extraordinary waves. Furthermore, the rms values of the heights of the roughness were found in the range of nanometers for all three films.

The ISR was utilized to investigate the thickness non-uniformity of the films. It was found that the films are uniform in thickness along the substantial parts of their area; slight thickness non-uniformity was observed only in the regions near the edges.