Color Simulation of Multilayered Thin Films Using Python

Abstract

Physical insight into a material can be first gained by its color, as the reflectance spectrum of an object reflects its microstructure and complex refractive indices. Here, we present a comprehensive overview of electrodynamics and optics related to reflectance spectra and color. We provide an open-source Python code for simulating reflectance spectra and extracting color values. The validity and applicability of the code are demonstrated through a comparative analysis with both the literature and experimental data. For SnO₂ and ZnO thin films deposited on SiO₂/Si substrates using rf sputtering, the Python code and simulation predict color variations with the film thickness and effectively capture their angular dependence. This code will help in understanding and making use of color-related phenomena. It can be further used and developed for various purposes, particularly machine learning, which requires extensive spectral and color data for model training.

1. Introduction

The color of an object is one of the first characteristics that captures our attention. In the 17th and 18th centuries—the Age of Enlightenment—vigorous efforts were made to understand natural phenomena, and the study of light and color was systematized by luminaries such as Isaac Newton and Augustin-Jean Fresnel. With the advancement of computational techniques and theoretical frameworks, the simulation of color based on electrodynamics and optics has become an effective tool for predicting and elucidating experimental results in condensed matter physics and materials science [1,2,3]. As thin-film interference and multilayer coatings have recently garnered significant attention due to their ability to manipulate light at the nanoscale, color simulation techniques have been employed increasingly frequently. These techniques aim to describe or predict colors based on various factors, such as surface oxidation [4], crystallinity [5], and composition [6]. They have also been applied to explain the color variations generated in Mie resonators based on liquid crystals [7].

Despite the availability of commercial software and web-based tools used for optical simulations [8,9,10], access to internal processes is often restricted due to proprietary algorithms, limiting customization. Nevertheless, researchers and engineers frequently look for flexible and customizable tools to adjust their simulations to meet specific experimental conditions and material properties.

Recently, many individuals have used programming languages to carry out their own tasks, such as data analysis [11,12], fitting [13,14], and simulation [15,16]. Among the different programming languages, Python stands out because of its open-source license, user-friendly interface, and various libraries which include functions essential for simulation, including calculation, visualization, and interpolation [17,18,19]. In addition to these features, advanced algorithms such as machine learning algorithms have been developed and utilized [20,21,22]. Consequently, there is a growing demand for fundamental theoretical simulations that can complement and work with modern programming languages and algorithmic approaches.

Here, we introduce a Python-based approach for simulating the color of thin-film multilayers, while providing ample background on its underlying physics and coding—from basic electrodynamics, including the Fresnel equations, to the conversion of reflectance spectrum to RGB values based on a color matching function. Our code requires only information on the complex refractive indices (within the visible range) and thickness of each layer and a specific incident angle. To validate this code, we compare the simulated colors of several thin films with other results: SiO₂ films from the existing literature and SnO₂ and ZnO thin films with a thickness gradient, which were prepared for this study. The code is straightforward and allows individuals with basic Python knowledge to easily modify or add conditions, facilitating color simulation under various conditions and for integration with other applications. We hope that this code helps students and researchers working on materials science and AI-based research in various ways.

2. Materials and Methods

2.1. Theoretical Background

2.1.1. Reflection from a Bulk Surface

When light is incident on an object, it undergoes multiple absorptions, reflections, and transmissions at the surface, and the color of that object is determined by the spectrum of light reflected from its surface. Fresnel equations describe the reflected and transmitted amplitudes of light with respect to that of the incident wave, which can be derived from the boundary conditions of electric and magnetic fields at an interface.

We define the plane of incidence that contains the surface normal vector and the propagation vector of incident light and determine p- or s-polarization from the direction of E oscillation with respect to the plane of incidence (Figure 1). The p-polarization indicates that E oscillates parallel to the plane of incidence (Figure 1a). Here, we first derive the Fresnel equations for the case of polarization parallel to the plane of incidence (i.e., p-polarization). Since the parallel component of E and the normal component of B with respect to the interface need to be continuous, we arrive at the following equations:

(1)$E_{\mathrm{i}}\cos {\mathit{\theta }}_{\mathrm{i}}-E_{\mathrm{r}}\cos {\mathit{\theta }}_{\mathrm{r}}=E_{\mathrm{t}}\cos {\mathit{\theta }}_{\mathrm{t}}$

