## Abstract

Fresnel reflectance and transmittance coefficients of a thin film system consisting of an arbitrary number of layers are expressed explicitly in the form of a power series of Fresnel coefficients for individual boundaries and phase terms for the individual films. The series is based on the evaluation of all possible paths light can pass through the system. However, the series is written as consolidated, i. e. all paths corresponding to the same powers are represented using a single term in the series, with multiplicity which is a simple product of binomial coefficients. This result is used to express the normal reflectance of a thin film system with arbitrarily correlated randomly rough boundaries and it is shown that such approach can be computationally efficient in practice.

© 2014 Optical Society of America

## 1. Introduction

Interaction of light with a layered medium, such as a system of thin films, can be viewed as a series of interactions (reflections and refractions) at layer boundaries and propagation through the individual films. However, this view is mostly relegated to introductory optics textbooks where the enumeration of all the interactions may be used for the derivation of reflectance and transmittance of a single thin film. In most practical calculations in optics of thin films a matrix formalism is used instead. Several such formalisms exists, e. g. a 2 × 2 or 4 × 4 formalisms [1] based on the formalism developed by Jacobson [2], or an admittance-based approach [3], nevertheless, they all express the relation between incoming and outgoing light wave amplitudes at one end and the other end of the thin film system using a matrix. Another frequently used approach to multi-layer systems is the recursive or inductive application of formulae for a single layer [3–6].

Both the matrix formalism and recursive formulae are excellent for the calculation of reflectances and transmittances of multi-layer systems with smooth boundaries. However, in the case of systems with rough boundaries the situation is different. The only theoretical approach that fits well into the matrix formalism is the effective medium approximation (EMA) [7–9]. This approach is very popular, also due to its simplicity, however, it is valid only for very fine roughness with high spatial frequencies [1]. More advanced approaches used in optical characterisation of rough thin films and their systems include in particular the scalar diffraction theory (SDT) [10–17] and perturbation theories [18–26] such as the Rayleigh–Rice theory (RRT) [27–32]. An overview of the theoretical approaches can be found in [23, 33, 34]. New methods of calculation of light scattered from rough surfaces and thin films are still being actively developed. However, the SDT and second-order RRT remain the state of the art in the field of ellipsometric and spectrophotometric characterisation of rough thin film systems using specularly reflected light. Both SDT and second-order RRT are able to describe roughness outside of the validity range of EMA, each for a different range of magnitudes and spatial frequencies of the roughness. But both also express the reflectances and transmittances globally for the entire thin film system and are not factorable to matrices. This means that they cannot be included within the matrix formalism. Note that a simple matrix approximation of the RRT, applicable in certain special cases, has been formulated [1, 30]. In the case of SDT, an approximate approach was published, based on averaging of individual elements of the matrices entering the matrix formalism [14, 15]. However, this approximation is too rough and often leads to erroneous results. An improved but still approximate method was based on the approximation 〈*a/b*〉 ≈ 〈*a*〉/〈*b*〉 in the recursive formulae for reflectance of a multi-layer thin film [35, 36] (〈...〉 denotes averaging).

In the SDT, it is advantageous to express the reflection coefficients of a layered system using a power series in the Fresnel coefficients and phase terms for individual boundaries and layers instead [11,12,37,38]. However, neither the matrix formalism nor the recursive formulae express the Fresnel coefficients of multi-layer systems as power series. They express the coefficients as relatively complicated rational functions of Fresnel coefficients of individual boundaries and phase terms describing the propagation inside the layers. The sequence-of-interactions view, on the other hand, leads to a power series suitable for the SDT. As noted above, writing the series for a single layer is trivial. Utilisation of the series within the SDT leads to the well-known expressions for the reflectance of rough thin films [12]. Multiple layers, however, represent a challenge.

All possible sequences of reflections and refractions that start with light wave falling onto the uppermost boundary from the ambient and ending with the light wave leaving the uppermost or lowermost boundary towards the ambient were enumerated by Butler in a mathematical journal probably not read by many physicists [39]. The connection to Catalan paths was also noted. Catalan paths for systems from two to six layers are enumerated for instance in the on-line encyclopaedia of integer sequences [40–44]; they converge for an infinite number of layers to Catalan numbers [45].

For any multi-layer system the total number of possible ‘paths’ grows exponentially with the number of interactions. However, many of the paths lead to the same combinations of Fresnel coefficient and phase term powers. An illustration of two simple sequences of reflections and refractions that produce equivalent terms in the series is shown in parts (a) and (b) of Fig. 1. Therefore, it would be useful to write a consolidated power series which contains each such term only once, with the corresponding multiplicity as a mere combinatorial factor. And it would be useful to group these terms further by the powers of phase terms as this is also useful in SDT calculations. A path that corresponds to different powers of Fresnel coefficients but the same phase term powers is illustrated in Fig. 1(c).

Both consolidations are possible and the construction of such series will be presented in this work. Furthermore, it will be shown how this series expansion can be applied within SDT to calculation of normal-incidence reflectance of a thin film system with rough boundaries.

