# Dry powders reflectance model based on enhanced backscattering: case of hematite α-Fe2O3

Morgane Gerardin, Pauline Martinetto, Nicolas Holzschuch

► To cite this version:

Morgane Gerardin, Pauline Martinetto, Nicolas Holzschuch. Dry powders reflectance model based on enhanced backscattering: case of hematite α-Fe2O3. Journal of the Optical Society of America. A Optics, Image Science, and Vision, 2023, 40 (9), pp.1817-1830. 10.1364/JOSAA.487498. hal-04316639v2

HAL Id: hal-04316639

https://inria.hal.science/hal-04316639v2

Submitted on 1 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Dry powders reflectance model based on enhanced back-scattering: case of hematite

α – Fe₂O₃

MORGANE GERARDIN,¹,²,* PAULINE MARTINETTO,² AND NICOLAS HOLZSCHUCH¹

¹Univ. Grenoble Alpes, CNRS, Inria, Grenoble INP, LJK, 38000 Grenoble, France

²Univ. Grenoble Alpes, CNRS, Institut Néel, 38000 Grenoble, France

*gerardin.morgane@gmail.com

**Abstract:** By performing Bidirectional Reflectance Distribution Function (BRDF) measurements, we have identified back-scattering as the main phenomena involved in the appearance of dry nano-crystallized powders. We introduce an analytical and physically-based BRDF model that relies on the enhanced back-scattering theory to accurately reproduce BRDF measurements. These experimental data were performed on optically thick layers of dry powders with various grains' morphologies. Our results are significantly better than the ones obtained with previous models. Our model has been validated against the BRDF measurements of multiple synthesized nano-crystallized and mono-disperse α – Fe₂O₃ hematite powders. Finally, we discuss the ability of our model to be extended to other materials or more complex powder morphologies.

© 2023 Optica Publishing Group

## 1. Introduction

Accurate reproduction of powdered material appearance is an important issue in many domains, such as environmental sciences, teledetection, digital prototyping, entertainment and art preservation. Powder materials are ubiquitous in real life. They can be used with a binder like pigments in most painting techniques and in some cosmetics e.g. lipstick. But they are also used without any binder. Dry powders are found as spices, dry soils, dust, and dry pigments for example.

The latter are particularly important as they are involved in parietal paintings from paleolithic to neolithic era. In some cases, the pigments were used pure and dry on the wall. Other application techniques required the pigments to be dispersed in an organic binder [1,2]. However, only few traces of these binders are left nowadays, and they are not the main component responsible for the painting color. The parietal paintings that are visible nowadays can therefore be described as a dry pigment layer on top of a rock substrate. These cultural heritage artefacts have been degraded over the years by both environmental and anthropogenic causes, so that some of them are nowadays almost illegible [1]. Images of such heavily deteriorated artefacts are recorded using image enhancement techniques [3,4]. Time, weathering and human activity are responsible for the peeling off, the erasure and possible color changing of the original paintings.

Computer Graphics is a real asset for the preservation of parietal paintings. Photo-realistic image renderings using accurate appearance models would allow the reproduction of such artefacts and the reconstruction of the most erased ones. It will also promote their broadcasting to a large audience.

The colour of a powder material depends on both its chemical composition and on its morphometry, which refers to the size and the shape of its grains. This hypothesis is well known and particularly well justified in the case of pure α – Fe₂O₃ hematite powders that can show various shades from light red to purple [5,6]. The chemical composition of the material is involved in the absorption phenomena, while the morphology of the grains influences the way light is scattered by the material. Such light/matter interactions at the nanometric scale result in
the macroscopic material appearance an observer perceives.

At the macroscopic scale, light scattering by a material is modeled using a *Bidirectional Scattering Distribution Function* (BSDF) model. When there is no light transmission and only light scattered in reflection is at stake, the BSDF is reduced to *Bidirectional Reflectance Distribution Function* (BRDF). Most models assume that the light is mainly scattered at the surface of an object, around the specular direction ie. in the direction opposite to the incident direction. Depending on the roughness of the surface, a diffuse look is obtained. According to these surface models, the color of the rendered materials only depends on their chemical composition through the optical constants. These models are not suitable in the case of dry powders for two reasons. First, powders are a medium composed of a collection of randomly oriented scatterers, so light is scattered multiple times in the volume of the medium. To account for volume scattering, one may use the Mie scattering theory [7, 8]. This is relevant for powders diluted into a binder, and takes into account volume scattering depending on the size and shape of the grains. However, in the absence of a binder, the grains composing the powder are too close to each other, making the Mie scattering model unsuitable as shown by [6]. The second issue is that only few of surface scattering models account for scattering in the backward direction [9].

In this paper, we show that contrary to the prediction of most BRDF models, the major contribution to a dry powder BRDF is not the specular peak, but the back-scattering peak. We propose an analytical and physically-based BRDF model involving six parameters able to simulate the appearance of optically thick layers of dried pigments. We separate the model into two independent contributions: volume scattering and surface scattering. The former decomposes into single and multiple scattering that are responsible for the sharp back-scattering lobe, while the latter accounts for the specular scattering. This model has been developed and experimentally validated for seven hematite $\alpha - Fe_2O_3$ powders of different morphologies. Although our model still tends to slightly underestimate the back-scattering contribution to the BRDF, it gives a good intuition on the color variations of the studied samples. It is ready-to-use to reproduce the color of mutliple studied hematite powders through renderings.

In Section 2 we will review previous works on BRDF models suitable for the description of dry powder materials, and more particularly on models that take into account the back-scattering contribution. The theoretical background necessary for the modeling of light scattering in a particle medium is introduced in Section 3. The $\alpha - Fe_2O_3$ samples and their measurements are described in Section 4. The description of our appearance model and its optimization procedure is done on Section 5 whereas in Section 6 we present its results, showing that it fits well the experimental BRDF of the studied samples. We conclude and give some future work directions in Section 7.

## 2. Previous work

Dry powders and pigments are complex materials composed of a dense collection of particles with nanometric to micrometric grains size. In the absence of a binder, the dry grains are too close to each other to be modeled using Mie scattering theory [6]. Their appearance is the result of absorption and scattering phenomena, so its description suggests the full understanding of the interactions of light with the material, at both its surface and in its volume.

We have identified two models that could be used for dry powders' BRDF simulations: the spectralon model and the RTLSR model. Firstly, they both take into account the *back-scattering*, a phenomenon which has already been observed on dry powders such as meteorite soils [10]. Then, although these models have not been used for reproducing the appearance of dry powders or pigments, they have been introduced for describing particulate media without binder. These two reasons make them worthy to study in the case of dry powders.
## 2.1. Spectralon

An empirical model has been developed by [11] to reproduce the BRDF of the Spectralon®, a material that resembles a powder since it is made out of a collection of compressed polymer particles. This model proposes a decomposition into four independent scattering phenomena. These contributions are by order of importance: Deep scattering D, which accounts for volume scattering at small incidence and observation angles, Back-scattering B and Specular scattering S, to account for the backward and forward surface scattering, and Subsurface scattering SS, which accounts for volume scattering at more grazing incidence/observation geometries. Each of these contributions is modeled using very simple analytical expression.

$$BRDF(\theta_i, \theta_o) = D(\theta_i, \theta_o) + B(\theta_i, \theta_o) + S(\theta_i, \theta_o) + SS(\theta_i, \theta_o) \quad (1)$$

One should note that this model has been developed in the incident plane. To model the complete BRDF of a material some extension is necessary to account for the different azimuth angles. We propose such an adaptation in Supplementary work.