(2)${\mathit{B}}_{\mathrm{i}}+{ \mathit{B}}_{\mathrm{r}}={ \mathit{B}}_{\mathrm{t}}$

where the subscripts i, r, and t represent the incident, reflected, and transmitted waves, respectively. Since $\mathit{B}=E·n/c$ (n: refractive index,$n \equiv c/v$) where $v$ is the speed of an electromagnetic wave in matter [23], Equation (2) indicates that

(3)$n_{\mathrm{i}}\left(E_{\mathrm{i} }+E_{\mathrm{r}}\right)=n_{\mathrm{t}}{E}_{\mathrm{t}}$

where $n_{\mathrm{i}}$ and $n_{\mathrm{t}}$ are the refractive indices of the media before and after transmission. Given the laws of reflection (${\mathit{\theta }}_{\mathrm{r}}={\mathit{\theta }}_{\mathrm{i}}$) and refraction (i.e., Snell’s law; $n_{\mathrm{i}}\sin {\mathit{\theta }}_{\mathrm{i}}=n_{\mathrm{t}}\sin {\mathit{\theta }}_{\mathrm{t}}$) and by solving Equations (1) and (2), we obtain the reflection and transmission coefficients for p-polarized light:

(4)$r_{\mathrm{p}} \equiv {\displaystyle \frac{E_{\mathrm{r}}}{E_{\mathrm{i}}}}={\displaystyle \frac{n_{\mathrm{t}}\cos {\mathit{\theta }}_{\mathrm{i}}-n_{\mathrm{i}}\cos {\mathit{\theta }}_{\mathrm{t}}}{n_{\mathrm{t}}\cos {\mathit{\theta }}_{\mathrm{i}}+{ n}_{\mathrm{i}}\cos {\mathit{\theta }}_{\mathrm{t}}}}$

(5)$t_{\mathrm{p}} \equiv {\displaystyle \frac{E_{\mathrm{t}}}{E_{\mathrm{i}}}}={\displaystyle \frac{2n_{\mathrm{i}}\cos {\mathit{\theta }}_{\mathrm{i}}}{n_{\mathrm{t}}\cos {\mathit{\theta }}_{\mathrm{i}}+n_{\mathrm{i}}\cos {\mathit{\theta }}_{\mathrm{t}}}}.$

For s-polarization (Figure 1b), the boundary conditions for E and B become

(6)$E_{\mathrm{i}}+E_{\mathrm{r}}=E_{\mathrm{t}}$

(7)${-\mathit{B}}_{\mathrm{i}}\cos {\mathit{\theta }}_{\mathrm{i}}-{\mathit{B}}_{\mathrm{r}}\cos {\mathit{\theta }}_{\mathrm{r}}={-\mathit{B}}_{\mathrm{t}}\cos {\mathit{\theta }}_{\mathrm{t}}.$

The reflection and transmission coefficients for s-polarized light can be obtained by the same procedure.

(8)$r_{\mathrm{s}} \equiv {\displaystyle \frac{E_{\mathrm{r}}}{E_{\mathrm{i}}}}={\displaystyle \frac{n_{\mathrm{i}}\cos {\mathit{\theta }}_{\mathrm{i}}-n_{\mathrm{t}}\cos {\mathit{\theta }}_{\mathrm{t}}}{n_{\mathrm{i}}\cos {\mathit{\theta }}_{\mathrm{i}}+n_{\mathrm{t}}\cos {\mathit{\theta }}_{\mathrm{t}}}}$

(9)$t_{\mathrm{s}} \equiv {\displaystyle \frac{E_{\mathrm{t}}}{E_{\mathrm{i}}}}={\displaystyle \frac{2n_{\mathrm{i}}\cos {\mathit{\theta }}_{\mathrm{i}}}{n_{\mathrm{i}}\cos {\mathit{\theta }}_{\mathrm{i}}+n_{\mathrm{t}}\cos {\mathit{\theta }}_{\mathrm{t}}}}.$

The reflectance (R) and transmittance (T) are defined by the ratio of the reflected and transmitted intensity (average power per unit area) to the incident intensity ($I=1/2ϵv{E}^2$); thus,