Note that the path summation approach was used for the expression of reflection and transmission coefficients of individual inhomogeneous layers [46], although only the lowest order terms were expressed explicitly. In this work the full infinite series containing paths of arbitrary lengths will be formulated. The complete series is equivalent to the expressions obtained using the matrix formalism or recursive formulae.

## 2. Construction of the series

Since any practical calculation employing summation of the series may include only a finite number of terms, the series will be constructed from the beginning as truncated to a certain number of interaction at boundaries. Letting the number of interactions grow to infinity in order to obtain the infinite series is then trivial, whereas the limiting the infinite series to a specific number of interactions could be considerably less straightforward. It should be noted the ordering by the number of interactions is used because it is convenient for the construction of the infinite series. In practical calculations, the series can be truncated to a finite number of terms differently as will be discussed at the end of paragraph 2.2.

#### 2.1. Parametrisation

The media in the layered system will be indexed *j* = 0, 1,...,*L*, with 0 denoting the ambient, *L* denoting the last medium which is presumably the substrate. Index *j* will be referred to as medium depth for brevity. The reflection and transmission coefficients and phase terms will be then denoted *r̂ _{j}*,

*r̂′*

_{j}*t̂*,

_{j}*t̂′*and exp(i

_{j}*X̂*) as displayed in Fig. 2, where

_{j}*N̂*,

_{j}*h*and

_{j}*ϑ̂*denote the complex refractive index, thickness, and complex angle of incidence corresponding to medium

_{j}*j*; and

*λ*is the wavelength. No specific relations among the coefficients will be assumed because it would not bring any significant simplification while it would prevent application of the results in other situations, namely when the layer boundaries are not simple ideal boundaries.

A light path which is continuous and starts and ends in medium 0 contains an odd number *m* of interactions at the boundaries. It is more useful to introduce path length *p* using the formula

*p*= 0 corresponds to a single reflection at the 0–1 boundary. The total power of all Fresnel coefficients in the term corresponding to this path is 2

*p*+ 1 and the total power of all phase terms is 2

*p*. Beside path length we will also define path depth

*d*which is equal to the depth of the deepest medium the path passes through (so

*d*≤

*L*). Any two paths that differ in either

*p*or

*d*are obviously different.

A suitable set of parameters describing a path of depth *d* is then formed by non-negative integers *m*_{0}, *m*_{1}, . . . , *m _{d}* called medium counts and

*v*

_{0},

*v*

_{1}, . . .

*v*called visit counts.

_{d}The medium count *m _{j}* expresses how many times the light:

- passes through medium
*j*downwards, - passes through medium
*j*upwards, - is incident on the
*j*–(*j*+ 1) boundary from the side of medium*j*, - is incident on the (
*j*− 1)–*j*boundary from the side of medium*j*.

*m*is also equal to the power of exp(2i

_{j}*X̂*) in the term corresponding to paths described by this value of

_{m}*m*. The medium count can also be interpreted as the number of times the light is reflected off the entire subsystem from medium

_{j}*j*+ 1 onward, which will be useful below. The path length is related to

*m*as follows and for all paths and all layer systems it holds

_{j}*m*

_{0}= 1.

The visit count *v _{j}* expresses how many times the light:

- visits the subsystem from medium
*j*downward, - passes through the (
*j*− 1)–*j*boundary from the side of medium*j*− 1, - passes through the (
*j*− 1)–*j*boundary from the side of medium*j*.

*j*− 1 and segments contained in media from

*j*downward. Visit count

*v*is the number of the latter segments. For all paths and all layer systems it holds

_{j}*v*

_{0}= 1.

#### 2.2. Parameter ranges

Medium and visit counts are bound by a set of simple inequalities that follow from continuity of the path (i. e. that it cannot skip any layers). First, for all *j* ≤ *d* it holds

*j*>

*d*it holds

*v*=

_{j}*m*= 0. Once the light enters medium

_{j}*j*it must pass through it at least once which leads to the inequality

*v*≤

_{j}*m*. Finally, since visits of medium

_{j}*j*+ 1 must be separated by at least one pass through medium

*j*it must hold

*v*

_{j}_{+1}≤

*m*. The preceding two relations can be written summarily for

_{j}*v*or

_{j}*m*in a more compact form

_{j}*m*

_{1},

*m*

_{2}, . . . ,

*m*and accompanying set of

_{d}*v*

_{1},

*v*

_{2}, . . . ,

*v*bound by these conditions corresponds to at least one possible path.

_{d}As can be seen with the help of Fig. 3, if the path

- enters medium
*j*from medium*j*− 1 the same number of times as it leaves medium*j*to medium*j*− 1, - enters medium
*j*+ 1 from medium*j*the same number of times as it leaves medium*j*+ 1 to medium*j*, - leaves the (
*j*− 1)–*j*boundary downwards the same number of times as it is incident on the*j*–(*j*+ 1) boundary from the side of medium*j*, and - leaves the
*j*–(*j*+ 1) boundary upwards the same number of times as it is incident on the (*j*− 1)–*j*boundary from the side of medium*j*then it is always possible to ‘match the arrows’ in medium*j*close to the (*j*− 1)–*j*boundary with those close to the*j*–(*j*+ 1) boundary. Hence a path without any loose ends can be always formed provided that*m*−_{j}*v*and_{j}*m*−_{j}*v*_{j}_{+1}are non-negative – which corresponds exactly to conditions (5) (if the matching is done in order left to right no cycles are formed).