**Volume scattering contributions:** For small incident angles $\theta_i$, there is a higher probability for the light to penetrate deeply into the medium before it exits with a direction close to the surface normal, with a small observation angle $\theta_o$. In this case, the entry and the exit points of the light are quite close. Such a phenomenon is called Deep scattering D and is identified as the main contribution to the BRDF. Its modeling in the incident plane is inspired from an empirical expression from [12] and involves five material dependent fitting parameters $d_j$:

$$D(\theta_i, \theta_o) = d_0 + d_1(\theta_i^2 + \theta_o^2) + d_2(\theta_i^4 + \theta_o^4) + d_3\theta_i^2\theta_o^2 + d_4\theta_i^4\theta_o^4 \quad (2)$$

Deep scattering results in a large scattering lobe that describes the general shape of the measured BRDF.

For larger incident angles, light penetrates much less deeply and propagates right underneath the medium surface to an exit further away from the entry point which is modeled by the Subsurface scattering SS. Two models were proposed by [11] involving an exponential term so that the more grazing the angle, the more important this contribution is. For simplicity, we only give here the expression of the simplest, which involves two empirical parameters $ss_0$ and $ss_1$:

$$SS(\theta_i, \theta_o) = ss_0 \exp(ss_1|\theta_i + \theta_o|) \quad (3)$$

**Surface scattering contributions:** The surface scattering phenomena are responsible for the finer details of the measured BRDF. The Back scattering is modeled using a Lorentzian function with amplitude $b_0$ and width $b_1$, whereas the Specular scattering is described as a Gaussian function with width $s_1$ and whose amplitude depends on two other parameters $s_0$ and $s_2$.

As the Spectralon model involves a large number of empirical but material-dependent parameters optimized in the specific case of the Spectralon material, it can be adapted to reproduce with a high fidelity a large range of materials. Nevertheless, these parameters have no explicit relationship with the chemical composition nor the physical parameters of the considered material so they cannot be predicted, but only numerically optimized from measurement fittings.

## 2.2. RTLSR

The RTLSR (Ross-Thick Li-Sparse Reflectance) model [13, 14], is a semi-empirical model used in teledetection for simulating the BRDF of forest canopies. It decomposes the BRDF into two contributions: volume scattering contribution $R_{vol}$ deduced from an approximation of the radiative transfer theory (RTT) [15] and accounting for both single and multiple scattering between the leaves, and surface scattering $R_{geo}$ which is inferred according to geometrical
considerations [16] and accounts for shadowing and occlusions of leaves on each others. An optimization factor $\alpha$ allows to control the proportion of each of the contributions.

$$BRDF(\theta_i, \theta_o, \varphi_i, \varphi_o) = \alpha R_{vol}(\theta_i, \theta_o, \varphi_i, \varphi_o) + (1 - \alpha) R_{geo}(\theta_i, \theta_o, \varphi_i, \varphi_o) \quad (4)$$

Initially developed for thick and sparse canopies, extensions of this model have been proposed for thin and dense media [17]. [18] brought modifications to the model for the study of snow, and multiple works [16, 19, 20] agree that it can be used for bare soils although they do not explicitly introduce the changes that should be done in such a case. We propose an adaptation of this model for dry powders in Supplementary work. Its unique optimization factor makes this model too rigid and remains unable to reproduce fine appearance variations.

## 3. Background

As a dry powder layer involves a large number of scatterers with relatively high density, it is complex to describe the scattering of light in such a medium. When there is no interference phenomena at stake, it can be done using the RTT [21, 22]. As it relies on the assumption that the scattered fields are uncorrelated, RTT by itself is unable to predict the appearance of materials which exhibits a strong back-scattering contribution.

A deeper analysis of the back-scattering phenomena (see Section 3.5) occurring in a scattering medium refers to Enhanced Back-Scattering (EBS), which describes it as the result of constructive interferences [23, 24]. Contrary to the RTT, the scattered fields are supposed to be correlated. The EBS is therefore the sum of the two terms from RTT (single and multiple scattering) and a third one accounting for the constructive interferences in the backward direction. As the analytical results are deduced within the diffusion approximation, only isotropic particles are considered. Nevertheless, Tsang & Kong [23] extended the model to account for anisotropic scattering in the case of incident field at normal incidence.

This section gives an overview of the theoretical aspects and hypotheses used in the EBS theory.

### 3.1. Medium description

Let's consider a region 0 free of scatterers with wavenumber $k_0 = 2\pi/\lambda$, and a medium 1 characterized by an effective complex wavenumber $K = K' + iK''$ containing scattering spherical inclusions with a density $\mathcal{D}$. All scatterers have the same radius $b$ and scattering properties (cf. Figure 1). Note that for clearer reading, we adopt bold notations for vectors representation ($\vec{v} = \mathbf{v}$) in the following.

The medium is illuminated in the incoming direction $\mathbf{K}_i$ with an angle $\theta_i$ to the surface normal and azimuthal angle $\varphi_i$ by light coming from a point source at the position $\mathbf{r}_i$ in the region 0. The light scattered by the slab is observed at a point in region 0 in the direction $\mathbf{K}_o$, described by the angle $\theta_o$ to the normal, and $\varphi = |\varphi_i - \varphi_o|$ in the azimuthal plane at the position $\mathbf{r}_o$ which is far enough from the interface ($||\mathbf{r}_o|| \gg 1$). To ease the description of the scattering inside the layer, a point scatterer approximation is made: the scattering events occur at the center point of the spherical scatterers.

### 3.2. Foldy's approximation

For a sufficiently small density $\mathcal{D}$ of scatterers, the effective wavenumber $K$ can be expressed using the Foldy's approximation [25]:

$$K \approx \frac{2\pi}{\lambda} + j\lambda\mathcal{D}F \quad (5)$$

where $j^2 = -1$, $F$ is the complex scattering amplitude describing the scattering properties of the particles, and $\lambda$ is the wavelength at which the light propagates.
<figure><img src="image_2.png" /><figcaption>Fig. 1. N-th order multiple scattering in a medium composed of particles. The incident light follows a path going from the particle at $r_1$ to the one at $r_N$ (black path). Interferences occur with a second path that goes from $r_N$ to $r_1$ (gray path). The incoherent part of the resulting field (denoted as $\mathcal{L}_N$ in Section 5.1) contributes to the ladder term whereas its coherent part (denoted as $C_N$ in Section 5.1) accounts for the cyclical term.</figcaption></figure>

### 3.3. Light propagation

The propagation of light inside a medium containing scattering particles is illustrated by the Figure 1. The light wave coming from region 0 without being scattered is described by the incident field:

$$ \psi(r_i) = \exp(jK_i \cdot r_i) \qquad (6) $$

where $K_i$ is of magnitude K given by the Foldy's approximation, and oriented in the incident direction.