(10)$R \equiv {\displaystyle \frac{I_{\mathrm{r}}}{I_{\mathrm{i}}}}={\left|{\displaystyle \frac{E_{\mathrm{r}}}{E_{\mathrm{i}}}}\right|}^2={ \left|r\right|}^2$

(11)$T \equiv {\displaystyle \frac{I_t\cos {\mathit{\theta }}_{\mathrm{t}}}{I_{\mathrm{i}}\cos {\mathit{\theta }}_{\mathrm{i}}}}=\left({\displaystyle \frac{n_{\mathrm{t}}\cos {\mathit{\theta }}_{\mathrm{t}}}{n_{\mathrm{i}}\cos {\mathit{\theta }}_{\mathrm{i}}}}\right){\left|{\displaystyle \frac{E_{\mathrm{t}}}{E_{\mathrm{i}}}}\right|}^2=\left({\displaystyle \frac{n_{\mathrm{t}}\cos {\mathit{\theta }}_{\mathrm{t}}}{n_{\mathrm{i}}\cos {\mathit{\theta }}_{\mathrm{i}}}}\right){\left|t\right|}^2.$

The reflectance depends on the polarization of the light, and we consider natural, or unpolarized, light in this article. The reflectance of unpolarized light (R_n) can be given by the averaged reflectance of all incident waves with different polarization angles:

(12)$R_{\mathrm{n} }={\displaystyle \frac{1}{2\pi }}{\int }_{-\pi }^{\pi }{R}_{\mathsf{\psi }}d\psi ={\displaystyle \frac{\left(R_{\mathrm{p}}+R_{\mathrm{s}}\right)}{2}}$

where $\psi$ is the polarization angle [24].

2.1.2. Reflection from a (Multilayered) Thin Film

For a thin film or multilayered structure, we need to consider the absorption and optical interference of multiple reflections, transmissions, and absorptions occurring at interfaces underneath the surface.

This absorption can be formulated by incorporating a complex refractive index into the wave function (Figure 2). The wave function of a plane wave traveling in the x direction can be written in complex notation as follows:

(13)$\stackrel{\mathrm{~}}{\mathbf{E}}={\stackrel{\mathrm{~}}{\mathbf{E}}}_0\exp [i(\omega t-\tilde{K}x\left)\right]$

where ${\stackrel{\mathrm{~}}{\mathbf{E}}}_0$, $\omega$, t, $\tilde{K}$, and x are the complex wave amplitude’s, angular frequency, time, complex wavenumber, and position, respectively. The complex wavenumber can be written as

(14)$\stackrel{\mathrm{~}}{K} \equiv {\displaystyle \frac{2\pi \stackrel{\mathrm{~}}{N}}{\lambda }}$

where $\stackrel{\mathrm{~}}{N}$ is the complex refractive index. $\stackrel{\mathrm{~}}{N}$ is defined as

(15)$\stackrel{\mathrm{~}}{N} \equiv n-\mathit{ik}$

where n and k are the refractive index and extinction coefficient, respectively. By substituting Equations (14) and (15) into Equation (13), we obtain an expression for the electromagnetic wave that takes attenuation (i.e., exponential decay) due to the absorption into account:

(16)$\stackrel{\mathrm{~}}{\mathbf{E}}={\stackrel{\mathrm{~}}{\mathbf{E}}}_0\exp \left(-{\displaystyle \frac{2\pi \mathit{kx}}{\lambda }}\right)·\exp \left[i\left(\omega t-{\displaystyle \frac{2\pi n}{\lambda }}x\right)\right].$

In addition to absorption, we need to consider optical interference. Figure 3 shows the presence of reflected light from the interfaces between (1) air and the thin film (upper) and (2) the thin film and its substrate (lower). The reflected light waves will interfere either constructively or destructively, depending on the phase difference between the two waves. The path length difference between the two waves is $\stackrel{-}{\mathrm{AB}}+\stackrel{-}{\mathrm{BC}}-\stackrel{-}{\mathrm{AD}}$. The paths $\stackrel{-}{\mathrm{AB}}$ and $\stackrel{-}{\mathrm{BC}}$ traverse a medium with a refractive index of $N_1$, while path $\stackrel{-}{\mathrm{AD}}$ traverses a medium with a refractive index of $N_0$. Therefore, the phase difference ($\alpha$) is obtained by applying the different wavenumbers ($K=2\pi N/\lambda$) as follows:

(17)$\alpha ={\displaystyle \frac{2\pi N_1}{\lambda }}\left(\stackrel{-}{\mathrm{AB}}+\stackrel{-}{\mathrm{BC}}\right)-{\displaystyle \frac{2\pi N_0}{\lambda }}\stackrel{-}{\mathrm{AD}}.$

Considering the incident and transmitted angles, we obtain $\stackrel{-}{\mathrm{AD}}=\stackrel{-}{\mathrm{AC}}\sin {\mathit{\theta }}_{\mathrm{i}}$ and $\stackrel{-}{\mathrm{AC}}=2d\tan {\mathit{\theta }}_{\mathrm{t}}$, which lead to

(18)$\stackrel{-}{\mathrm{AD}}=2d\tan {\mathit{\theta }}_{\mathrm{t}}·\sin {\mathit{\theta }}_{\mathrm{i}}.$

By incorporating Snell’s law ($\sin {\mathit{\theta }}_{\mathrm{i}}=\sin {\mathit{\theta }}_{\mathrm{t}}{N}_1/N_0$),

(19)$\stackrel{-}{\mathrm{AD}}=2d{\displaystyle \frac{{\sin }^2{\mathit{\theta }}_{\mathrm{t}}}{\cos {\mathit{\theta }}_{\mathrm{t}}}}·{\displaystyle \frac{N_1}{N_0}}.$

Using Equation (19) and $\stackrel{-}{\mathrm{AB}}=\stackrel{-}{\mathrm{BC}}=d/\cos {\mathit{\theta }}_{\mathrm{t}}$, we arrive at

(20)$\alpha ={\displaystyle \frac{4\pi d{N}_1}{\lambda }} · \cos {\mathit{\theta }}_{\mathrm{t}}.$

The plane wave, when considering the phase difference, becomes

(21)$\exp \left[i\left\{\omega t-\left(\mathit{Kr}+\alpha \right)\right\}\right]=\exp \left[i\left(\omega t-\mathit{Kr}\right)\right] · \exp \left(-i\alpha \right).$

Equation (21) indicates that the phase difference can be imposed by multiplying $\exp \left(-i\alpha \right)$. To account for the optical interference, we need to find the superposition of the light waves which are reflected from both upper and lower boundaries. The primary wave (i.e., that directly reflected from the film surface) becomes $E_1={r_{01}E}_{\mathrm{i}}$, which is the reflection coefficient determined by the Fresnel equations for the boundary between the air and the thin film. The secondary wave, which is transmitted at the boundary between the air and the thin film, is then reflected from the substrate and transmitted at the boundary between the thin film and air and can be written by imposing the coefficients corresponding to each stage as follows:

(22)$E_2=t_{01}{t}_{10}{r}_{12}\exp \left(-i\alpha \right)E_{\mathrm{i}}.$

As each additional order is considered, the wave undergoes two reflections (r₁₀ and r₁₂) and acquires a phase difference. Each subsequent term is multiplied by a factor of ${r_{10}r}_{12}\exp \left(-i\alpha \right)$. The total reflected wave (E₀₁₂) can be expressed as the sum of the primary wave (E₁), the secondary wave (E₂), the third wave (E₃), and so on:

(23)$\begin{array}{c}{E}_{012}=(r_{01 }+t_{01}{t}_{10}{r}_{12}\exp \left(-i\alpha \right)+{t_{01}{t}_{10}{r}_{10}{r}_{12}^2\exp \left(-i2\alpha \right)+\dots )E}_{\mathrm{i}}\\ =[r_{01 }+t_{01}{t}_{10}{\displaystyle \sum _{\mathrm{m}=1}^{\infty }}{r}_{10}^{m-1}{r}_{12}^{\mathrm{m}}\exp \left(-\mathrm{m}i\alpha \right)]E_{\mathrm{i}}=r_{012}{E}_{\mathrm{i}}.\end{array}$