Using the second inequality in Eq. (5), ranges of visit counts can be easily expressed if the medium counts are given. The overall summation scheme for paths of positive depth can thus be constructed as follows:

where*P*denotes the maximum path length which can be either specified explicitly or determined implicitly by a termination condition, e. g. by terminating the series once the absolute value of the entire length-

*p*term falls below a certain small value. The vectors

**m**and

**v**were introduced as shorthand for

*m*

_{1},

*m*

_{2},...,

*m*

_{d}_{−1}and

*v*

_{2},

*v*

_{3},...,

*v*. The

_{d}**m**-summation is then

**v**summation is

*M*denotes the sum of medium counts for media from 1 to

_{j}*j*

*m*=

_{d}*p*−

*M*

_{d}_{−1}which follows from Eq. (3). No summation is thus carried out over

*m*because it is uniquely determined by the path length and other

_{d}*m*. Similarly, no summation is carried out over

_{j}*v*

_{1}because

*v*

_{1}≤ min(

*m*

_{0},

*m*

_{1}) = 1, hence

*v*

_{1}= 1.

It should be noted that the ordering by path length in Eq. (6) was the only reason why any relations between ranges of individual *m _{j}* appeared in Eq. (7). If this requirement is removed the upper limits of

*m*summations can be chosen differently, e. g. each can be chosen independently based on physical considerations. Therefore, a criterion based on the behaviour of individual reflection and refraction coefficients can be used, taking into account more reflections in specific selected layers.

_{j}The above scheme does not include the singular term corresponding to a single reflection off the uppermost boundary (*p* = 0, *d* = 0). It could be formally included in the scheme but only at the cost of impairing its clarity. Since it would still represent a special case, only made more obscure, the result would not translate well into a computer implementation. Hence it was excluded and must be added separately.

#### 2.3. Coefficient and phase term powers

The powers of Fresnel coefficients and phase terms can be expressed trivially with the help of Fig. 3. Layer *j* contributes the factor

Depth 0 corresponds to the single term ${P}_{0}={r}_{1}^{{m}_{0}}={\widehat{r}}_{1}$ representing the reflection off the uppermost boundary.

Depth 1 represents reflections within the uppermost layer, i. e. the classic series

*v*

_{1}= 1 was used to obtain the latter expression from the former which corresponds formally to Eq. (10). The following formulae will already be written with this simplification.

Adding the contribution of layer 2 results in the expression of depth-2 terms:

#### 2.4. Multiplicities

The only missing piece is the multiplicity factor for each term in the series, i. e. the number of possible paths that give raise to each unique combination of Fresnel coefficient and phase term powers. Again, we will refer to Fig. 3 from which it can be seen that in each layer two independent choices are made:

- distribution of the
*v*visits of medium_{j}*j*among the*m*_{j}_{−1}incidences of light on (*j*− 1)–*j*boundary from the side of medium*j*− 1, - distribution of the
*m*downward passes of light through medium_{j}*j*among the*v*visits of this medium._{j}

The first distribution can be done

ways while the second can be done ways. Therefore, layer*j*contributes to the term multiplicity with an independent factor that will be denoted

*F*(

*v*,

_{j}*m*,

_{j}*m*

_{j}_{−1}), i. e.

*C*for a depth-

_{d}*d*summation term can be obtained simply by multiplying these factors for all

*d*layers:

#### 2.5. Second consolidation and summation algorithm

The order of summation in Eq. (6) permits grouping all identical phase terms together as follows:

*Q̂*(

**m**) is a homogeneous polynomial of degree 2

*p*+1 in Fresnel coefficients. The ranges of

*m*, which are omitted here for brevity, are the same as in Eq. (6) where they were written out explicitly. The polynomial

_{j}*Q̂*(

**m**) is constructed using the tail part of Eq. (6), expression (10) with phase terms excluded, and formula (17):

*v*. Conversely, terms ${\widehat{r}}_{j+1}^{{m}_{j}-{v}_{j+1}}$ in Eq. (10) were pulled into one level deeper summations since they belong to layer

_{j}*j*conceptually but depend on

*v*

_{j}_{+1}, not on

*v*. The product from

_{j}*j*= 2 to

*d*should be understood as empty (equal to 1) for

*d*= 1.

The full factorisation of Eq. (20) means that if the Fresnel coefficients can be treated independently, the number of arithmetic operations necessary to calculate the sum (20) grows only as *p* instead of the *p ^{d}/d*! growth which corresponds to the fully expanded sum. If the further mathematical operations can be only carried out with complete terms, the expansion is straightforward: the summations over

*v*are performed sequentially, in a nested manner indicated in Eq. (6), instead of factoring the sums. In paragraph 2.6 where the number of terms is discussed, the latter, more pessimistic, variant will be assumed. It should be noted that both variants coincide for a two-layer system because only one

_{j}*v*summation (over

_{j}*v*

_{2}) is done.

