Semi-analytical technique for the design of disordered coatings with tailored optical properties

Disordered coatings, which consist of dielectric/metal nanoparticles dispersed randomly in a matrix, find their application in several fields such as solar thermal absorber coatings [1], solar reflecting coatings [2], color paints [3], translucent paints [4], tissue optics [5], daytime passive radiative cooling coatings [6, 7, 8], and many more. The main advantages that such disordered media offer that make them an attractive proposition for use in these applications are their cost-effective means of fabrication, and tunability of the desired optical properties of the coating - since the spectral position of Mie (plasmon) resonance of the embedded dielectric (metal) particles strongly depend on the size of the particles. The main challenging task in the design of such disordered media is the modelling of its optical properties. Techniques based on homogenization of the composite structures that predict effective permittivity and permeability of the disordered media, such as Maxwell-Garnett theory [9] and Bruggeman’s model [10], are valid only when the particle sizes are much smaller than the incident wavelength [11]. Doyle et al. [12] showed that the use of Mie coefficients in this effective medium theory provides good accuracy to the calculation of effective optical properties of metal spheres suspended in a polymer. However, the theory predicts absorption for a nonabsorbing particle [11] and thus cannot be used to predict the effective refractive index of disordered media for solar reflecting paint/coatings where non absorbing particles are utilized. In literature, other analytical techniques developed for this objective include those which consider diffusion of photons [13], and those which solve the radiative transfer equation under $N$-flux ($2\leq N\leq 4$) approximations [14, 15]. Of these methods Kubelka Munk (KM) theory [16] (for which $N=2$) and the four-flux (FF) method [17] are commonly used. In addition, simulation techniques such as the Monte Carlo method [18], and exact electromagnetic solvers are employed to model the optical/radiative properties of disordered coatings. However, these simulation techniques do not present a clear picture linking the microscopic properties of the particles, such as the scattering and absorption coefficients, to the macroscopic optical properties of the coating. Moreover, exact electromagnetic solvers which solve Maxwell’s equations numerically to obtain radiative properties of the coating, put a premium on computational resources and the time for design when the thickness of random media is in the order of tens/hundreds of microns - as is currently being deployed in these applications. Particularly when several parameters of the configuration are in play, such as that encountered in disordered media, analytical techniques such as KM and FF theories provide important means to arrive at an optimum combination of the parameters with minimal computational resources while also explicitly linking the properties of the micro constituents to the observed optical properties of the coating.

KM and FF theories have so far in literature been used in applications where the spectrum of interest has been limited to either the visible spectrum or IR separately. For example, KM theory has been used extensively in paints and coatings [1, 3], paper industry [19], tissue optics [5] among others. Similarly, the FF method has been used extensively by researchers to model, predict and optimize the optical properties of light scattering coatings [2, 20, 21, 22]. However, the applicability of these theories over a broad spectrum covering both visible and IR spectrum simultaneously has not been a subject of attention. This becomes important when designing coatings for applications such as day-time passive radiative cooling and solar thermal absorber plate coatings where tailored optical properties over a broad spectrum covering both the solar spectrum as well as far-infrared are crucial. For example, coatings for day-time passive radiative cooling [23] require high reflectivity in solar spectrum (0.3-3 $\mu$m wavelength range) and high emissivity in the infra-red spectrum (5-15 $\mu$m wavelength range). It is not obvious that the analytical techniques retain their accuracy over such a broad spectrum since with increasing wavelength there is a possibility that the nature of scattering transitions from the independent scattering regime (where scattering cross-section of particles can be superposed) to dependent scattering regime (where near-field effects, and interference between far-field scattered radiation become important). Previously reported relation [24] demarcating the two regimes has been obtained from experimental observations carried out in the visible spectrum only. Thus there is a need to explore the applicability of these analytical techniques over a broad spectrum in greater depth. In regimes where the predictions from these analytical techniques are not satisfactory, other possible methods of design which combine the accuracy of exact electromagnetic solvers with the minimal computational requirements of analytical methods are expected to be of pressing need to researchers interested in designing such coatings.