$r_{012}$ thus represents the reflection coefficient obtained when considering the optical interference. Using the formula for the sum of an infinite geometric series, we can simply write Equation (23) as

(24)$r_{012}=r_{01}+{\displaystyle \frac{t_{01}{t}_{10}{r}_{12}\exp \left(-i\alpha \right)}{1-r_{10}{r}_{12}\exp \left(-i\alpha \right)}}.$

Using $r_{10}=-r_{01}$ and $t_{01}{t}_{10}=1-r_{01}^2$, obtained by Equations (4), (5), (8), and (9), $r_{012}$ can be reduced further [25]:

(25)$r_{012}={\displaystyle \frac{r_{01}+r_{12}\exp \left(-\mathrm{i}\alpha \right)}{1+r_{01}{r}_{12}\exp \left(-\mathrm{i}\alpha \right)}}.$

Using Equation (25), the reflection coefficient for the air/thin film/substrate structure (Figure 3) can be obtained. Based on Equation (25), additional layers can be incorporated. When an additional layer is introduced (e.g., the buffer layer in an air/thin film/buffer layer/substrate structure), the reflection coefficient $r_{0123}$ is derived. To determine $r_{0123}$, Equation (25) is first applied to compute $r_{123}$:

(26)$r_{123}={\displaystyle \frac{r_{12}+r_{23}\exp \left(-i\alpha \right)}{1+{ r}_{12}{r}_{23}\exp \left(-i\alpha \right)}}.$

By substituting $r_{123}$, obtained from Equation (26), into the $r_{12}$ term in Equation (25), the expression for $r_{0123}$ is obtained as follows:

(27)$r_{0123}={\displaystyle \frac{r_{01}+r_{123}\exp \left(-i\alpha \right)}{1+{ r}_{01}{r}_{123}\exp \left(-i\alpha \right)}}.$

This recursive approach, based on Equation (25), allows us to determine the reflection coefficient of a multilayer structure with an arbitrary number of layers. Using Equations (10) and (27), we can obtain the reflectance of multilayered structures for a specific single wavelength. To simulate the color, we need to obtain the reflectance spectrum in the visible range, which thus requires us to use different refractive indices and extinction coefficients depending on the wavelength (i.e., dispersion) for the visible range.

2.1.3. Color Simulation

To convert a reflectance spectrum to a color (i.e., sRGB value), the spectrum is first mapped to the CIE 1931 XYZ color space using formulas based on color matching functions [26]:

(28)$N=\int \stackrel{-}{y}\left(\lambda \right)I\left(\lambda \right)d\lambda$

(29)$X={\displaystyle \frac{1}{N}}\int \stackrel{-}{x}\left(\lambda \right)S\left(\lambda \right)I\left(\lambda \right)d\lambda$

(30)$Y={\displaystyle \frac{1}{N}}\int \stackrel{-}{y}\left(\lambda \right)S\left(\lambda \right)I\left(\lambda \right)d\lambda$

(31)$Z={\displaystyle \frac{1}{N}}\int \stackrel{-}{z}\left(\lambda \right)S\left(\lambda \right)I\left(\lambda \right)d\lambda$

where $\lambda$ is the wavelength and $\overline{x}$, $\overline{y}$, and $\overline{z}$ are the color-matching functions, which are numerical descriptions of the chromatic and tristimulus response of cone cells [27]. The parameters $S\left(\lambda \right)$ and $I\left(\lambda \right)$ are the spectral reflectance and standard illuminant, respectively. In this study, we employed a CIE 1931 2° standard observer and D65 white reference as the color matching functions and the reference illuminant, respectively [28,29]. The CIE 1931 2° standard observer quantifies human color perception. Due to the distribution of cones in the human eye, the tristimulus values of colors vary depending on the observer’s field of view. The standard observer model was defined based on stimuli subtending a 2° visual angle, corresponding to the foveal region of the retina. The D65 white reference is the standard daylight illuminant, positioned using CIE values. The values X, Y, and Z can be converted to either CIE xy values or sRGB values. For CIE xy values,

(32)$x={\displaystyle \frac{X}{X+Y+Z}}$

(33)$y={\displaystyle \frac{Y}{X+Y+Z}}$

