4.12. Inversion method to derive complex-valued refractive index functions of transparent materials

4.12.1. Background and theory

As stated above, it is essential to have knowledge of the complex-valued refractive index (cRI) functions of all materials used in physically-based raytracing. While it is possible to retrieve models or empirical data for some “bulk” materials, e.g. of the refractiveIndex database (Polyanskiy, 2022), this data is unfortunately currently not available for many other materials. For building science applications, it is of particular importance to have models for coated and uncoated glass panes, to be able to model glazings accurately. While the necessary refractive index information is not available, Lawrence Berkley National Laboratories kindly maintains the freely accessible and continuously updated glazing database IGDB (LBNL Windows and Daylighting Group, 2021), containing spectral measurements for almost all relevant types of glass. The database provides transmittance and reflectance spectra for the global radiation range. The datasets typically include measured spectral functions for the frontside reflectance 𝑅𝑓, backside reflectance 𝑅𝑏 and transmittance 𝑇. These measurements are generally provided for normal, i.e. near-normal, incidence only. The data is provided for uncoated as well as for coated glass panes.
A method has been developed to generate the required cRI functions based on these measured spectral data by inversion, i.e. appropriate parametric models for the cRI functions are varied until the modelled reflectance and transmittance spectra match empirical spectra provided by IGDB. The two real-valued components forming the complex-valued refractive index function, the refractive index 𝑛 and the extinction coefficient , have to be determined for the entire global radiation spectrum. Due to the many resulting degrees of freedom, the inversion can only be performed using a multi-parameter optimisation algorithm.

slab model is applied to efficiently model a coated or uncoated glass pane’s front- and backsidereflectances and transmittance. As discussed in section 4.3.7, this model assumes a non-coherentsuperposition of many interreflections within the glass pane. Figure 53 provides a schematic overviewof the model. Since the measurements are performed with air contact, the media 0 and 2 are modelledas air (approximated as vacuum). Medium 1 is the actual glass material. Interfaces 01 and 12 aremodelled using the generalized Fresnel approach (see section 4.7). At interface 01, a thin-film stack consisting of N layers can additionally be considered. The interfaces’ reflectivity and transmissivityvalues are denoted with ordered indices to indicate the direction, even though most are identical forboth directions. All coefficients are determined using equations (53). Note that for the case of a thinfilm stack present at interface 01, any absorption within the thin-film will also be included in themodel, as for this case 𝑇01 + 𝑅01 < 1 will apply.
Beyond the Fresnel coefficients, relevant for the scattering processes at the interfaces, the transmittance 𝑇1 of the glass pane is required. It is obtained by applying Lambert-Beer’s law (see section 4.3.6):

Whereas the attenuation coefficient  is calculated using equation (40), and the angle 𝜃1 is calculated based on Snell’s law (46).
The derivation of the total reflectance and transmittance, being the superposition of all interreflections, can easily be comprehended if the coefficients R and T are understood as probabilities. In this view, the probabilities for each path have to be added, while the required probabilities forming a path have to be multiplied. A similar derivation, but with slightly different components, is, e.g., provided by Stenzel (2016, p. 131).
The infinite summation of all reflection paths (see Figure 53) then leads to a geometric series, which can easily be evaluated into a simple fractional term:

In complete analogy, the total reflectance for the front side and, by flipping the direction and exchanging the interface indices correspondingly, the backside of the glass pane can be derived. The frontside reflectance 𝑅𝑓,𝑝𝑎𝑛𝑒 and the backside reflectance 𝑅𝑏,𝑝𝑎𝑛𝑒 are:

Note that the expressions could further be simplified by expressing transmittance values by reflectivities and applying identities for transmittance values in reverse directions. However, since the coefficients are readily provided by the algorithms and in order to average out numerical inaccuracies, the expressions were implemented as stated in equations (76) and (77). Note also that all quantities Light, sun and optics – applied principles, models and methods RadiCal, D. Rüdisser 95 are functions of wavelength and incidence angle; for reasons of clarity, the corresponding arguments were omitted in the above equations.