With this in mind, the manuscript has been arranged as follows. In Section 2 we have compared the reflectivity and emissivity predictions of KM and FF techniques with results from exact numerical solvers for different degrees of absorption in the particles (imaginary index of particle = 0.0, $10^{-2}$, and $10^{-1}$) and in the matrix (imaginary index of matrix = 0.0, $10^{-4}$, and $10^{-2}$) for different thickness of the coating (10 $\mu$m and 50 $\mu$m) in the wavelength range $0.3-15~{}\mu m$. We show that these techniques are accurate over the entire spectrum when particles are in the limit of independent scattering and under low absorption conditions but fail when volume fraction of particles is high such that interaction among particles is no longer negligible or when absorption in the matrix/particles is high. For such conditions where analytical techniques fail to predict the optical properties accurately we propose an alternative technique which combines the use of exact numerical solver and KM theory which we show can predict optical properties with accuracy as well as with minimal computational requirements. This ‘semi-analytical’ technique is detailed in Section 3. In the end, as an example to showcase the applicability of this semi-analytical technique, we predict the properties of a disordered coating suitable for the application of passive radiative cooling and compare these with experimental measurements previously reported in literature.

We start with the expressions for reflectivity and transmissivity of the coating as obtained from KM and FF theories which we use in the work to analyze the optical properties of the disordered coating. Detailed derivations of these expressions can be found in several references [25, 17, 26]. The optical coating considered in this work is a plane-parallel slab of particulate composite on a substrate as shown in Fig. 1. The composite is considered to be of finite thickness and infinite extension in the lateral direction. The randomly distributed spherical particles embedded within the host medium (also called the matrix) act as inhomogeneities to the propagating EM wave, thereby causing its scattering (and absorption, in case the particle is lossy). The objective is to predict the optical properties of this coating including the total reflectance, transmittance and absorption.

The expression for the reflectivity ($R_{\text{KM}}$) and transmissivity ($T_{\text{KM}}$) from KM theory is given by [25, 3]:

$$R_{\text{KM}}=\frac{(1-\gamma)(1+\gamma)(\exp(AL)-\exp(-AL))}{(1+\gamma)^{2}\exp(AL)-(1-\gamma)^{2}\exp(-AL)}$$

(1)

$$T_{\text{KM}}=\frac{4\gamma}{(1+\gamma)^{2}\exp(AL)-(1-\gamma)^{2}\exp(-AL)}$$

(2)

where, $L$ is the thickness of the layer, the coefficients $\gamma$ and $A$ are given by [3, 27]:

$$\gamma=\sqrt{k/(k+s^{\prime})};\,\,A=\sqrt{2k(2k+s^{\prime})}$$

(3)

with

$$s^{\prime}=\frac{3s(1-g)-k}{4};\,\,$$

(4)

and the factors $s$ and $k$ got using Mie theory [1]

$$s=\frac{3fQ_{\text{sca}}}{4r};\,\,\,\,\,\,k=\frac{3fQ_{\text{abs}}}{4r},$$

(5)

where $f$ is the volume fraction, $r$ is the radius of the sphere, $Q_{\text{sca}}$ $(Q_{\text{abs}})$ is the Mie scattering (absorption) efficiency of a single particle embedded in a host medium of index $n_{h}$, and $g$ is the asymmetry parameter. Expressions for $Q_{\text{sca}}$, $Q_{\text{abs}}$, and $g$ in terms of standard Mie coefficients can be found in Ref. [28]. It should be pointed out that the relations between the coating properties $\gamma$ and $A$, and the particle properties $s$ and $k$ given in Eq. 3 are not unique - several other relations [29, 30, 31, 32, 4, 5, 33] have been proposed over the years. The expressions in Eq. 3 and Eq. 4, taken from Ref. [3, 27], is representative and have been chosen for demonstrative reasons. As we will see in Sec. 3 the semi-analytical method being proposed in this work does not depend on such relations and hence do not affect the central results of this work.

In the limit of low absorption, $kL\rightarrow 0$, the reflectivity in Eq. (1) can be shown to reduce to [25]:

$$R_{\text{KM}}=\dfrac{s^{\prime}L}{s^{\prime}L+1}.$$

(6)

It must be noted that Eq. (1) and Eq. (2) do not take into account surface reflection of incident radiation at interface (1). Modified reflectance $R_{0}$ and transmittance $T_{0}$ which take into account surface reflection correction are calculated using [34]:

$$R_{\text{0}}=R_{\text{c}}+\frac{(1-R_{\text{c}})(1-R_{\text{i}})R_{\text{KM}}}{1-R_{\text{i}}R_{\text{KM}}};\,\,\,\,\,T_{\text{0}}=\frac{(1-R_{\text{c}})T_{\text{KM}}}{1-R_{\text{i}}R_{\text{KM}}}$$

(7)