#### 2.6. Number of terms

As noted in the introductory remarks, the number of possible paths grows exponentially with path length. Specifically for paths of depth 2, 3 and 4 it grows as 2^{p−1}, *F*_{2p+1} and (3*p* − 1)/2+1, respectively, where *F _{n}* denotes the

*n*-th Fibonacci number [40–42].

In contrast to the number of possible paths, the number of terms in the consolidated series grows polynomially with *p*. For paths of depth 2 and *p* ≥ 2 it immediately follows from parameter ranges given by Eqs. (4) and (5) that the number of terms is

*p*

^{4}/72. Therefore, the number of terms in the series including all paths with length up to

*p*, inclusively, grows as

*p*

^{3}/12 for a two-layer system and as

*p*

^{5}/360 for a three-layer system. Generally, for an

*L*-layer system it grows as

*p*

^{2}

^{L}^{−1}. This means that while the combinatorial explosion may be much slower asymptotically in the consolidated series than in the expanded series, it is still fast if the system consists of many layers.

Furthermore, it is not possible to terminate the series after a small number of terms (small *p*) if *L* is large and if any contribution of the deeper layers is to be included because the *j*-th medium contributes only to paths of length *j* or more.

Both properties limit the usability of the series expansion, even consolidated, in the case of many-layer systems. However, in the case of systems consisting only of several layers it has a certain potential.

The reduction of the number of terms by consolidation can be observed in Table 1 which compares the number of terms corresponding to all paths of length *p* in systems consisting of two to five layers (paths of any depth are counted, not just those that of maximum possible depth). It can be seen that twofold reduction occurs if paths of length at least 5 are considered so the presented scheme does not lead to a significant simplification if only a handful of terms is taken into account. For paths longer than 5, however, the gain from consolidation is considerable and it increases rapidly with increasing *p*.

#### 2.7. Transmission

The summation scheme formulated for reflection can be adapted for transmission in a fairly straightforward manner. All paths starting with light wave falling onto the uppermost boundary from the ambient and ending with the light leaving the lowermost boundary to the ambient have depth *d* = *L* so summation over *d* is avoided in formulae such as (6) and *d* is replaced with *L*. The number of terms in the consolidated series grows somewhat less rapidly, though still as *p*^{2}^{L}^{−1}.

The meaning of *m _{j}* and

*v*changes slightly: they still denote the numbers of downward passes and refractions, however, the numbers of upwards passes and refractions are now

_{j}*m*− 1 and

_{j}*v*− 1, respectively. Figure 3 remains largely unchanged, only with upward arrow counts replaced with

_{j}*v*− 1 and

_{j}*v*

_{j}_{+1}− 1. Definition (3) of path length

*p*can be retained, however, the total number of interactions is 2

*p*+ 1 −

*L*now and the total power of phase terms is similarly 2

*p*−

*L*. Condition (4) holds for all 0 ≤

*j*≤

*L*+ 1 and condition (5) remains unchanged.

The contribution of one layer to powers of Fresnel coefficients (10) changes to

*r̂*

_{1}disappears.

Concerning multiplicities, the reasoning leading to expression (15) remains valid, whereas Eq. (14) needs to be modified. To obtain a path leaving the layer system through the lowermost boundary, the last incidence of light must result in a refraction, not reflection. Therefore, only the remaining *v _{j}* − 1 visits can be freely chosen among the remaining

*m*

_{j}_{−1}− 1 incidences and the factor

*F*(

*v*,

_{j}*m*,

_{j}*m*

_{j}_{−1}), expressed by Eq. (16) for reflection, becomes

## 3. Application to scalar diffraction theory for normal incidence

The summation scheme presented above can be applied in a relatively straightforward manner to SDT for the normal reflectance of a thin film system with rough boundaries exhibiting arbitrary cross-correlations. It was shown [35,47,48] that the normal reflectance *R* is expressed

**u**=

**u**(

*x,y*) is the vector with

*L*+ 1 components representing the heights of irregularities of boundaries 1, 2,...,

*L*+1 at the Cartesian coordinates (

*x,y*) in the sample plane;

*u*

_{1}is the first component of

**u**. The symbol

*r̂*(

**u**) denotes the Fresnel reflection coefficient of the corresponding system of smooth films with thicknesses

*h*replaced with

_{j}*h*+

_{j}*u*−

_{j}*u*

_{j}_{+1}(putting

*u*

_{0}≡

*u*

_{L}_{+1}≡ 0 to simplify the notation). The thickness

*h*must be now understood as the mean thickness of the

_{j}*j*-th film. The factor

*v*is equal to where

*N*

_{0}is refractive index of the ambient. The averaging over

**u**, denoted with angle brackets 〈...〉, is carried out within the illuminated area. Under the assumption that the roughness corresponds to a stationary random process (in the wide sense), the sample-plane averaging can be replaced by statistical averaging with the multi-dimensional probability density function of

**u**, which will be denoted

*ρ*(

**u**):

**u**is

^{T}and

^{−1}denote transposition and matrix inversion, respectively. The matrix

**S**is a symmetrical positive semi-definite matrix with components where is the rms value of height irregularities of the