4.12.1.2. A parametric model for coated and uncoated glass panes

After some testing, the following models were identified as best suited to determine the cRI functions by inversion. As is always the case for such optimisation problems, it is critical to provide a physically appropriate model with as many variable parameters as required – but only as few as necessary, to avoid overtraining the model. In the implemented method, mainly natural cubic spline functions (G., Press et al., 1991) with N knots are used to individually approximate the real-valued functions 𝑛 and that constitute the complex-valued refractive index function. Since often ranges over several orders of magnitude, a logarithmic mapping of the spline was tested but was found not to perform better than the linear approach.
The following constraints were applied in the optimisation process:

• the wavelength (“x-coordinate”) of the first knot of the spline was constrained to range within 0 nm and 300 nm.
• the wavelength (“x-coordinate”) of the last knot of the spline was constrained to range within 2500 nm and 3000 nm.
• all values for glazings were constrained in the range of [1E-8, 10].
• all values for coatings were constrained in the range of [1E-8, 100].
• all 𝑛 values for glazings were constrained in the range of [1E-5, 4].
• all 𝑛 values for coatings were constrained in the range of [1E-5, 6].

Note that the actual 𝑛 and functions for the coatings do not have to reflect the true physical properties if taken individually and, therefore, can seem entirely out of range. A two-layer substitute model is applied to approximate the behaviour of the thin-film stack, potentially consisting of a double-digit number of different materials, e.g. for a low-E coating (Ding and Clavero, 2017). The actual optical behaviour of the thin-film assembly is determined as the result of the multiplications of the individual characteristic matrices of each layer (see equations (50) – (52)). This matrix multiplication relies on the cRI functions of the individual layers, as well as on their thicknesses 𝑑𝑖 . Therefore, physically not meaningful values for 𝑛, and d achieved for the single layers of the twolayer model can still be suitable to resemble the correct physical behaviour of a multi-layer assembly. While spline functions are used to model the refractive index and extinction coefficient functions of all coating layers, the refractive index of the glass can be modelled based on a less complex model. The following dispersion formula (see TRIsimpleDispersionAndKspline of Table 5) has proven to be suitable to model the relevant uncoated glass pane types:

As can be seen, this function relies on three parameters only, while the natural cubic spline function has N degrees of freedom.
Since IGDB contains data for coated and uncoated glass panes, the uncoated panes’ spectra are always used to determine the model parameters for the plain glass material. The given thickness of the glass pane is provided as a constant model parameter.

In order to determine the model for coated glass panes, the model parameters of the uncoated glass, comprising three constants for the refractive index (eq. (79)) and a series of N spline knots for the extinction coefficient, are used as constant values. Only the parameters relevant for the two coating layers, comprising the thickness and the two series of knots coordinates for each coating layer, are the free parameters in the optimisation process to replicate the measured spectra of the coated glass pane.

4.12.1.3. An ideal optimisation cost function

Ideally, the optimisation should be controlled by applying transmittance and reflectance spectra for multiple incidence angles as the target quantity of the optimisation. Based on the experience gained when working with the optimisations, it seems likely that a number of only two to four different incident angles should be sufficient for this purpose.
The objective of the optimisation process is defined by providing a cost function that rates the optimality of a solution. For the present case, the following definitions are proposed: first, the deviations (∆𝑇, ∆𝑅𝑓 and ∆𝑅𝑏) of the model vs. the measurement for a specific incidence angle 𝜃𝑛 and the wavelength are defined as:

Where 𝑇, 𝑅𝑓 and 𝑅𝑏 are the model results, according to (76) – (78), and the corresponding values marked with the hat accent represent measurement values provided by the IGDB spectrum files. The global radiation spectrum is sampled in two dimensions. First, by using 𝑀 distinct, energetically equidistant wavelengths (see section 4.4.2), and second, for 𝑁 available incidence angles 𝜃𝑛. The loss functions to assess the quality of the fitting are defined as:

The weighted loss functions finally determine the cost function for the optimisation process:

The coefficients 𝒄𝒊 can be used to balance the three components of the cost function. Unfortunately, this approach could not be applied and tested fully, as only spectra for one angle (near-normal incidence) were available. However, based on the presented definitions, an extended cost function is introduced below.

4.12.1.4. The pragmatic optimisation cost function

It turned out that applying spectral data only for one angle (near-normal incidence) was insufficient to perform a meaningful determination of the model parameters for the refractive index functions. As the model offers a relatively large number of free parameters, overfitting occurs if insufficient fitting data is provided. Therefore, additional data or constraints are needed. Since no measurements for varying angles are available, two additional constraints are imposed:

• 1. the normal-incidence reflectance and transmittance spectra should also be valid for a specified low, but non-zero, incidence angle 𝜃1.
• 2. the ratio of the total transmittance for a specified higher oblique angle 𝜃2 to the total transmittance at normal incidence should match a specified value.

Both constraints are physically reasonable. The first condition helps to avoid spectrally sensitive thin-film solutions that frequently occur if the constraint is not applied. Such solutions will provide the correct reflectance behaviour for the normal incidence but will show strong spectral and angular dependences for oblique angles, e.g. as the colourful “rainbow” interference known from thin oil or soap films. Such effects are, in fact, avoided in the production of glazings, as this would impair the visual quality of the glazing. Since the angular characteristics of glazing generally show very little change at incidence angles less than 10 to 20 degrees (Rubin et al., 1999; Roos et al., 2001), an angle 𝜃1 within this range can be chosen.
To further restrict the optimisation results to physically plausible solutions, a target value for the decrease of the transmittance with higher angles is specified in the second constraint. While the entire spectrally resolved reflectance and transmittance data is provided in the first constraint, here, only the total (=spectrally integrated) value of the transmittance is considered. Since the sampling wavelengths 𝑚 are provided for energetically equidistant intervals, the total transmittance can simply be calculated by averaging the values. Hence, the additional loss function can be written as:

Where 𝑓𝑇 represents the target value for the decay of the total transmittance. The other loss functions remain the same as stated above, however two angles are considered there now: 𝜃0 representing near-normal incidence, and 𝜃1 as low angle for which the normal-incidence spectra should still be valid.

In the practical implementation, the relevant angles and weighting-factors were defined as: 𝜃0 = 10−3 ° (nonzero for numerical stability), 𝜃1 = 15 ° , 𝜃2 = 70 ° and 𝑐1 = 𝑐2 = 𝑐4 = 1, while 𝑐3 = 0.5.

4.12.2. Implementation

The implementation is contained in the module rc_fittrrspectrum3, whereby the version number indicates the evolution of the module. The class TfitTRRSpectrum contains all functions necessary for the optimisation process and additional tools used to import, export and report the optimisation results.
The progress of the potentially multi-hour-long optimisation can be monitored in text form, as well as in the form of an updating chart.

The implementation is an application of the simulated annealing algorithm (see section 5.2), where equation (84) is used as a cost function. The relations for the thick slab model derived above, see (76) – (78), are used to compute the analytic transmittance and reflectance values used for the fitting to the empirical data during the optimisation process.
In the implementation, the ipfisn and ipfspn of the NumLib library (Free Pascal Team, 2021) are used to solve and evaluate the linear equation systems for the applied natural cubic spline fitting.

4.12.3. Testing and validation

In order to demonstrate the potential of the method, the complex-valued refractive index (cRI) functions of three different materials are determined from their transmission and reflection spectra. The objective of the test is two-fold:

• 1.) If applied in the raytracing (or thick-slab model), the determined cRI function should, with good accuracy, result in transmission and reflection spectra matching the measured data they were generated on. The match is quantified with three different parameters that are also applied during the optimisation process:
  • a. Mean error of samples: as the wavelength samples are distributed according to the global radiation spectral density, the total mean error also coincides with the integrated total error considering the entire spectrum.
  • b. Root-mean-square error of samples: this value is the best indicator for the average spectral deviation of the three spectra.
  • c. Max. error of samples: the maximum deviation for any sample in any of the three spectra (transmission, reflection-frontside, reflection-backside) is determined.