where x and y are the coordinates of the CIE 1931 color space chromaticity diagram, which is widely used to specify colors in practice. CIE XYZ coordinates are linearly transformed to sRGB coordinates using the following conversion matrix [30]:

(34)$\left[\begin{array}{c}r\\ g\\ b\end{array}\right]=\left[\begin{array}{ccc}3.24045& -1.53714& -0.49853\\ -0.96927& 1.87601& 0.04156\\ 0.05564& -0.20403& 1.05723\end{array}\right]\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]$

where r, g, and b are linear RGB. The sRGB color space defines its values in a nonlinear form that corresponds to human visual perception. However, the conversion from XYZ to RGB derived above produces linear RGB values. To obtain sRGB values, a gamma-compressed transformation is used [30]:

(35)$V=\left\{\begin{array}{cc}12.92v& v \le 0.0031308\\ 1.055v^{1/2.4}-0.055& \mathrm{otherwise}\end{array}\right.$

where $v\in \left\{r, g, b\right\}$ and $V\in \left\{R, G, \mathit{B}\right\}$. Since the values obtained using Equation (35) are in the nominal range (0–1), we multiply each component by 255 to represent conventional sRGB values in the range of 0–255.

2.2. Python Coding

Based on the aforementioned theoretical background, we developed a Python 3.9.13 code to simulate the colors of multilayered thin films. We imported libraries—NumPy, pandas, SciPy, and matplotlib—for numerical computations, data manipulation, interpolation, and visualization, respectively. A wavelength array in the visible range with a fixed step size was generated using the arange function in the NumPy library (line 2 in Listing 1). To use the dispersion relation of each layer, we needed to import the refractive indices (n) and extinction coefficients (k) of the wavelengths in the array (lines 4–10 in Listing 1). Using this code, we aimed to simulate the color of SnO₂ layers on a Si substrate, which has a native oxide layer; thus we used the n and k values of SnO₂, SiO₂, and Si from the literature [31,32,33,34]. We imported a standard illuminant (D65) and applied it to the obtained reflectance spectra to simulate the spectra from the sample, and the color-matching functions (CMFs) were imported and applied to extract RGB values from the simulated spectra (lines 35–42 in Listing 1). Here, missing n, k, and CMFs values, which are required but absent in the literature, can be estimated by interpolation using the interp1d function in the SciPy library. Since both the extinction coefficients and dispersion of air in the visible range are negligible, we set the refractive index of air as the constant, 1.0003, and defined the imaginary unit as j (lines 12–15 in Listing 1).

Listing 1. Python code for data preprocessing.

The reflection coefficients of the p- and s-polarized waves and the reflectance are calculated in lines 1–10 and lines 11–45 of Listing 2, respectively. Here, N_j and N_k indicate the complex refractive indices of media j and k, respectively. ang_j and ang_k are the angles of incidence and transmission, respectively. In line 11 of Listing 2, wl, N_0, N_1, N_2, and N_3 are the wavelength and the complex refractive indices of the atmosphere (i.e., air), top, middle, and bottom layers (here the bottom layer corresponds to the substrate), respectively. d_1, and d_2 are the thicknesses of the top and middle layers, respectively. The mathematical expression corresponding to each line of code is given above it using the equation numbers described in Section 2.1 (Theoretical Background).

Listing 2. Python code used to calculate the reflection coefficients and the reflectance.

By modifying the values of d_1, d_2, and the angle in line 11 of Listing 2, we can obtain reflectance spectra for various thicknesses and incident ray angles. To convert the reflectance spectra into RGB values, we apply Listing 3; in line 1, R_n is the calculated reflectance spectrum, while D65, CMFs_X, CMFs_Y, and CMFs_Z are the standard illuminant and the color matching functions for X, Y, and Z, respectively, which were imported and defined earlier (in line 38–42 of Listing 1).

Listing 3. Python code used to convert the reflectance spectrum to RGB values.

3. Results and Discussion

To validate the Python code used for simulating the color of multilayered structures, we simulated the colors of three different oxide thin films (SiO₂, SnO₂, and ZnO) deposited on a Si substrate as a function of the films’ thickness and compared them with the literature and our own experiment. We first checked the color of the SiO₂ thin film on a Si substrate using a function of the film thickness. It can be clearly seen that the obtained color chart is very consistent with the reported simulation result (Figure 4a) [35]. For the comparative analysis with our experiment, we fabricated a SnO₂ film with a thickness gradient on a Si wafer using a sputtering method and measured its thicknesses by performing line profile measurements across a masked area using noncontact (tapping mode) atomic force microscopy (AFM; Icon-PT-PLUS, Bruker, Germany). For each point, three independent line profile measurements were taken to obtain average thickness and estimate the standard deviation. The step size of the line profile measurements was 1.5 μm, and the standard deviation of the thickness values was generally around 0.3 nm. The base pressure was ≈3 × 10⁻⁸ Torr, and the working pressure of 10 mTorr was adjusted using Ar (99.9999%). The films were grown at room temperature. We assumed the presence of a 2 nm thick native oxide layer (SiO₂) between the film and substrate [36,37]. The two white dashed lines in Figure 4b indicate the experimentally measured thickness range (i.e., 76–385 nm). Although the thickness gradient of the experimental samples is nonlinear with respect to the position of each measured point, the colors and the thicknesses measured at the nine points are in good agreement with the theoretically predicted colors and thicknesses.

Figure 5a shows the experimental colors of a ZnO thin film with a thickness gradient on a Si substrate, taken at different incident angles (5, 30, and 45°). Similarly to the SnO₂ thin film, we fabricated the ZnO film using the sputtering method and measured its thicknesses using AFM. As the incident angle increases, the overall color gradient shifts toward the thicker region, and as the thickness increases, the gradient becomes more pronounced. For instance, the color of the film changes from blue to magenta at 240 nm and from cyan to pink at 410 nm as the angle varies. To visualize the change in the color gradient, we simulated the colors of the ZnO thin film on a Si substrate, with a 2 nm thick SiO₂ layer in between the two, as a function of the thickness and angle (Figure 5b). In this simulation, we used the n and k values of ZnO from the literature [38]. The simulation results obtained using the Python code describe the experimental color variations as depending on both thicknesses and incident angles. The three white boxes in Figure 5b represent the color ranges expected for the given thicknesses and angle variations in the experimental sample. The theoretical results clearly demonstrate the color transition seen due to the change in angle: from cyan to pink in the thicker regions and with an expansion of the yellow region in the thinner regions.

So far, we have obtained CIE XYZ values from simulated reflectance spectra and have converted them into sRGB for more intuitive color representations. Additionally, using Equations (32) and (33), we can transform CIE XYZ into CIE xy, allowing us to map the colors onto the CIE 1931 color space. This two-dimensional color space facilitates the identifying of color differences and supports a device- and observer-independent color representations suitable for various display and printing applications. Figure 6 shows the simulated colors of the SnO₂ and ZnO thin films represented in both the sRGB and CIE xy color spaces. In the sRGB representation, the fact that the color variations are a function of thickness is easily observable. In contrast, the CIE xy representation provides insight into the color gamut achievable from a given material. For example, although the colors of SnO₂ and ZnO appear similar in the RGB color space, the CIE color space reveals that SnO₂ can be used to produce a slightly more vivid blue.

The results above demonstrate the ability of the Python code to create accurate reflectance spectra and corresponding color simulations. This code could be developed further to estimate the thickness of a film based on known n and k values, and vice versa. Moreover, recent advances in materials research have increasingly made use of AI-driven approaches, such as optics-based inverse design [39,40,41,42,43,44,45] and active learning [46]. In this context, this code, implemented in Python, which has been primarily used in AI research, serves as an accessible and easily modifiable, yet powerful, tool.

4. Conclusions

We demonstrated that our Python script can simulate the color of multilayered structures, with a detailed theoretical background and annotations provided for readers as well as users. The developed code was validated by applying it to well-known oxide materials (SnO₂ and ZnO). We believe that anyone can adjust its structural and material parameters for their own system and research by following the instructions presented here. It is also worthwhile noting that the code generates a substantial quantity of data very rapidly—for example, simulating the 54,000 colors displayed in Figure 5b took merely 52 s (using CPU: Intel i7-1255U; GPU: Intel Iris Xe Graphics). This code can be further used and developed for various purposes—especially for machine learning, which requires a large amount of spectral and color data for model training.

Author Contributions

Conceptualization, S.L.; methodology, D.L. and S.L.; validation, D.L.; formal analysis, D.L.; investigation, D.L.; writing—original draft preparation, D.L. and S.L.; writing—review and editing, D.L. and S.L.; visualization, D.L.; supervision, S.L.; project administration, S.L.; funding acquisition, S.L. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data and their implementation details are available on GitHub (https://github.com/CQMLab/Color-simulation-python-code-and-examples, accessed on 24 April 2025).

Conflicts of Interest

The authors declare no conflicts of interest.

Abbreviations

The following abbreviations are used in this manuscript:

SnO2	Tin dioxide
ZnO	Zinc oxide
SiO2	Silicon dioxide
Si	Silicon
RGB	Red, Green, and Blue
AI	Artificial intelligence
sRGB	Standard Red, Green, and Blue
CIE	Commission International de l’Eclairage
CMFs	Color-matching functions

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1 Incident wave at angle θ_i hits the interface between two media with refractive indices n₁ and n₂: (a) the electric field (E) parallel to the plane of incidence (p-polarization) and (b) E perpendicular to the plane of incidence (s-polarization). E and the magnetic field (B) are mutually perpendicular (i.e., the propagation vector k = E × B). The subscripts i, r, and t represent the incident, reflected, and transmitted waves, respectively.

Figure 2 Wave propagation in a medium with (a) zero and (b) a nonzero extinction coefficient (k). When the extinction coefficient is zero (a), the amplitude (E_t) remains constant. In contrast, with a nonzero extinction coefficient (b), the amplitude exhibits exponential decay.

Figure 3 Optical interference in a thin film. The interference between the reflected light from both the upper and lower boundaries needs to be taken into account. θ_i and θ_t are the angles of incidence and refraction, respectively. θ_i and θ_t are equal to the ∠ACD and half of the ∠ABC, respectively, according to the reflection law. d is the thickness of the film. Each arrow represents the path of a reflected or transmitted wave, and the adjacent coefficients $r_{mn}$ ($t_{mn}$) indicate the amplitude reflection (transmission) coefficient for the light incident from medium N_m to medium N_n. The topmost coefficients along each ray represent the product of the corresponding amplitude coefficient and the phase difference accumulated along that path.

Figure 4 Comparison of simulated colors with literature and experiment: (a) SiO₂ thin film on a Si substrate: simulated colors as a function of the SiO₂ film thickness from literature (upper) [35] and by Python code in this study (lower). (b) SnO₂ thin film on a Si substrate: the photograph of the SnO₂ thin film with thickness gradient (upper) and simulated color as a function of the SnO₂ film thickness by the Python code (lower). The grid in the photograph is a masked area for line profile measurements—the value below each dot is the film thickness from the line profile measurements nearby. The thickness of a native SiO₂ layer was set to 2 nm in the simulation. The two white dashed lines correspond to the identified thickness range of the experimental samples (i.e., 76–385 nm).

Figure 5 Comparison of experimental colors and simulated colors of the ZnO thin film. (a) Experimental colors of the ZnO thin film on a Si substrate: photographs of the ZnO thin film with a thickness gradient taken from various angles (5, 30, and 45°). The dots marked on the Si substrate are approximately 1 cm apart and the values below each dot correspond to the thicknesses obtained from line profile measurements nearby. (b) Simulated colors of the ZnO thin film on a Si substrate, using the Python code, as a function of the thickness and incident angle. The thickness of the native SiO₂ layer was set to 2 nm in the simulation. The three white dashed boxes correspond to the color ranges expected from the experimental samples (thickness range: 200–410 nm; incident angle: 5, 30, and 45°).

Figure 6 Simulated color gamuts of SnO₂ and ZnO thin films on a Si substrate, presented as sRGB (upper, colors as a function of thickness) and CIE 1931 xy chromaticity diagrams (lower): (a) SnO₂ and (b) ZnO. The thickness of the films was simulated in the range of 0 to 2000 nm, with an incident angle of zero. The points in the CIE xy diagram correspond to the obtained sRGB values demonstrated in the color bar.