*j*-th boundary and the cross-correlation coefficients between

*i*-th and

*j*-th boundary are equal to [35]

The integral (26) can be evaluated for each term in the consolidated series (19) using the well known formula for multi-dimensional Laplace transform of a general Gaussian function

**f**is an arbitrary vector) and the following series is obtained:

*ϑ̂*= 1 and thus formula (1) reduces to The factor

_{j}*Ĥ*(

*m*) is given by

*m*

_{0}≡ 1 and

*m*= 0 for

_{j}*j*>

*d*to include the edge terms

*D̂*

_{1}and

*D̂*

_{d}_{+1}in a single expression. For a single layer-system the sum (32) coincides with the well-known sum [12].

Formula (32) is applicable to thin film systems exhibiting arbitrary correlations between the Gaussian boundary roughness of individual layers. Several specific types of rough thin film systems have been discussed in the literature, for some of them simplifications can be achieved due to a special form of **S** [38, 49].

## 4. Numerical example

Alternative methods of evaluation of the SDT for normal reflectance of rough thin film systems exist. One possibility is the utilisation of the approximate formulae [35]. For a small number of layers, one can also try to evaluate the integral in Eq. (24) by a direct integration. Furthermore, it is possible to attempt an explicit tracing of the largest contributions to the result instead of the indiscriminate inclusion of all contributions, whether large or small, in the series.

To illustrate the trade-offs between precision and calculation efficiency of the various approaches, their performance was compared for a simple example of a three-layer thin film system. In specific, the following procedures were compared:

- Direct evaluation of Eq. (24) by means of a numeric quadrature. The change of variables
**u**=**Uv**, where**UU**^{T}=**S**, transforms the probability density function to a spherically symmetrical multi-dimensional Gaussian with unit dispersion in**v**. The resulting integral was evaluated using a Cartesian product of one-dimensional Hermite–Gauss quadratures [50], which represents a fairly efficient integration method. The reflectance coefficient*r̂*(**u**) in the integrand was expressed using the 2 × 2 matrix formalism. Note that using the recursive formulae would lead to practically identical results because these two approaches are very similar from the point of view of computational efficiency. - Utilisation of consolidated series (32) presented in this work.
- Ray tracing. This means the light paths were constructed recursively. All possible reflections and transmissions were raced, stopping the recursion when the estimated contribution to the reflection coefficient of the entire system dropped below a prescribed threshold. This stopping criterion permits the exclusion of combinations of reflections and transmissions that contribute little to the result, which is not easy using the series. Formula (34) was then used for each generated path leaving the upper boundary to include the effect of roughness. The exclusion of many small terms that contribute little to the result might be expected to lead to a better computation efficiency than the utilisation of the series. On the other hand, the effect of the combinatorial explosion can be much stronger in this approach.
- Application of the approximate formulae [35].
- Disregarding the roughness, i. e. using the formulae for a system of thin films with smooth boundaries.

The simulated thin film system is schematically depicted in Fig. 4. It consists of three layers in the following order (from top): Si_{3}N_{4} (150 nm), SiO_{2} (130 nm) and Si_{3}N_{4} (100 nm). The substrate is formed by crystalline silicon. Optical constants of all materials were taken from standard tables [51–53].

The boundary roughness of a three-film system is described by ten parameters, namely the ten independent elements of matrix **S** or, equivalently, by four *σ _{j}* and six

*C*. For the simulated thin film system, they were constructed based on the notion of a gradual built-up of roughness during deposition. It was assumed that the

_{jk}*k*-th layer adds independent roughness contribution

*ε*(

_{k}*x*,

*y*) of the same rms value

*σ*on top of the

*k*+ 1-th layer. The heights of irregularities were thus related:

*σ*

_{L}_{+1}and

*σ*denote the rms of heights of the lowermost boundary and the roughness added at each subsequent film boundary, respectively. It can be seen that this results in the following values of

**S**matrix elements: Values

*σ*

_{L}_{+1}= 10nm and

*σ*= 5nm were used in the calculation.

The calculated spectral dependence of reflectance in the range from 200 to 800 nm is plotted in Fig. 5 for all procedures. Of course, the series, numeric quadrature and ray tracing represent different numerical methods for the calculation of the same result. Therefore, they all permit, in principle, calculation with an arbitrarily high precision (within the SDT). Hence one curve is used to represent all in the figure as they are indistinguishable when calculated to a sufficient precision. The approximate formulae provide reflectance values that are relatively close to the values obtained by precise SDT. Nevertheless, deviations up to several percents are noticeable, namely in the short-wavelength region. As can be also seen from the figure, the formulae for a system of thin films with smooth boundaries would be inadequate for a thin film system similar to the simulated system.

It was already noted above that the series, numeric quadrature and ray tracing are not fixed approximations (unlike the approximate formulae) and allow a trade-off between the precision and speed of computation. It is, therefore, interesting to compare them also from this point of view. Choosing a specific precision is relatively straightforward for the series because it suffices to terminate the series at path length *p* for which the entire *p*-length term falls below a prescribed value.