where, $R_{\text{c}}$ is the specular reflectance of incident light got from Fresnel reflection which for normal incidence from a medium of index $n_{\text{surr}}$ reads:

$$R_{\text{c}}=\left(\frac{n-1}{n+1}\right)^{2}$$

(8)

with $n=n_{\text{h}}/n_{\text{surr}}$ and $R_{\text{i}}$ is the diffuse reflectance of internal radiation at interface (1), marked in Fig. 1, which is calculated using:

$$R_{i}=2\int_{0}^{\pi/2}\rho(\theta)\,\text{sin}\theta\,\text{cos}\theta\,\text{d}\theta.$$

(9)

where, from Fresnel’s coefficients:

$$\rho(\theta)=\frac{1}{2}\left[\left(\dfrac{\sqrt{n^{2}-\text{sin}\theta}-\text{cos}\theta}{\sqrt{n^{2}-\text{sin}\theta}+\text{cos}\theta}\right)^{2}+\left(\dfrac{n^{2}\text{cos}\theta-\sqrt{n^{2}-\text{sin}\theta}}{n^{2}\text{cos}\theta+\sqrt{n^{2}-\text{sin}\theta}}\right)^{2}\right].$$

(10)

The expression for $R_{\text{i}}$ from Eq. 9 can be used even the limit of low diffuse scattering since the contribution from the product $R_{i}R_{\text{KM}}$ will be negligible in this regime.

Many configurations developed for radiative cooling application [7, 35, 36] and solar absorber plates [37, 38] involve use of a substrate. In the presence of a substrate, the net reflectance and transmittance from Eq. 7 will have to be further modified as [39]:

$$R=R_{\text{0}}+\frac{T_{\text{0}}^{2}R_{\text{g}}}{1-R_{\text{0}}R_{\text{g}}};\,\,\,T=\frac{(1-R_{\text{g}})T_{\text{0}}}{1-R_{\text{0}}R_{\text{g}}}$$

(11)

Here $R_{\text{g}}$ is the diffuse reflectance at interface (2) obtained from Eq. 9 with $n=n_{\text{h}}/n_{\text{g}}$. The substrate index $n_{g}$ is taken to be 1.5 in this work.

The derivation of reflection and transmission coefficients from KM theory assumes that the incident light is diffuse. When the incident radiation is collimated, alternate methods such as the four-flux theory, which take into account the propagation of both collimated and diffuse radiation across the interfaces in two directions, are expected to be more accurate. This careful consideration of both collimated and diffuse components leads to expressions for the optical properties being far more complicated than in KM theory. The net reflection and transmission coefficients when incident radiation is fully collimated can be expressed in terms of a summation over collimated-collimated reflectivity $(R_{\text{cc}})$, collimated-diffuse reflectivity $(R_{\text{cd}})$, collimated-collimated transmissivity $(T_{\text{cc}})$, and collimated-diffuse transmissivity $(T_{\text{cd}})$ as:

$$R=R_{\text{cc}}+R_{\text{cd}}+R_{\text{dd}};\,\,\,\,T=T_{\text{cc}}+T_{\text{cd}}+T_{\text{dd}}$$

(12)

Expressions for $R_{\text{cc}}$, $R_{\text{cd}}$, $T_{\text{cc}}$ and $T_{\text{cd}}$ are quite elaborate and have been included in the supplementary document (Section S1) for reference.

In Sections 2.1, 2.2, and 2.3, we use the expressions for $R$ and $T$ given in Eqs. 11 for KM theory and Eqs. 12 for FF to predict the optical properties of disordered coatings and compare these with the results obtained from Lumerical FDTD solver [40]. We analyze situations where both the particles and host medium are absorbing as well as non-absorbing, and also consider the effect of different thickness of the coating. The degree of absorption in particles considered in this work are relevant for dielectric inclusions typically included in coatings for use in radiative cooling and solar thermal applications. In addition, to facilitate the parametric study we assume non-dispersive form of refractive index for both the particles as well as the host matrix. We first confine our analysis to the independent scattering regime in Sections 2.1 and 2.2, and extend the analysis to dependent scattering regime later in Sec. 2.3. The FDTD simulations were set up in ANSYS Lumerical. Periodic boundary conditions were applied in the lateral $x$ and $y$ directions, and coating is illuminated with a plane wave source from $z$ direction. A mesh size of 30 nm was used which we find is sufficient for convergence (mesh convergence study is shown in supplementary Fig. S3).