• 2.) The determined cRI functions should be physically correct, i.e. represent the “true” cRI functions.

4.12.3.1. Case 1 – determination of cRI functions for uncoated window pane

In this first case, the cRI of an uncoated pane used for insulated glazing units (IGU) is determined. The result is later used in the full-system validation (chapter 8). The pane is of type Planibel Clearlite by the manufacturer AGC Glass Europe. The nominal thickness is 4 mm, while the actual thickness is 3,85 mm. The related spectra were exported from the IGDB database (LBNL Windows and Daylighting Group, 2021).
For the inversion process to determine the cRI of the uncoated glass, the dispersion relation stated in equation (79) was used to model the refractive index function (𝑛), while a spline function with 10 knots was used to model the extinction coefficient function. The knots are moved in both dimensions during the optimisation process. Only the wavelength coordinates of the first and last knot are held constant at the boundaries of the global radiation spectrum. The optimisation process comprised 35 million iteration steps and lasted a few hours. The finally determined cRI is depicted in Figure 55. The dots on the blue line indicate the knots of the spline function.

Figure 56 shows the cRI as determined by Rubin (1985) for unspecified “clear glass”. As Rubin notes in his publication, the absorptance and reflectance of the window glass are affected by the FeO concentration and further additives. This means that significant variations can be expected among different clear glasses. However, as can be seen by comparison of Figure 55 and Figure 56, the match between the values for a modern glass, determined with the novel method presented here and the values determined by Rubin in 1985, is remarkably high.
The match of the transmission and reflection spectra is depicted in Figure 57. The precision of the reconstruction is very high, as the only visible deviation occurs in the transmission spectrum at approx. 2.2 µm. However, the deviation is still small, and the energetic impact of the deviation is even less significant, as this part of the global radiation spectrum exhibits a low spectral flux density.
The maximum sample error is 1.1% for the transmittance and less than 0.08% for the reflectance, and the root mean squared error is 0,4% for all transmittance and reflectance samples. The mean sample error is 0.09% for the transmittance and + 0.001% for the reflectance samples. Since mean sample errors are spectrally weighted, the errors are consistent with the total error for the entire global radiation spectrum.

Figure 58 finally shows the angular dependence of the transmittance and reflectance of the glass pane for global radiation. As discussed in the state-of-the-art analysis (section 3.2), even the more advanced, commonly used modelling approaches apply simplified models in order to determine the angular dependence of transmittance and reflectance. The depicted angular dependencies are relevant for global radiation determined by spectrally weighted integration.

4.12.3.2. Case 2 – determination of cRI functions for low-E coated window pane

In this example, a low emissivity (low-E) coating (of type iplus Top 1.1) is applied to the previously considered uncoated glass pane. As discussed in the theory section, the coating, consisting of multiple layers, is approximated using a two-layer substitute model. The already determined cRI for the 3.85 mm thick glass pane is held constant in the optimisation process. The optimisation involves the 𝑛 and functions of both layers, where each function is modelled using spline functions with 10 knots. Beyond that, the two thicknesses d of the layers are variables during the optimisation process. As discussed above, an additional parameter for the total angular weakening of transmittance is required due to the lack of angular spectral measurements. The value was taken from measurement data contained in a publication by Rubin et al. (1999). The transmission of the specimen ‘double Ag 3mm clear’ was found to decrease to 63.8% of its near-normal value at 70°. This value was used as a target for the optimisation process. The optimisation terminated after 12.3 million iterations. The layer thicknesses for the substation model were determined as 9.0 nm for the first (=glass-sided) layer and 37.6 nm for the second (=top) layer, respectively.

be assumed that it comprises several different materials, including one to three silver layers. Due to production and durability aspects, the actual layering is usually complex and involves additional materials (Ding et al., 2017).
The comparison of the measured vs. the modelled spectra (Figure 61) shows that the applied twolayer model is able to replicate the transmittance and reflectance behaviour of the more complex layering very well. This is also supported by the error measures: transmission: max. sample error: 2.8%, mean sample error: <0.01%; reflection-frontside: max. sample error: 0.85%; mean sample-error: <0.01%; reflection-backside: max. sample error: 3.2%; mean sample-error: <0.01%;).