In the case of Hermite–Gauss quadratures, the precision control is more difficult as it is not obvious how the number of quadrature points translates to precision of the result. An adaptive multi-dimensional quadrature could provide an error estimate. For the purpose of this illustration, a simpler approach was used: The reflectance was evaluated using quadratures with different numbers of points from 3 to 25, the computation time was measured and the precision was determined ex post by finding the maximum difference from a reference result obtained by evaluation of the series to full double precision. For consistency, the computation time and precision of the other methods was measured in the same manner.

The result of this benchmark is shown in Fig. 6. The time plotted in this figure is the mean computation time per one spectrum point for the three-layer system discussed above, measured on author’s Athlon 64 PC. The approximate formulae occupy the lower right corner of the graph, corresponding to short computation time but low precision. The performance of numeric quadrature may be perhaps acceptable in applications requiring only about three significant digits, even for this precision range the consolidated series represents a considerably better choice. A notable property of the series is the slow increase of computation time required for higher numbers of significant digits. This is illustrated by extending the abscissa to precision 10^{−12} that would be, of course, absurd for characterisation of rough thin film systems in practice. The slow computation time growth is in a stark contrast with the behaviour of ray tracing. Although ray tracing offers a reasonable speed–precision trade-off at lower precisions (10^{−3}), the computation time grows rapidly with increasing precision due to the combinatorial explosion discussed in paragraph 2.6. For the number of significant digits that would be required for instance for numerical differentiation in data fitting, the computation time is already unacceptable. Furthermore, precision control is more difficult for ray tracing since the recursion stopping threshold does not translate directly to precision of the result. Therefore, it can be concluded that of the five approaches tested, the consolidated series preformed best. In fact, it turned out to be the only method capable of calculation of the reflectance to the full double precision in a reasonable time.

Finally, a peculiar feature of the series should be noted. It becomes more efficient with increasing roughness, whereas the opposite holds for the application of the matrix formalism and integration. This is because the exp(−*x*^{2})-like factors in the terms decrease more rapidly with increasing roughness. For integration, however, increased roughness means larger ranges of **u** and this implies more abscissa points are required to achieve the same precision of the numeric quadrature. This effect is illustrated in Fig. 7, where the computation efficiencies of both methods are compared in dependence of the roughness magnitude. The computations were performed for rough thin film systems that were identical as in the preceding example, except for the overall roughness magnitude. The overall roughness magnitude was varied by scaling all *σ* values by the same factor. The value of the lowermost roughness *σ _{L}*

_{+1}, which is plotted on the abscissa, was varied in the interval [1.25, 15]nm. The mean computation times plotted in the figure correspond to times needed for obtaining the entire reflectance spectra with precision at least 10

^{−6}. Again, the precision was determined as the difference to the values calculated using the series to the full double precision. From Fig. 7 it can be seen that although the efficiency of the consolidated series is decreasing with decreasing roughness, the series is more efficient for all roughness values above approximately 2 nm for this thin film system. Of course, this threshold value depends on the structure and optical constants of the thin film system and would differ for other systems. For zero roughness, i. e. a smooth system, the series is not very efficient. However, for such system, one would use the matrix formalism or recursive formulae in practice.

## 5. Conclusion

Fresnel reflectance and transmittance coefficients of the system of an arbitrary number of thin films were expressed in the form of an explicit power series in Fresnel coefficients for individual boundaries and phase terms for the individual films. The series was constructed by considering all possible paths of light through the multiple layers and it is ordered by the length of this path. It was, however, constructed as consolidated, i. e. all paths leading to the same powers are represented using a single term in the series and a numeric multiplicity. Notably, the multiplicities were found to be simple products of binomial coefficients and their evaluation does not require any advanced number-theoretic methods. This fundamental result is interesting from a mathematical point of view, however, the series expression is not a mere mathematical curiosity. Since the expression of Fresnel coefficients of a thin film system as a power series is advantageous for the calculation of optical quantities of rough thin film systems using the scalar diffraction theory, the series was utilised to obtain formulae for the normal reflectance of a thin film system with arbitrarily correlated Gaussian rough boundaries. Using a numerical example, it was illustrated for a three-layer system that such approach is viable and can be used for an efficient calculation of optical properties of rough thin film systems using the scalar diffraction theory.

## Acknowledgments

This work was supported by the projects “CEITEC – Central European Institute of Technology” (CZ.1.05/1.1.00/02.0068) and “R&D center for low-cost plasma and nanotechnology surface modifications” (CZ.1.05/2.1.00/03.0086) from European Regional Development Fund.

## References and links

**1. **I. Ohlídal and D. Franta, “Ellipsometry of thin film systems,” in *Progress in Optics*, E. Wolf, ed. (Elsevier, 2000), vol. 41, pp. 181–282. [CrossRef]

**2. **R. Jacobson, “Light reflection from films of continuously varying refractive index,” in *Progress in Optics*, E. Wolf, ed. (North-Holland, 1966), Vol. 5, pp. 249–286.

**3. **Z. Knittl, *Optics of Thin Films* (John Wiley, 1976).

**4. **P. Rouard, “Etudes des propriétés optiques des lames métalliques très minces,” Ann. Phys. **7**, 291–384 (1937).