The comparison with FDTD simulations shown in Section 2 demonstrate the failure of KM and FF analytical methods in configurations where dependent scattering is not negligible and when matrix/particles are absorbing. This failure can be attributed to the actual scattering and absorption coefficients of these coatings diverging from the values calculated using Mie scattering coefficients of the individual particles. At present no single analytical technique exists that can correctly predict the optical properties of particulate media in the presence of dependent scattering effects as well as correctly account for the absorption in matrix/particles. One can then resort to using exact numerical solvers to accurately estimate the optical properties of the coating in such cases. However, as Fig. 6 shows, the computational time required to simulate such structures increases exponentially with thickness of the coating. For coatings of thickness in the range 100-500 microns which are currently being adopted in literature for the radiative cooling application [6, 49, 46, 36, 8] the design time is clearly prohibitive.

In such cases it becomes imperative to develop alternate techniques which can combine the accuracy power of exact FDTD solvers with the simplicity and minimal computational requirements of the analytical techniques. Particularly when multiple parameters are involved in design - such as that observed for disordered media - such a method will prove to be useful in reducing the design time to find the optimum combination of parameters necessary to obtain the required optical properties of the coating. In order to obtain a better estimate for the absorption and scattering coefficients of such media where dependent scattering effects are non-negligible, researchers have previously [27, 51, 52, 53] relied on experimental measurements of the optical properties of a fabricated coating and then using the KM theory results from Section 2 to extract the required coefficients. Instead of relying on experimental measurements which is not always feasible especially at the initial state of design, we modify this technique and instead propose the following two-step semi-analytical method to estimate the optical properties of random media of thickness $L$ when usage of exact numerical solvers to simulate the properties of such a thick coating is prohibitive.

• •

  Step 1: Use a numerical solver to obtain the optical properties $R$ and $T$ of a similar coating but with much smaller thickness $t_{s}\ll L$ and extract the $\gamma$ and $A$ parameters by inverting Eq. 1 and 2. Care must be taken at this step to ensure that the configuration set up in the solver considers incident light to be in the same medium as the index of matrix i.e., $n_{\text{surr}}=n_{\text{h}}$ in order to ensure that reflection from surfaces and substrates are not included in this step. In case the host matrix is absorbing then only the real part is considered i.e., $n_{\text{surr}}=\text{Re}(n_{\text{h}})$. Care must also be taken to ensure that when scattering efficiency of the particles is high, the value of $t_{s}$ should be chosen such that $t_{s}\gg l_{s}$ where $l_{s}\approx 1/(N\sigma_{s})$ is the scattering mean-free path with $N$ being the particle number density and $\sigma_{s}$ the scattering cross section. At the other limit when scattering efficiency is low the optical properties of the coating are primarily determined from surface reflection and transmission which are accounted for in step 2. Thus the choice of $t_{s}$ is determined from the scattering mean-free path calculated in the high-scattering regime.

• •

  Step 2: From the $\gamma$ and $A$ parameters extracted from step 1 use the analytical expressions from KM theory i.e. Eqs. 1, 2, 7 and 11 to predict the optical properties of the coating of the required thickness $L\gg t_{s}$. Specular reflection at the surfaces as well as at the substrate are accounted for here.

A more elaborate procedure, along with details of a supporting convergence test which may need to be incorporated in some cases to arrive at the value of thickness $t_{s}$ is included in Section S4 of supplementary.

We now apply this technique for the cases considered in Section 2 where the predictions from analytical methods deviated significantly from those of FDTD solver, such as for the dependent scattering regime, as well as when the absorption in the particles/host matrix is significant. Fig. 7 (Fig. 8) shows the comparison between the predictions from the semi-analytical technique and from FDTD simulations when absorption in particles (host matrix) is varied. In both these cases the semi-analytical technique uses the results of exact FDTD simulations of a $10~{}\mu$m thick coating to predict the optical properties of a larger $50~{}\mu$m thick coating. A volume fill fraction $f=0.3$ is maintained in both these cases where dependent scattering effects are known to be dominant, while keeping other parameter values same as that analysed for the monodisperse case of Sec. 2.1.