After it has been scattered by a particle, light travels in the region 1 towards another scatterer, as if it was in a space free of scatterers. Under this hypothesis, the propagation of light from a position $r'$ towards r is given by the Green's function $G(r, r')$, that satisfies the Helmholtz equation:

$$ (\Delta + k^2)G(r, r') = \delta(r - r') \qquad (7) $$

By averaging over the particles, the mean Green's function referring to the propagation of light from an inclusion to another is given as follows:

$$ G_{11}(r_n, r_m) = \frac{\exp(jK \cdot |r_n - r_m|)}{4\pi|r_n - r_m|} \qquad (8) $$

where $r_m$ and $r_n$ are the positions of two distinct scatterers in the medium 1. In the same way,

$$ G_{01}(r_n, r_o) = \frac{\exp(jk_0 r_o)}{4\pi r_o} \exp(-jK_o \cdot r_n) \qquad (9) $$

describes the propagation from a scatterer at position $r_n$ in region 1 towards the camera at $r_o$ in region 0. The vector $K_o$ is of magnitude K and oriented towards the camera.
### 3.4. Light scattering by one particle

The scattering events caused by the encountering of light with the n-th scatterer is described using the transition operator $T_n$ that we introduced as :

$$T_n(r_n) = 4\pi F(\Theta_n)\delta(r_{in} - r_{on}) = 4\pi F(\Theta_n) \quad (10)$$

where $\Theta_n = \theta_{on} - \theta_{in}$ is the scattering angle between the incident direction on the n-th scatterer $\theta_{in}$ and the outgoing direction $\theta_{on}$.

The complex scattering amplitude $F(\Theta_n)$ is linked to the scattering phase function $p(\Theta_n)$ using the relationship:

$$\rho p(\Theta_n) = \frac{4\pi N}{\kappa} |F(\Theta_n)|^2 \quad (11)$$

where $\rho$ is the single scattering albedo and $\kappa$ the extinction rate.

Although any phase function can be used, a common one would be the Henyey-Greenstein phase function [26] (Eq. 12). It involves one parameter $g$ known as the anisotropy and whose value lies between -1 and 1 describing either backward or forward scattering. Isotropic scattering would then be described with $g = 0$.

$$p(\Theta) = \frac{1}{4\pi} \frac{1 - g^2}{(1 + g^2 - 2g\cos\Theta)^{3/2}} \quad (12)$$

### 3.5. Enhanced back-scattering

The EBS theory offers a physical description of the back scattering phenomenon observed in dense media. It is an extension of RTT [21, 22]: the so called ladder terms accounting for both single scattering $BRDF_{ebs, L_1}$ and some multiple scattering $BRDF_{ebs, L_{2+}}$ take into account the contribution of uncorrelated scattered fields already described by the RTT. This model also accounts for coherent scattered fields that induce interference effects as light undergoes multiple scattering events in the medium. It is described by a cyclical term $BRDF_{ebs, C}$ and results in a sharp peak in the back-scattering direction.

The BRDF is then composed of three contributions. Since the angular and spectral dependency of these three components is the same than the overall BRDF $(\theta_i, \theta_o, \varphi_i, \varphi_o, \lambda)$, it is not reminded in the following equations for readability purposes.

$$BRDF_{ebs}(\theta_i, \theta_o, \varphi_i, \varphi_o, \lambda) = BRDF_{ebs, L_1}(\kappa, g) + BRDF_{ebs, L_{2+}}(\kappa) + BRDF_{ebs, C}(\kappa) \quad (13)$$

This model has been developed for thick layers of a medium composed of spherical scatterers of radius $b$ embedded in a host medium [23, 24, 27, 28].

Analytical expressions of the multiple scattering terms are obtained by considering that for optically thick media with high albedo the resulting overall scattering statistically tends to be isotropic. So the phase function used in multiple scattering is supposed to be constant $p(\Theta) = \frac{1}{4\pi}$. This is known as the diffusion approximation and leads to results in Eqs. 14.

$$BRDF_{ebs, L_1}(\kappa, g) = \frac{\rho p(\Theta)}{\mu_o + \mu_i} \left[ 1 - \exp\left(-\kappa d\left(\frac{1}{\mu_o} + \frac{1}{\mu_i}\right)\right) \right] \quad (14a)$$

$$BRDF_{ebs, L_{2+}}(\kappa) = \frac{3\rho^2}{1 - \mu_i^2 c^2} \left[ \frac{2Q}{1 - \mu_o c} - \frac{\mu_i^2}{\mu_i + \mu_o} \right] \quad (14b)$$

$$BRDF_{ebs, C}(\kappa) = \frac{3\rho^2[c + v(1 - \exp(-2c\tau_o))]}{2\mu_i\mu_o c v[(c + v)^2 + u^2]} \quad (14c)$$
with $\mu_i = \cos\theta_i$, $\mu_o = \cos\theta_o$, d the thickness of the powder layer, and $\tau_o$ the extinction ratio whose expression is given in Eq. 16.

Multiple parameters involved in Eqs 14 are introduced to make the equation more legible and depend on the geometry under which the medium is illuminated and observed. These parameters have been introduced by [24] for isotropic scatterers such as:

$$ u = \frac{k_0}{\kappa}(\mu_i - \mu_o) \qquad (15a) $$

$$ v = \frac{1}{2}\left(\frac{1}{\mu_i} + \frac{1}{\mu_o}\right) \qquad (15b) $$

$$ \alpha = \frac{k_0}{\kappa} f \qquad (15c) $$

$$ c = \sqrt{3(1 - \rho) + \alpha^2} \qquad (15d) $$

$$ Q = \frac{[1 + \mu_i c - (1 - \mu_i c)\exp(-2c\tau_o(\lambda))]}{4c} \qquad (15e) $$

where f is the norm of the vector f defined as difference between the projections of the normalized incoming and outgoing direction vectors (Figure 2) which allows the BRDF modelization in any azimuthal plane.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Definition of the vector f.</figcaption></figure>

Only the single scattering term involves the phase function in such a case, and therefore the anisotropy parameter g if the Henyey-Greenstein phase function is at use (see Eq. 12). Tsang & Kong [23] introduced the consideration of anisotropic scatterers but only with normal incidence. We will discuss the extension of the model for anisotropic scattering in the Section 5. The last parameter involved in the BRDF expressions (Eqs. 14) is the extinction ratio $\tau_o$ that is defined as intrinsic to the scatterers' material. It accounts for internal reflections that may occur between the interface of the medium. A definition is given by [29]:

$$ \tau_o(\lambda) = \frac{2(1 - \bar{R}(\lambda))}{3(1 + \bar{R}(\lambda))} \qquad (16) $$

$$ \text{with } \bar{R}(\lambda) = \frac{3C_2 + 2C_1}{3C_2 - 2C_1 + 2} \quad \text{and} \quad C_i = \int_0^{\frac{\pi}{2}} R(\theta, \lambda) \cos^i \theta \sin \theta d\theta $$

In the previous equations, $R(\theta)$ is the Fresnel's angular reflection coefficient [30]. The parameter $\tau_o$ thereby depends on the optical constants of the material the particles are made of and on the wavelength.

## 4. Sample description and BRDF measurements

The macroscopic physical phenomena at stake in the appearance of a material can be inferred through BRDF measurements.
242

## 4.1. Sample description

243

244

245

246

For the study of dry powder appearance, we focused on synthetic nano-crystallized hematite $\alpha$-Fe$_2$O$_3$ powders with various colors, that we obtained using different routes allowing to control the nucleation and growth of $\alpha$-Fe$_2$O$_3$ crystallites in solution. These samples are presented in Figure 3.

247

248

249

250

251

252

The samples denoted from 1 to 5 have been prepared following methods described by Schwertmann & Cornell [31]. 1, 2 and 5 have been obtained by forced hydrolysis of Fe$^{III}$ salts solutions in acidic conditions, while 3 and 4 are the results of the transformation of a ferrihydrite precipitate in aqueous suspension under basic conditions. We generated peanut-like hematite particles 6 following the routes proposed by Sugimoto et al. [32], NaH$_2$PO$_4$ salts. Sample 7 is a commercial Iron(III) oxide powder sample purchased from Puratronic®.

253

254

255

256

257

258

259

260

The color and the morphometry of each sample are respectively deduced from naked eye observations completed with optical microscopy using a Nikon Eclipse LV100ND microscope, and by scanning electron microscopy observations performed with a Zeiss Ultra+ microscope. Various grain shapes and sizes ranging from 10 nm to 10 µm are obtained. The small grains are characterized by a lighter red color (samples 1 to 5) whereas a darker purple hue is observed for the pigments with larger grains (samples 6 and 7). We may also note that samples 5, 6 and 7 can be described by two grain sizes: either two different size populations (sample 5), or large grains that are agglomerates of smaller ones (samples 6 and 7).

<figure><img src="image_4.png" /><figcaption>Fig. 3. $\alpha$-Fe$_2$O$_3$ hematite samples under study, showing various colors and grains morphology. Upper row: The photographs are calibrated using an X-rite 24 patches Classic Colorchecker. Middle row: Observations using optical microscopy. Lower row: Samples morphology observations using Scanning Electron Microscopy. The samples 1 to 4 show similar grain shape, with slightly different sizes from 10 to 50 nm in diameter, and light red hue. Sample 5 presents two size populations and shows a red-pinkish hue. Samples 6 and 7 have a more complex morphology with grains that seem to be agglomerations of smaller ones. The grains are larger and a darker purple hue is observed. According to optical microscopy, 7 is even composed of different types of grains, with different colors.</figcaption></figure>
## 4.2. BRDF measurements

The whole BRDF of each sample is measured by combining two distinct measurements, performed with two different setups depending on which scattering geometry is at stake. The first one does not account for the back-scattering, while the second one is specifically performed to account for this contribution. All the measurements are done on $d = 1$ mm thick powder layers.

### 4.2.1. Measurement out of the back-scattering direction

The main part of the BRDF is measured using the spectro-goniophotometer SHADOWS [33], for over 180 geometries: the scattered light is observed for $\theta_o \in [-70°, 70°]$ with a step of $10°$, for incidences $\theta_i \in [0°, 20°, 40°, 60°]$ and azimuthal angles $\varphi \in [0°, 30°, 60°, 90°, 135°]$. The spectral range has been limited to 360 nm - 860 nm with a step of 20 nm. The reflectance measurements are calibrated with a Spectralon. Every point is measured with an integration time of 300 ms and averaged over 10 spectra. Such measurement takes about 18 hours.

When it comes to a geometry where incident and observation directions are identical, the instrument design causes the sensor to be obstructed by the light source. A specific device needs to be designed to evaluate the amount of light scattered in these directions.

### 4.2.2. Measurement in the back-scattering direction

Our design for measuring the back-scattering contribution to the BRDF is inspired from the work of [9] and [11]. As these only deal with vertically held samples, they cannot be used for studying dry powders. We propose a new design to allow the measurement of any horizontally held sample.

We chose a halogen light source to illuminate at normal incidence the sample after being reflected by a beam splitter at $45°$. The size of the illumination spot is minimized by focusing the light with a 50 mm focal lens. An additional diaphragm is added to clean the focused beam. The direction of the incident light $\theta_i$ is adjusted by positioning the beam splitter at an angle of $45° + \theta_i/2$. As the light interacts with the powder, it is scattered in every direction and preferably back into the incident direction. It is then partially transmitted through the beam splitter and collected by a Thorlabs F950FC-A collimator. The latter is fixed on a vertical rotating plate so that we can scan different observation directions around the exact back-scattering one in the incident plane. Although the incident light may be partially polarized after reflecting on the beam splitter, for small enough incident angles the diffuse interaction with the sample has the property of depolarizing the light. For large incidence angles, the reflection on the surface becomes more important and so are the polarization effects. To neglect these polarization effects, we will only consider incident angles smaller than $40°$. We use a fiber spectrometer OceanOptics USB650 along with the SpectraSuite software in order to detect and visualize the measured signal.

In the incident plane ($\varphi = 0°$), we scan the observation angles in the range [$\theta_i - 10°, \theta_i + 10°$] with a step of $2°$, for incidences $0°, 20°$ and $40°$. For each incidence angle of interest, the observation angles the furthest from the exact back scattering direction that are scanned are redundant with geometries reached by the SHADOWS instrument. This allows for verification that the two distinct measurements are consistent with one another. The reflectance measurements are calibrated with a Spectralon. The integration time is chosen to maximize the detected signal and an averaging over 5 spectra is performed to minimize the noise. Less than 20 seconds is necessary to measure one spectrum.

## 5. Powder BRDF Model

As shown by the measurements on Figure 6, powders' BRDFs mainly involve back-scattering and almost no specular scattering which makes the usual specular BRDF models ineffective to reproduce their appearance.
<figure><img src="image_5.png" /><figcaption>Fig. 4. Back-scattering measurement device. The beam of a white source is focused on the sample surface at an incidence $\theta_i$ using a 50 mm focal lens L and a beam splitter. The reflected light is collected through a collimator C in different observation directions $\theta_o$.</figcaption></figure>

Although our model is derived from a study on pure nano-crystallized hematite powders, we chose a physically-based approach so that it can be applied to any other pigment. It involves both chemical composition and morphological parameters to accurately account for the color variations observed on such pigments. Moreover it allows to consider either isotropic or anisotropic scatterers. Inspired by both RTLSR and Spectralon's models, our model is divided into two contributions, accounting for volume and surface scattering respectively.

$$BRDF(\theta_i, \theta_o, \varphi_i, \varphi_o, \lambda) = A_vBRDF_v(\theta_i, \theta_o, \varphi_i, \varphi_o, \lambda) + A_sBRDF_s(\theta_i, \theta_o, \varphi_i, \varphi_o, \lambda) \quad (17)$$

These two contributions are weighted with the scalar amplitude parameters $A_v$ and $A_s$. Because the angular and spectral dependency of the two contributions $BRDF_v$ and $BRDF_s$ is the same as the complete BRDF, it will not be reminded in the following equations for readability purposes. The main difference with the previous models is that, based on the EBS theory the back scattering is considered as part of the volume contribution. The surface contribution is then reduced to specular scattering caused by surface roughness.

## 5.1. Volume scattering

As EBS theory is established within the diffusion approximation, the average scattering is assumed to be isotropic. This greatly simplifies the multiple scattering contribution by considering the scattering phase function as constant over the whole scattering hemisphere. The anisotropic scattering of each particle is accounted for by choosing the phase function parameter g involved in the single scattering contribution. The extension of the model for taking into account this anisotropy in multiple scattering has only been done for a normal incidence by Tsang et al. [23].

In the following we expose the calculations we made in order to consider anisotropic scattering in both single and multiple scattering terms for any $\theta_i/(\theta_o, \varphi_o)$ geometry. The final model only involves two variable parameters (phase function anisotropy parameter g, extinction rate $\kappa$) that rely on the physical and chemical properties of the material that is under study.

**First order ladder term $L_1$**: This contribution accounts for light that is only scattered once before being observed and is similar to the one introduced in Section 3.5. The anisotropy is taken care of by the scattering phase function involved in the analytical expression of $BRDF_{ebs, L_1}$ in Eq. 14.
**Multiple ladder term $\mathcal{L}_{2+}$:** The phase function is involved at every scattering event that light encounters as it travels the particulate medium. We denote $P_1$ and $P_n$ respectively the first and the last particles that are responsible for light scattering along a light path.

At the second scattering order, the light is first incoming on $P_1$, scattered with an angle $\Theta_1$ towards the particle $P_n$, and finally scattered towards the observer with an angle $\Theta_n$ as illustrated by the Figure 5.

For the n-th scattering order, light undergoes n scattering events. As for the second order, it undergoes a first one described by the angle $\Theta_1$ while encountering the particle $P_1$ and last one with the particle $P_n$ and the angle $\Theta_n$. But it also interacts with n-2 other particles in between, resulting in n-2 scattering events respectively described by the scattering angles $\Theta_k = \Theta_{ik} - \Theta_{ok}$, difference between the incoming direction $\theta_{ik}$ and outgoing direction $\theta_{ok}$ on the k-th particle, with $k = \ll 2; n-1 \gg$ as it is shown on the Figure 5.

<figure><img src="image_6.png" /><figcaption>Fig. 5. Multiple scattering with two and n particles. The first scatterer that light encounters as it travels inside the medium is denoted as $P_1$ and involves a scattering event described by the angle $\Theta_1$. The last scattering event before light goes out of the medium is caused by the scatterer $P_n$ with an angle $\Theta_n$. In between these two scatterers, light is scattered with an angle $\Theta_k$ at each interaction with a particle k.</figcaption></figure>

By averaging over all the possible paths that light can follow inside the medium while encountering n scatterers, we can write the n-th order ladder term as follows:

$$ \mathcal{L}_n = \frac{(4\pi)^{n-2}(\rho\kappa)^n}{r^2} \times \int |p(\Theta_n)|^2 |G_{11}(\mathbf{r}_n, \mathbf{r}_{n-1})|^2 \dots |p(\Theta_1)|^2 \exp\left(\kappa\left(\frac{z_1}{\mu_o} + \frac{z_n}{\mu_i}\right)\right) d\mathbf{r}_1 \dots d\mathbf{r}_n \quad (18) $$

with r the distance from the sample to the observation point.

The total ladder term accounting for multiple scattering $k = \ll 2; \infty \gg$ can be written in a recursive way, giving:

$$ \mathcal{L}_{2+} = \iint \exp\left(\kappa\left(\frac{z_i}{\mu_o} + \frac{z_o}{\mu_i}\right)\right) A(\mathbf{r}_i, \mathbf{r}_o) d\mathbf{r}_i d\mathbf{r}_o \quad (19) $$

with $A(\mathbf{r}_i, \mathbf{r}_o)$ defined such as:

$$ \begin{aligned} A(\mathbf{r}_i, \mathbf{r}_o) = & \left(\frac{\rho\kappa}{4\pi}\right)^2 p(\Theta_1)p(\Theta_n) \frac{\exp(\kappa|\mathbf{r}_i - \mathbf{r}_o|)}{|\mathbf{r}_i - \mathbf{r}_o|^2} \\ & + \int \left(\frac{\rho\kappa}{4\pi}\right) p(\Theta_k) \frac{\exp(\kappa|\mathbf{r}_i - \mathbf{r}_k|)}{|\mathbf{r}_i - \mathbf{r}_k|^2} A(\mathbf{r}_k, \mathbf{r}_o) d\mathbf{r}_k \end{aligned} \quad (20) $$
Under the assumption of the diffusion approximation, this term can be simplified to average the contributions of all scatterers by extending the work from [23]. The resulting BRDF of the ladder contribution is therefore reduced to the following analytical expression:

$$BRDF_{L_2+}(\kappa, g) = \frac{3\rho^2}{1 - \mu_i^2 c^2} \left[ \frac{2Q}{1 - \mu_o c} - \frac{\mu_i^2}{\mu_i + \mu_o} \right] \quad (21)$$

with

$$c = \sqrt{3(1 - \rho)(1 - \rho g) + \alpha^2} \quad (22)$$

and parameters $\alpha$ and $Q$ defined as in Eqs 15.

One should note that the ladder term expression is similar to the one defined by the Eq. 14, except that it accounts for possible anisotropic average scattering in the medium through a varying phase function of parameter $g$ (Eq. 22). If one wants to describe isotropic scattering, a constant phase function such as $p(\Theta) = p$ (similarly $g = 0$) should be used which leads back to the expression originally introduced in Section 3.5.

**Multiple cyclical scattering C:** As for the ladder term, we are aiming to take into account the anisotropic scattering of particles in the cyclical term. We will use the same assumptions as illustrated by the Figure 5. We consider a beam of light that is scattered by $n$ particles as it follows a path inside the medium. By accounting for the correlation between this path and its reversal and averaging over all the possible paths, the $n$-th order cyclical term can be written:

$$C_n = \frac{(4\pi)^{n-2}(\rho\kappa)^n}{r^2} \int |p(\Theta_n)|^2 |G_{11}(\mathbf{r}_n, \mathbf{r}_{n-1})|^2 \dots |G_{11}(\mathbf{r}_2, \mathbf{r}_1)|^2 |p(\Theta_1)|^2 \\ \times \cos \left( k(\mu_i - \mu_o)(z_1 - z_n) \right) \exp \left( \frac{\kappa}{2} \left( \frac{1}{\mu_o} + \frac{1}{\mu_i} \right) (z_1 + z_n) \right) dr_1 \dots dr_n \quad (23)$$

The complete cyclical contribution is obtained by summing over the $n$-th orders and can be simplified replacing $G_{11}$ by its definition (Eq. 8). By using the new expression of $A(\mathbf{r}_i, \mathbf{r}_o)$ given by Eq. 20:

$$C = \iint \cos(2k(\mu_i - \mu_o)(z_o - z_i)) \exp\left(\frac{\kappa}{2}\left(\frac{z_i}{\mu_o} + \frac{z_o}{\mu_i}\right)\right) A(\mathbf{r}_i, \mathbf{r}_o) d\mathbf{r}_i d\mathbf{r}_o \quad (24)$$

Finally, by using the diffusion approximation, the cyclical term contributes to the complete appearance of the dry powder through the partial BRDF:

$$BRDF_C(\kappa, g) = \frac{3\rho^2[c + v(1 - \exp(-2c\tau_o))]}{2\mu_i\mu_o c v[(c + v)^2 + u^2]} \quad (25)$$

where $c$ is defined as in Eq. 22, and $u, v$ as in Eqs. 15.

Once again, this expression is similar to the one obtained assuming isotropic averaged scattering in the medium (Eq. 14). However it involves the phase function parameter $g$ to account for anisotropic scattering. The use of a constant phase function $p(\Theta) = p$ (similarly $g = 0$) brings back to the isotropic case already discussed in Section 3.5.

Our model's complete EBS contribution accounting for anisotropic scattering in the powder medium is the sum of the two ladder components given by first equation of Eq. 14 and Eq. 21, and the cyclical component given by Eq. 25.
## 5.2. Surface scattering

In our model, surface scattering is reduced to the specular contribution. The light that is incoming on a rough surface such as the interface between air and a dry powder is scattered in a lobe around the specular direction.

Two options have been considered to include such a phenomenon in our model. The former is to model the specular peak by a gaussian function. According to the measurements performed in Section 4 the specular scattering is very weak so this simple method is reliable enough to reproduce the specular contribution of dry powders, although it is empirical. By assuming a uniform distribution of the surface normals, it requires only one parameter which is the gaussian standard deviation $\sigma$. With this modeling the surface scattering BRDF is written as:

$$BRDF_s(\sigma) = \frac{1}{\sigma\sqrt{2\pi}}\exp\left(-\frac{\Theta^2}{2\sigma^2}\right) \qquad (26)$$

The second option, which is the one we chose to use, relies on the physical approach from Harvey-Shack [34] that accounts for diffraction effects caused by the surface asperities. Using this theory, the incoming light on a rough surface is then scattered around the specular direction according to the following BRDF:

$$BRDF_s(l_c) = \frac{4\pi^2}{\lambda^4}(\cos \theta_i + \cos \theta_o)^2 \text{PSD}(\mathbf{f}) \qquad (27)$$

where the power spread function (PSD) models a rough surface topography. We chose a gaussian power spread function, involving a roughness parameter $\sigma_s$ and a correlation length $l_c$ that quantifies the mean distance for which the surface roughness of two points on the surface are completely independent.

$$\text{PSD}(\mathbf{f}) = \pi \sigma_s^2 l_c^2 \exp(-\pi^2 l_c^2 (f_x^2 + f_y^2)) \qquad (28)$$

## 5.3. Complete model

Our complete model is the weighted sum of volume and surface contributions according to Eq. 17 respectively given by the adapted enhanced back scattering model and a specular model. Note that the roughness parameter involved in the surface contribution has been included in the amplitude parameter $A'_s = A_s \sigma_s^2$ in order to reduce the number of optimization parameters.

$$BRDF(\theta_i, \theta_o, \varphi_i, \varphi_o, \lambda) = A_v[BRDF_{L_1}(\kappa, g) + BRDF_{L_2}(\kappa, g) + BRDF_C(\kappa, g)] + A'_s BRDF_s(l_c) \qquad (29)$$

The model depends on a total of 2 empirical ($A_v$ and $A_s$) and 3 physical ($g, \kappa, l_c$) parameters to accurately simulate the BRDF of dry powder. As it is shown in the Section 6, our model works even better while including another physical parameter to refine which is the *effective radius* of the grains, defined as the mean size of the grains actually seen by the incoming light and therefore responsible for the scattering. The influence of such parameter is shown on Figure 7. The parameters involved in our model are summarized in Table 1.

## 5.4. Optimization procedure

In order to reproduce the BRDF of a powder sample, the model is fitted to measurements by optimizing the variable parameters of Table 1. We used the least squares minimization of scipy along with the sequential least square programming method (SLSQP) to fit our BRDF model to the measurement data $y_m$ (Eq. 30). A first option is to optimize while using a fixed radius.
<table><caption>Table 1. Parameters involved in our BRDF model. The model can be optimized by refining only 5 parameters while the grain radius *b* is fixed, or by refining 6 parameters including the same 5 previous ones and an effective radius.</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>Volume scattering</td><td></td></tr><tr><td>A<sub>v</sub></td><td>Intensity of volume scattering</td></tr><tr><td>g</td><td>Scattering anisotropy</td></tr><tr><td>κ</td><td>Extinction rate</td></tr><tr><td>b</td><td>effective radius (optional)</td></tr><tr><td>Surface scattering</td><td></td></tr><tr><td>A'<sub>s</sub></td><td>Intensity of surface scattering</td></tr><tr><td>l<sub>c</sub></td><td>Correlation length</td></tr></tbody></table>

The average dimensions of the grains of a powder sample is deduced from SEM images and reduced to only one dimensional parameter by computing the radius of a sphere that is equivalent to the grain. To better account for the shape anisotropy this radius is computed so that the ratio of the surface over the volume of the equivalent sphere and of the grain remains constant. The final mean radius is obtained by averaging on the radii of multiple grains that are in the field of view of a single SEM image. A second option is to optimize on the grains' effective radius, which both releases constraints on the optimization process and gets rid of the uncertainty on the experimental grain size definition.

The optimization is performed for all measurement geometries $\gamma = (\theta_i, \theta_o, \varphi_i, \varphi_o)$ with $\theta_i < 60^{\circ}$.

For larger incidences, we were not able to measure the back-scattering contribution so we do not account for these geometries in the optimization process.

$$E = \min\left( \sum_{\gamma, \theta_i \le 40^{\circ}} |\text{BRDF}_{\gamma}(A_v, g, \kappa, b, A'_s, l_c) - y_{\gamma,m}|^2 \right) \quad (30)$$

We only account for the measurements at the wavelength 700nm. The spectral dependency of the model will be discussed later on.

# 6. Results

## 6.1. Complete measurement

The complete BRDF measurements of dry hematite $\alpha - Fe_2O_3$ powders were performed in two independent steps described in Sections 4.2. By comparing the results for geometries reached by both instruments, the measurements happen to be shifted in amplitude from one another. This may be caused by the underestimation of the back scattering of the spectralon used for photometric calibration. Since the difference is smaller than 5%, a re-scaling of the back-scattering measurement is done so it matches the goniometer results.

A typical experimental BRDF is shown on Figure 6 after merging the measured data. The measurements show that the BRDF of dry hematite powders is composed of an overall diffuse contribution, and two directional scattering phenomena including mainly back-scattering and some specular scattering (Figure 6). The BRDFs of characterized samples vary in intensity and width of the specularly and back scattered lobes. We give results for the six other samples under study in Supplementary Work.
<figure><img src="image_7.png" /><figcaption>Fig. 6. Complete BRDF measurement of the hematite powder sample 1 at the wavelength 700 nm for an incidence $\theta_i$ = 40°. The measured geometries are spotted with black dots, and the intermediate values are interpolated. Light is mostly reflected in the backward direction, and a very weak specular scattering peak may be observed as well. The spreading of the back scattering peak with the azimuth angle is due to the colormap interpolation between the measured points.</figcaption></figure>

## 6.2. Model validation

The ability of the model to simulate the BRDF of pigments is validated against the various hematite powders measurements by following the procedure of Section 5.4. In the following we chose to define the grains' effective radius as an optimization parameter as it leads to better results. Indeed, as it is shown in Table 2 it allows a nicer modeling of the back scattering peak amplitude, while keeping values pretty close to the radius deduced from microscopy observations. Although the two populations are well distinguishable in samples 5 and 7 (see Figure 3), it is harder to spot on sample 6. In that case, the smaller grains actually refer to the roughness at the surface of the observed ellipsoidal grains. For illustration, the results from two optimizations involving or not the radius performed on sample 1 are shown on Figure 7.

The same optimization process is performed using both Spectralon and RTLSR models adapted to be relevant in the case of dry powders. Such adaptations to the original models are detailed in Supplementary work. To compare our results to the one obtained through these two other models allows a better evaluation of our model.

Our model is able to reproduce the measured BRDF of hematite powders with an accuracy comparable to that of the Spectralon model (< 8% of the measured reflectance factor), but with a much lower complexity: only 6 parameters related to physical properties of the samples are required against 12 empirical ones for the Spectralon's model (cf. Figures 9, 10 and 11). Our model nicely fits the back-scattering peak, and is even more accurate than the other ones for reproducing the slope of the measured BRDF (cf. Figure 8). The mean squared error (MSE) from the fitting procedure is also of same magnitude for both Spectralon and our model (MSE < 0.05). In comparison, the RTLSR gives much poorer results (MSE > 0.1). Detailed values are given for each samples in Table 2. Our model is particularly efficient in the case of samples whose morphology is easily characterized, that is to say for the samples with well-crystallized single-crystals. Both the back-scattering and the specular contributions are detected and modeled,
even though the latter is sometimes too weak for the modeling to be reliable. Errors are mostly spotted around the back scattering peak because of the insufficient measurement resolution and interpolation at these geometries.

One should note that we have encountered difficulties in the modelization of samples 5 and 7's BRDFs that we attributed to the more complex morphology of these powders. These powders are not monodisperse since they show two size populations and agglomerations of grains. An extension of our model is proposed in Supplementary work to deal with powders presenting two grain size populations.

<table><caption>Table 2. Results from optimizations on radii. The radius deduced from SEM observations are compared with the effective radius obtained by the optimization of our model. For samples 5 to 7, the model extension for two size populations presented in Supplementary work is used. As it takes two different radii to optimize, two effective radii are obtained. The resulting MSE of such optimizations are also given and can be compared to MSE of optimizations over the Spectralon model and the RTLSR model.</caption><thead><tr><th>Sample</th><th>Radius SEM (nm)</th><th>Effective radius (nm)</th><th></th><th>Ours MSE</th><th>Spectralon MSE</th><th>RTLSR MSE</th></tr></thead><tbody><tr><td>1</td><td>12 ± 1.6</td><td>16</td><td></td><td>0.035</td><td>0.043</td><td>0.17</td></tr><tr><td>2</td><td>13 ± 3.0</td><td>17</td><td></td><td>0.052</td><td>0.052</td><td>0.21</td></tr><tr><td>3</td><td>22 ± 3.7</td><td>16</td><td></td><td>0.044</td><td>0.046</td><td>0.25</td></tr><tr><td>4</td><td>25 ± 4.1</td><td>15</td><td></td><td>0.037</td><td>0.039</td><td>0.29</td></tr></tbody><thead><tr><th>Sample</th><th>Large radius SEM (nm)</th><th>Effective radius (nm)</th><th>Small radius SEM (nm)</th><th>Effective radius (nm)</th><th></th><th></th><th></th></tr></thead><tbody><tr><td>5</td><td>83 ± 34</td><td>83</td><td>47 ± 5.3</td><td>46</td><td>0.063</td><td>0.074</td><td>0.71</td></tr><tr><td>6</td><td>687 ± 64.0</td><td>687</td><td>10 ± 0.4</td><td>43</td><td>0.041</td><td>0.043</td><td>0.49</td></tr><tr><td>7</td><td>2222 ± 1054</td><td>2222</td><td>49 ± 13</td><td>49</td><td>0.019</td><td>0.019</td><td>0.30</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Influence of the grain radius optimization on the results of sample 2, in the incident plane ($\varphi = 0^\circ$) for incidences $\theta_i = 0^\circ, 20^\circ$ and $40^\circ$. The effective radius found by optimization ($b = 17$ nm) is pretty close to the theoretical result deduced from SEM observations ($b = 13 \pm 3.0$ nm).</figcaption></figure>

Since the refined parameters permit our model to fit quite precisely with the measurements in any geometry, it appears that none of them depends on the incoming, outgoing and azimuthal
<figure><img src="image_9.png" /><figcaption>Fig. 8. Comparison of the optimization results of sample 2 from the three BRDF models: Our model, Spectralon model, RTLSR model. (a) Results in the incident plane ($\varphi = 0°$) for incidence $\theta_i = 20°$. (b) Residual error between the measurements and the fits. Our model slightly underestimates the back-scattering peak compare to the Spectralon model. However, it fits the slope of the measured BRDF much better. The residual error of these two models are similarly small. In comparison, the RTLSR model is less suitable to fit the measured data.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Upper row: Fitting of the measured BRDF of sample 1 using our model, spectralon and RTLSR models adapted to dry powders, at incidence $\theta_i = 0°$. Lower row: difference between the measurement and the fit. Our physically-based model (6 optimized parameters) gives better results than the two other ones as it is more performant than the RTLSR model (1 optimized parameter) and as performant as the Spectralon model (12 optimized parameters) with far less optimization parameters.</figcaption></figure>

<footnote n="480">480 directions. The model is then able to predict the appearance of the studied powder in any other geometry that hasn't been considered during its optimization. The results are convenient even for more grazing angles such as incidence $\theta_i = 60°$ as shown on Figure 7.</footnote>
<figure><img src="image_11.png" /><figcaption>Fig. 10. Upper row: Fitting of the measured BRDF of sample 1 using our model, spectralon and RTLSR models adapted to dry powders, at incidence $\theta_i = 20^\circ$. Lower row: difference between the measurement and the fit. Our physically-based model (6 optimized parameters) gives better results than the two other ones as it is more performant than the RTLSR model (1 optimized parameter) and as performant as the Spectralon model (12 optimized parameters) with far less optimization parameters.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Upper row: Fitting of the measured BRDF of sample 1 using our model, spectralon and RTLSR models adapted to dry powders, at incidence $\theta_i = 40^\circ$. Lower row: difference between the measurement and the fit. Our physically-based model (6 optimized parameters) gives better results than the two other ones as it is more performant than the RTLSR model (1 optimized parameter) and as performant as the Spectralon model (12 optimized parameters) with far less optimization parameters.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 12. Optimization of our model on sample 4, at different wavelengths using the optical constants from [35]. For smaller wavelengths, our model is not able to accurately fit the measurements.</figcaption></figure>

## 7. Discussion & Conclusion

From BRDF measurements, we highlighted that in the particular case of dry powders, the light is mainly back scattered. We proposed a new BRDF model accounting for three different scattering interactions: single scattering, incoherent multiple scattering and coherent multiple scattering responsible for the sharp peak in the backward direction. It is a semi-empirical and analytical model based on the propagation of light inside a granular medium that involves the physico-chemical properties of the powders connected to only six refinement parameters.

The model has been validated against the BRDF measurements of seven $\alpha$-Fe$_2$O$_3$ hematite powders, and compared to two adapted existing models. Ours is able to reproduce the measured BRDF with an accuracy comparable to the one of the Spectralon's model, but with a much lower complexity. It can be used in its original form to reproduce the appearance of a mono disperse powder with spherical grains, or in its extended form presented in Supplementary work to account for more complex morphology involving multiple grains size populations.

Our model still has some limitations. The results obtained for specular scattering are quite questionable since the specular contribution is very weak and almost unnoticeable on the BRDF measurements. A more accurate model should be established by studying powders showing shiny and sparkling reflections such as eye shadow.

For now, our model can only be used as a fitted model. The parameters describing the appearance of a specific powder are retrieved from the fitting of the model to the measured BRDF of such powder. In order to predict the BRDF of any other powder of known chemical composition and morphology that has not been measured, it would be necessary, as a future work, to establish a relationship between the optimization parameters of the model and the physico-chemical parameters of the material the powder is composed of. This requires to study much more samples with more variability in sizes.

The validity of the model wavelength dependency is also a main issue that couldn't be solved throughout this work. At smaller wavelengths, the optimization is not robust enough to lead to nice fitting results (see Figure 12). This might be caused by a poor knowledge of the optical constants of the studied material which is the only fixed parameter known as spectrally varying. Indeed, it is mainly described by the optical constant at stake whose value appears to be quite doubtful according to the discrepancies between the results in the literature [35, 36]. Such uncertainties are often ignored in both material science and computer graphics whereas it surely affects the results.

As a future work we would like to implement our dry powder BRDF model in a physically based renderer so that photo-realistic renderings of the powders can be made. This would first allow the validation of the model by comparing the rendered image with photographs of the actual
520

521

samples. It would also permit, by letting the user choose the values of the variable parameters to either reproduce the appearance of other powder materials, or to create new ones.

522

523

524

525

526

Although the model has been tested for only one material, we believe it can be used for any other pure pigment. In the case of parietal paintings, ochers (mixture of one iron-based coloring pigment and colorless powders such as clays) and mixtures of different pigments were more commonly used than pure and almost pure iron oxide. In order to reproduce photo-realistic appearance of such paintings, we would like to extend our model to more complex powders.

527

528

529

530

The assumption of optically thick layer of powder is also too restrictive. An extension of the enhanced back scattering theory for thin films and isotropic scatterers has been proposed by [24] and may be a good starting point for our model adaptation, leading to a better reproduction of the appearance of painting layers.

531

532

533

534

535

536

537

538

539

**Acknowledgments:** The authors would like to thank Alain Ibanez and Simon Giraud (Institut Néel) for the hematite samples syntheses. We also thank Corinne Félix (Institut Néel) for her help in developing the back-scattering measurement device, as well as Bernard Schmitt and Olivier Brissaud (IPAG) for their help in the reflectance measurements using SHADOWS gonio-spectro-photometer. At last, we thank Romain Pacanowski (INRIA) for the proofreading of this paper. This work was financially supported by the French National Research Agency in the framework of the Investissements d'Avenir program (ANR-15-IDEX-02, Cross Disciplinary Program Patrimalp).

541

543

544

545

**Disclosures:** The authors declare no potential conflicts of interest.

546

**Data availability:** Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

**Supplemental document:** See Supplement 1 for supporting content.

## References

547

548

1. C. Vignaud, H. Salomon, E. Chalmin, J.-M. Geneste, and M. Menu, "Le groupe des «bisons adossés» de lascaux. étude de la technique de l'artiste par analyse des pigments," *L'anthropologie* **110**, 482-499 (2006).

549

2. A. Casoli, "Research on the organic binders in archaeological wall paintings," *Appl. Sci.* **11**, 9179 (2021).

550

3. L. M. Brady, "Documenting and analyzing rock, paintings from torres strait, ne australia, with digital photography and computer image enhancement," *J. Field Archaeol.* **31**, 363-379 (2006).

551

552

4. I. Domingo, B. Carrión, S. Blanco, and J. L. Lerma, "Evaluating conventional and advanced visible image enhancement solutions to produce digital tracings at el carche rock art shelter," *Digit. Appl. Archaeol. Cult. Herit.* **2**, 79-88 (2015).

553

554

555

556

5. E. Chalmin, M. Menu, M.-P. Pomiès, C. Vignaud, N. Aujoulat, and J.-M. Geneste, "Les blasons de lascaux," *L'anthropologie* **108**, 571-592 (2004).

557

558

6. M. Gerardin, N. Holzschuch, A. Ibanez, B. Schmitt, and P. Martinetto, "Influence of micro-structural features on the colour of nanocrystallised powders of hematite and visible-nir reflectance spectra simulations," *J. international colour association* **26**, 41-48 (2021).

559

560

7. C. F. Bohren and D. R. Huffman, *Absorption and scattering of light by small particles* (John Wiley & Sons, 2008).

561

8. J. R. Frisvad, N. J. Christensen, and H. W. Jensen, "Computing the scattering properties of participating media using lorenz-mie theory," in *ACM SIGGRAPH 2007 papers*, (2007), pp. 60-es.

562

563

9. L. Belcour, R. Pacanowski, M. Delahaie, A. Laville-Geay, and L. Eupherte, "Bidirectional reflectance distribution function measurements and analysis of retroreflective materials," *JOSA A* **31**, 2561-2572 (2014).

564

10. P. Beck, A. Pommerol, N. Thomas, B. Schmitt, F. Moynier, and J.-A. Barrat, "Photometry of meteorites," *Icarus* **218**, 364-377 (2012).

565

566

11. M. Lévesque and M. Dissanka, "Measurement and modeling of the spectralon spectro-polarimetric bidirectional reflectance distribution function (brdf)," *Tech. rep.*, Defence Research and Development Canada (2016).

567

568

12. C. Walthall, J. Norman, J. Welles, G. Campbell, and B. Blad, "Simple equation to approximate the bidirectional reflectance from vegetative canopies and bare soil surfaces," *Appl. Opt.* **24**, 383-387 (1985).

569

570

13. J.-L. Roujean, M. Leroy, and P.-Y. Deschamps, "A bidirectional reflectance model of the earth's surface for the correction of remote sensing data," *J. Geophys. Res. Atmospheres* **97**, 20455-20468 (1992).

571

572

14. W. Lucht, C. B. Schaaf, and A. H. Strahler, "An algorithm for the retrieval of albedo from space using semiempirical brdf models," *IEEE Transactions on Geosci. Remote. sensing* **38**, 977-998 (2000).

573
575

15. J. Ross, "The radiation regime and architecture of plant stands. the hague: Dr. w," Junk Publ. doi **10**, 978-94 (1981).

576

16. X. Li and A. H. Strahler, "Geometric-optical bidirectional reflectance modeling of a conifer forest canopy," IEEE Transactions on Geosci. Remote. Sens. pp. 906-919 (1986).

577

17. W. Wanner, X. Li, and A. Strahler, "On the derivation of kernels for kernel-driven models of bidirectional reflectance," J. Geophys. Res. Atmospheres **100**, 21077-21089 (1995).

578

579

18. Z. Jiao, A. Ding, A. Kokhanovsky, C. Schaaf, F.-M. Bréon, Y. Dong, Z. Wang, Y. Liu, X. Zhang, S. Yin et al., "Development of a snow kernel to better model the anisotropic reflectance of pure snow in a kernel-driven brdf model framework," Remote. sensing environment **221**, 198-209 (2019).

580

581

582

19. T. Nilson and A. Kuusk, "A reflectance model for the homogeneous plant canopy and its inversion," Remote. Sens. Environ. **27**, 157-167 (1989).

583

584

20. W. Ni, X. Li, C. E. Woodcock, M. R. Caetano, and A. H. Strahler, "An analytical hybrid gort model for bidirectional reflectance over discontinuous plant canopies," IEEE Transactions on Geosci. Remote. Sens. **37**, 987-999 (1999).

585

586

21. S. Chandrasekhar, *Radiative transfer* (Courier Corporation, 2013).

587

22. A. Ishimaru, Wave propagation and scattering in random media, vol. 2 (Academic press New York, 1978).

588

23. L. Tsang and J. A. Kong, *Scattering of electromagnetic waves: advanced topics*, vol. 26 (John Wiley & Sons, 2004).

589

590

24. M. B. Van Der Mark, M. P. van Albada, and A. Lagendijk, "Light scattering in strongly scattering media: multiple scattering and weak localization," Phys. Rev. B **37**, 3575 (1988).

591

592

25. L. L. Foldy, "The multiple scattering of waves. i. general theory of isotropic scattering by randomly distributed scatterers," Phys. review **67**, 107 (1945).

593

26. L. G. Henyey and J. L. Greenstein, "Diffuse radiation in the galaxy," The Astrophys. J. **93**, 70-83 (1941).

594

27. B. Bret, "Multiple light scattering in porous gallium phosphide," Ph.D. thesis, Universiteit Twente (2005).

595

28. S. Lönn and D. Kroon, "Utilizing enhanced backscattering for determination of scattering properties in turbid media," Lund Reports At. Phys. (2011).

596

597

29. J. Zhu, D. Pine, and D. Weitz, "Internal reflection of diffusive light in random media," Phys. Rev. A **44**, 3948 (1991).

598

599

30. M. Born and E. Wolf, *Principles of optics: electromagnetic theory of propagation, interference and diffraction of light* (Elsevier, 2013).

600

601

31. U. Schwertmann and R. M. Cornell, *Iron oxides in the laboratory: preparation and characterization* (John Wiley & Sons, 2008).

602

32. T. Sugimoto, A. Muramatsu, K. Sakata, and D. Shindo, "Characterization of hematite particles of different shapes," J. colloid interface science **158**, 420-428 (1993).

603

604

33. S. Potin, O. Brissaud, P. Beck, B. Schmitt, Y. Magnard, J.-J. Correia, P. Rabou, and L. Jocou, "Shadows: a spectro-gonio radiometer for bidirectional reflectance studies of dark meteorites and terrestrial analogs: design, calibrations, and performances on challenging surfaces," Appl. optics **57**, 8279-8296 (2018).

605

606

607

34. A. Krywonos, "Predicting surface scatter using a linear systems formulation of non-paraxial scalar diffraction," Ph.D. thesis, University of Central Florida (2006).

608

609

35. M. R. Querry, "Optical constants," Tech. rep., Missouri Univ-Kansas City (1985).

610

36. A. H. M. J. Triaud, "Refractive indices of hematite," http://eodg.atm.ox.ac.uk/ARIA/data?Minerals/Hematite/(Triaud_2005)/hematite_Triaud_2005.ri (2005).

611

612