**5. **A. Vašíček, “The reflection of light from glass with double and multiple films,” J. Opt. Soc. Am. **37**, 623–634 (1947). [CrossRef]

**6. **A. W. Crook, “The reflection and transmission of light by any system of parallel isotropic films,” J. Opt. Soc. Am. **38**, 954–963 (1948). [CrossRef] [PubMed]

**7. **D. E. Aspnes, J. B. Theeten, and F. Hottier, “Investigation of effective-medium models of microscopic surface roughness by spectroscopic ellipsometry,” Phys. Rev. B **20**, 3292–3302 (1979). [CrossRef]

**8. **J. Szczyrbowski, K. Schmalzbauer, and H. Hoffmann, “Optical properties of rough thin films,” Thin Solid Films **130**, 57–73 (1985). [CrossRef]

**9. **L. Névot, B. Pardo, and J. Corno, “Characterization of X-UV multilayers by grazing incidence X-ray reflectometry,” Rev. Phys. Appl. **23**, 1675–1686 (1988). [CrossRef]

**10. **H. E. Bennett, “Specular reflectance of aluminized ground glass and the height distribution of surface irregularities,” J. Opt. Soc. Am. **53**, 1389–1394 (1963). [CrossRef]

**11. **I. Ohlídal, K. Navrátil, and F. Lukeš, “Reflection of light on a system of non-absorbing isotropic film–non-absorbing isotropic substrate with rough boundaries,” Opt. Commun. **3**, 40–44 (1971). [CrossRef]

**12. **I. Ohlídal, K. Navrátil, and F. Lukeš, “Reflection of light by a system of nonabsorbing isotropic film–nonabsorbing isotropic substrate with randomly rough boundaries,” J. Opt. Soc. Am. **61**, 1630–1639 (1971). [CrossRef]

**13. **J. M. Eastman and P. W. Baumeist, “Measurement of microtopography of optical surfaces using a scanning Fizeau interferometer,” J. Opt. Soc. Am. **64**, 1369 (1974).

**14. **J. M. Eastman, “Scattering in all-dielectric multilayer bandpass filters and mirrors for lasers,” in *Physics of Thin Films*, G. Hass and H. M. Francombe, eds. (Academic, 1978), Vol. 10, p. 167.

**15. **C. K. Carniglia, “Scalar scattering theory for multilayer optical coatings,” Opt. Eng. **18**, 104–115 (1979). [CrossRef]

**16. **J. Bauer, L. Biste, and D. Bolze, “Optical-properties of aluminum nitride prepared by chemical and plasmachemical vapor-deposition,” Phys. Status Solidi A **39**, 173–181 (1977). [CrossRef]

**17. **J. M. Zavislan, “Angular scattering from optical interference coatings: scalar scattering predictions and measurements,” Appl. Opt. **30**, 2224–2244 (1991). [CrossRef] [PubMed]

**18. **C. Amra, J. H. Apfel, and E. Pelltier, “Role of interface correlation in light scattering by a multilayer,” Appl. Opt. **31**, 3134–3151 (1992). [CrossRef] [PubMed]

**19. **C. Amra, C. Grèzes-Besset, and L. Bruel, “Comparison of surface and bulk scattering in optical multilayers,” Appl. Opt. **32**, 5492–5503 (1993). [CrossRef] [PubMed]

**20. **C. Amra, “Light scattering from multilayer optics I. tools of investigation,” J. Opt. Soc. Am. A **11**, 197–210 (1994). [CrossRef]

**21. **A. Duparré, J. Ferre-Borrull, S. Gliech, G. Notni, J. Steinert, and J. M. Bennett, “Surface characterization techniques for determining the root-mean-square roughness and power spectral densities of optical components,” Appl. Opt. **41**, 154–171 (2002). [CrossRef] [PubMed]

**22. **J. Meunier, “Optical reflectivity of thin rough films: Application to ellipsometric measurements of liquid films,” Phys. Rev. E **75**, 061601 (2007). [CrossRef]

**23. **A. A. Maradudin, ed., *Light Scattering and Nanoscale Surface Roughness* (Springer, 2010).

**24. **S. Schröder, T. Herffurth, H. Blaschke, and A. Duparré, “Angle-resolved scattering: an effective method for characterizing thin-film coatings,” Appl. Opt. **50**, C164–C171 (2011). [CrossRef] [PubMed]

**25. **J. E. Harvey, S. Schröder, N. Choi, and A. Duparré, “Total integrated scatter from surfaces with arbitrary roughness, correlation widths, and incident angles,” Opt. Eng. **51**, 013402 (2002). [CrossRef]

**26. **T. Herffurth, S. Schröder, M. Trost, and A. D. A. Tünnermann, “Comprehensive nanostructure and defect analysis using a simple 3D light-scatter sensor,” Appl. Opt. **52**, 3279–3287 (2013). [CrossRef] [PubMed]

**27. **J. W. S. Rayleigh, *Theory of Sound* (Macmillan, 1877), Vol. 2.

**28. **S. O. Rice, “Reflection of electromagnetic waves from slightly rough surfaces,” Commun. Pure Appl. Math. **4**, 351–378 (1951). [CrossRef]