For these cases we observe a close match in the predictions of the semi-analytical method with the FDTD results over the entire spectrum, with only a slight deviation observed for the higher wavelengths when absorption is high. The weighted-average reflectivity of the coating for solar spectrum, and emissivity over the infra-red spectrum for the cases considered in Figs. 7 and 8 are listed in Table 2 along with the deviation from FDTD simulations expressed in % error in brackets. Particularly illustrative of the effectiveness of the semi-analytical technique is the reduction in error (1.03 % in Sr. no. 1) as compared to those obtained from analytical techniques and reported in Table 1 ( 8.13 % using KM theory and 6.98 % using FF theory in Sr. no. 13) for the configuration: $n_{\text{p}}=2.5$, $n_{\text{h}}=1.5,$ $r=0.25~{}\mu\text{m},$ $f=0.3,$ and $L=50~{}\mu\text{m}$ where dependent scattering is expected to be dominant.

We now apply the semi-analytical technique described in Section 3 to predict the optical properties of fabricated coatings reported in literature which have been designed for radiative cooling application. We choose two such disordered coatings where dependent scattering is expected to be dominant so that analytical techniques are not applicable, and the thickness of the coating prohibits the use of exact electromagnetic solvers to predict the optical properties to good accuracy.

In Ref. [48], a hierarchically porous polymer (P(VdF-HFP)) coating of thickness 300 $\mu$m containing air voids with sizes ranging from 0.05-5 $\mu$m in P(VdF-HFP) matrix has been fabricated, and experimentally characterized to have solar reflectivity value of 0.96 and emissivity in the 8-13 $\mu$m wavelength range to be 0.97. In order to apply semi-analytical technique to predict the properties of this coating, we set up a simulation in FDTD solver with a smaller coating thickness $t_{s}=50~{}\mu$m (determined using the convergence test explained in Section S4 of supplementary). This thickness is chosen to ensure sufficient number of larger sized air voids ($r\approx 2.5~{}\mu$m) in this P(VdF-HFP) matrix. The size distribution of nano-micro air voids used in the simulation is given in supplementary (Fig. S4). Refractive index data of P(VDF-HFP) is extracted from Ref. [48]. The reflectivity data in the wavelength range $0.3-16~{}\mu$ m, predicted using the semi-analytical method for $L=300~{}\mu$m thickness, is compared with that reported in Ref. [48] in Fig. 9(a). While an appreciable match is noticed in the predicted values across the spectrum, small deviation observed in the reflectivity values can be attributed to our inability to incorporate exact size distribution of both micro and nano voids as present in the fabricated structure, in ANSYS Lumerical.

In Ref. [46] an ultrawhite BaSO₄ film of thickness 400 $\mu$m has been developed with 60 % volume fraction of BaSO₄ nanoparticles, and has been characterized to have reflectivity of 0.976 in the solar spectrum and emissivity of 0.96 in 8-13 $\mu$m wavelength range. In order to apply the semi-analytical technique to predict the properties of this coating, we set up a simulation in FDTD solver with structure thickness $t_{s}=15~{}\mu$m and BaSO₄ spherical particles randomly distributed with volume fraction 60 %. The particles are taken to be of uniform size distribution with diameters spread over the range $398\pm 130$ nm to match that reported in Ref. [46]. Matrix is considered to be air for BaSO₄ film. Refractive index data of BaSO₄ is extracted from Ref. [54]. The emissivity data in the wavelength range $0.3-16~{}\mu$m, predicted using the semi-analytical method for $L=400~{}\mu$m thickness, is compared with that reported in Ref. [46] in Fig. 9(b). While we again observe an appreciable match across the spectrum, some deviation observed particularly around wavelength of $2~{}\mu$m is suspected to be due to difference in the refractive index of the fabricated film and that calculated from first-principles in Ref. [54].

In this study we have analyzed the applicability of well-known analytical techniques of KM and FF theories to predict optical properties of a disordered metamaterial coating over a broad spectrum ranging from 300 nm to $15~{}\mu\text{m}$ wavelength. Recent advancements in the use of disordered coatings in applications such as radiative cooling and solar thermal absorber plates which require tailored optical properties over this wavelength range necessitates such a study. Based on deviations observed between the predictions of these analytical techniques and exact FDTD solver in the dependent scattering regime, a two-step semi-analytical technique has been proposed which can be used to predict optical properties of such coatings with good accuracy and minimal computational resources. Such a method is expected to be resourceful for designing coatings with specific optical properties where several parameter combinations need to be investigated to arrive at an optimal combination. Small deviations observed when absorption in host matrix is high warrants further research to improve this technique.

B.R.M. acknowledges support from Prime Minister’s Research Fellowship (PMRF). K.S. acknowledges support from La Fondation Dassault Systèmes and SERB Grant No. SRG/2020/001 511.

The authors declare no conflicts of interest.

See Supplement 1 for supporting content.