Figure 62 shows the angular dependencies for transmittance, frontside- and back-side reflectance as well as absorptance. Additionally, two references for the angular transmittance are depicted. The first is taken from Rubin et al. (1999) and represents measured data for a silver double-coated pane with a nominal thickness of 3 mm. Since the transmission ratio at 70° was an optimisation target, it shows a perfect match at this angle. However, a good match is also visible regarding the general shape of the curve. The second reference was taken from Roos et al. (2001). Figure 10 of the cited publication provides a model fit for a silver-coated Planitherm pane with unspecified thickness and few angular measurements. The angular dependence modelled by Roos et al. shows slightly higher values than are determined here.

4.12.3.3. Case 3 – determination of cRI functions for ETFE membrane

In order to include a non-glass example of practical relevance, the cRI function of a clear ETFE membrane, as used for skylights in buildings, is determined. Transmission and reflection spectra for the relevant global radiation range were extracted from a chart provided by a manufacturer (Buitink Technology, 2022) (“standard foil” curves). The thickness of the sample that was used to determine the spectra is specified as 200 µm.
While an appropriate dispersion relation was used to model the 𝑛 function for the uncoated glass, the cRI of the ETFE membrane is now determined by optimising two spline functions. The number of available knots was set to 9 for the 𝑛-component and 19 for the -component, as more fluctuation in the absorption can be expected. Even though the optimisation is offered a very high number of degrees of freedom for the spline functions, remarkably smooth and plausible 𝑛 and  functions are determined in the optimisation process (see Figure 63). As can be seen in the 𝑛 function, the optimisation process concentrated most knots below 1 µm to capture the rise towards shorter wavelengths accurately.

Again literature was reviewed to assess the validity of the inversion result. A reference containing measurement data has been found in a publication on Optical properties of materials for concentrator photovoltaic systems (French et al., 2009). The publication includes results of spectroscopic ellipsometer and spectrophotometer measurements performed on a clear ETFE membrane.

The comparison of the measured 𝑛-function (blue curve in Figure 64 left) vs. the function determined by inversion (Figure 63) shows that for wavelengths longer than 1,3 µm, both functions quickly converge to an almost constant value around 1.40. In the UV range, a sharp rise towards shorter wavelengths, typical for polymers, is visible. While the shape of the rise is very similar, a quantitative deviation is clearly visible in this region. However, the products are of different manufacturers, and limitations of the measurement devices in the UV region can distort the measurement results at the end of the measurement range.
The limited vertical resolution of the absorption coefficient chart (blue curve in Figure 64 right) allows only a crude assessment of the validity of the -function, containing the extinction coefficients. As can be seen in the chart, three main absorption regions can be identified (indicated with blue arrows). These correspond to three absorption regions determined in the inversion, see Figure 63. Note that Figure 64 uses a linear scale, while the extinction coefficient in Figure 63 is depicted using a logarithmic scale. By converting the estimated absorption values using Beer-Lamberts law, eq. Figure 64 Measured index of refraction and absorption coefficients of ETFE membranes (blue lines) source: French et al. (2009) Light, sun and optics – applied principles, models and methods RadiCal, D. Rüdisser 104 (26), it can be shown that the values match, at least in order of magnitude, being the relevant scope here.
Figure 64 Measured index of refraction and absorption coefficients of ETFE membranes (blue lines) source: French et al. (2009) Light, sun and optics – applied principles, models and methods RadiCal, D. Rüdisser 104 (26), it can be shown that the values match, at least in order of magnitude, being the relevant scope here. Figure 65 demonstrates that the match between the measured and modelled transmission and absorption spectra (the optimisation process’s target) is excellent: transmission: max. sample error: 1.8%, mean sample error: <0.01%; reflection: max. sample error: 0.4%; mean sample-error: <0.01%. The derived angular transmittance and reflectance for global radiation are depicted in Figure 66.