**29. **G. R. Valenzuela, “Depolarization of EM waves by slightly rough surfaces,” IEEE Trans. Antennas Propag. **15**, 552–557 (1967). [CrossRef]

**30. **R. Schiffer, “Reflectivity of a slightly rough surface,” Appl. Opt. **26**, 704–712 (1987). [CrossRef] [PubMed]

**31. **V. Freilikher, M. Pustilnik, I. Yurkevich, and V. I. Tatarskii, “Polarization of light scattered from slightly rough dielectric film,” Opt. Lett. **19**, 1382–1384 (1994). [CrossRef] [PubMed]

**32. **D. Franta and I. Ohlídal, “Ellipsometric parameters and reflectances of thin films with slightly rough boundaries,” J. Mod. Opt. **45**, 903–934 (1998). [CrossRef]

**33. **A. Krywonos, “Predicting surface scatter using a linear systems formulation of non-paraxial scalar diffraction,” Ph.D. thesis, University of Central Florida, Orlando (2006).

**34. **S. Schröder, A. Duparré, L. Coriand, A. Tünnermann, D. H. Penalver, and J. E. Harvey, “Modeling of light scattering in different regimes of surface roughness,” Opt. Eng. **19**, 9820–9835 (2011).

**35. **I. Ohlídal and F. Vižd’a, “Optical quantities of multilayer systems with correlated randomly rough boundaries,” J. Mod. Opt. **46**, 2043–2062 (1999). [CrossRef]

**36. **M. Šiler, I. Ohlídal, D. Franta, A. Montaigne-Ramil, A. Bonanni, D. Stifter, and H. Sitter, “Optical characterization of double layers containing epitaxial ZnSe and ZnTe films,” J. Mod. Opt. **52**, 583–602 (2005). [CrossRef]

**37. **I. Ohlídal and F. Lukeš, “Ellipsometric parameters of rough surfaces and of a system substrate-thin film with rough boundaries,” Opt. Acta **19**, 817–843 (1972). [CrossRef]

**38. **I. Ohlídal, F. Lukeš, and K. Navrátil, “Rough silicon surfaces studied by optical methods,” Surf. Sci. **45**, 91–116 (1974). [CrossRef]

**39. **J. T. Butler, “On the number of propagation paths in multilayer media,” Fibonacci Q. **28**, 334–339 (1990).

**40. **L. D. Killough, “The on-line encyclopedia of integer sequences,” http://oeis.org/A011782. Sequence A011782.

**41. **N. J. A. Sloane, “The on-line encyclopedia of integer sequences,” http://oeis.org/A001519. Sequence A001519.

**42. **C. Mallows, N. J. A. Sloane, S. Plouffe, and R. G. Wilson, “The on-line encyclopedia of integer sequences,” http://oeis.org/A007051. Sequence A007051.

**43. **H. Bottomley, “The on-line encyclopedia of integer sequences,” http://oeis.org/A080937. Sequence A080937.

**44. **N. J. A. Sloane, “The on-line encyclopedia of integer sequences,” http://oeis.org/A024175. Sequence A024175.

**45. **N. J. A. Sloane, “The on-line encyclopedia of integer sequences,” http://oeis.org/A000108. Sequence A000108.

**46. **M. Kildemo, O. Hunderi, and B. Drévillon, “Approximation of reflection coefficients for rapid real-time calculation of inhomogeneous films,” J. Opt. Soc. Am. A **14**, 931–939 (1997). [CrossRef]

**47. **I. Ohlídal, “Reflectance of multilayer systems with randomly rough boundaries,” Opt. Commun. **71**, 323–326 (1989). [CrossRef]

**48. **I. Ohlídal, “Approximate formulas for the reflectances, transmittances, and scattering losses of nonabsorbing multilayers systems with randomly rough boundaries,” J. Opt. Soc. Am. A **10**, 158–170 (1993). [CrossRef]

**49. **I. Ohlídal, K. Navrátil, and M. Ohlídal, “Scattering of light from multilayer systems with rough boundaries,” in *Progress in Optics*, E. Wolf, ed. (Elsevier, 1995), Vol. 34, pp. 249–331. [CrossRef]

**50. **A. H. Stroud and D. Secrest, *Gaussian Quadrature Formulas* (Prentice-Hall, 1966).

**51. **H. R. Philipp, “Silicon dioxide (SiO_{2}) (glass),” in *Handbook of Optical Constants of Solids*, E. Palik, ed. (Academic, 1985), Vol. I, pp. 749–763. [CrossRef]

**52. **H. R. Philipp, “Silicon nitride (Si_{3}N_{4}) (noncrystalline),” in *Handbook of Optical Constants of Solids*, E. Palik, ed. (Academic, 1985), Vol. I, pp. 771–774. [CrossRef]

**53. **C. M. Herzinger, B. Johs, W. A. McGahan, J. A. Woollam, and W. Paulson, “Ellipsometric determination of optical constants for silicon and thermally grown silicon dioxide via a multi-sample, multi-wavelength, multi-angle investigation,” J. Appl. Phys. **83**, 3323–3336 (1998). [CrossRef]