1. INTRODUCTION

Multiple wave scattering by discrete random media has been studied for decades due to its relevance in a wide array of fundamental and technological applications with several textbooks covering its physics and principles [1–3]. The role of spatial correlations, in particular, is drawing increasing attention for the rich array of fundamental mesoscopic physics [4–10] and potential applications that it could impact, ranging from quantum optics to energy harvesting and radiative cooling [11–15]. Despite the constant advancements, however, the study of light transport through an ideally simple system of scattering spheres in a homogeneous medium still does represent an overwhelmingly complex vector problem, hindering our ability to answer many basic questions lying at the core of this field. Arguably the most relevant of these questions regards the turbidity of a scattering layer. Indeed, increasing particle density enhances opacity only up to a certain level, above which the overall scattering efficiency is suppressed by the close vicinity between neighboring particles. An optimal density naturally arises from this interplay, yet its value is only vaguely expected to lie somewhere between a volume fraction of 20% and 50% [16], and even such broad limits are often debated, with evidence supporting both lower [17–20] and higher [21,22] optimal densities values.

Other similarly inconsistent claims are frequently found in the literature, most of which are impossible to compare directly due to different experimental conditions. A typical example is related to the use of titanium dioxide (${\mathrm{TiO}}_2$) nanoparticles in coating and paint applications, which we will use as a reference test case in this work due to its generality, technologic relevance, and long-standing research efforts on the topic [23–27]. Indeed, while there is general agreement on the fact that scattering nanoparticles for an optimal white paint formulation should have a diameter of roughly 200 nm [18,28,29], little else is agreed upon when it comes to the optimal density, spatial correlations, and degree of clustering. The latter is another exemplary issue of debate: while flocculation and particle clustering is typically seen as a detrimental factor [29–33], numerical and experimental evidence exist showing that light scattering is either insensitive to clustering [34] or even enhanced by it [35]. Similarly, even the onset of the so-called “dependent” scattering regime is often questioned, being generally set around 1% volume fraction [18–20,36] with experimental evidence suggesting values well above that [28,37], possibly up to 20% [38].

Reasons for such diverging conclusions arise from the large number of experimental parameters involved in sample preparation, including the choice of nanoparticles size, inorganic coatings, polydispersity, resin and binder formulations, drawdown system, chemical dispersants, solvents, and extenders. Adding to this, conflicting interests between industrial players fostering an increased use of either scatterering particles or inter-particle spacers also seem to bias findings and claims in part of the literature [33]. Finally, even when spatial correlations are explicitly considered, e.g., exploiting either steric or electrostatic repulsion between particles, the actual effect on the three-dimensional (3D) point statistics is hard to evaluate or quantify, being at best only indirectly inferred by time-consuming analysis of two-dimensional (2D) micrographs [29]. As can be imagined, similar inconsistencies affect the optical characterization of all discrete random media beyond the illustrative case of ${\mathrm{TiO}}_2$ in paint coatings.

The aim of this work is to provide a comprehensive framework for rigorously surveying scattering properties in discrete random media with a focus on the role of spatial point statistics. To date, in fact, state-of-the-art methods consist of combining different approximations depending on the average particle density. Even so, the collective scattering approximation based on the structure factor is shown to fail above a particle density of $\sim 20\%$, with attempts aimed at including near-field contributions improving the agreement only qualitatively and in a limited density range [39,40]. Here, we propose a phenomenological approach that is not based on any simplifying hypothesis and, therefore, is not subject to validity ranges, providing a useful benchmark for further modeling. To this purpose, we adopt a numerical approach allowing us to define, quantify, and independently vary all physical parameters of the scattering medium one at a time, encompassing for the first time, to the best of our knowledge, the whole parameter space in terms of both volume concentrations and spatial correlations. Here, we focus on the specific type of correlations induced by the packing of hard spheres, but our approach can be applied to any kind of spatial correlation in discrete systems, including deterministic aperiodic [41], potential-based [42], hyperuniform [8], or random fractal [43] structures. Owing to this richer characterization, novel insight is gained by analyzing the interplay between particle density and spatial correlations in terms of two key functionals of the 3D point statistics, namely the first moments of the nearest-neighbor distance and pore-size distributions. Based on these figures, we clearly identify the configurations where near-field effects weaken the overall scattering strength [40,44,45]. In addition, we reveal the existence of a second set of spatial arrangements, where, unexpectedly, a weakening of the scattering strength is provided by perfect particle segregation, shedding new light on the reason for many inconsistent results in the literature.

2. NUMERICAL RESULTS

A. Design of Sample Configurations

We consider cylindrical discrete random media containing nonabsorbing spherical particles with a diameter of $d=200\mathrm{nm}$ embedded in a uniform host environment. The refractive index is set to 2.67 and 1.5 for the particles and host matrix, respectively, in accordance with typical values used in the literature to describe paint film components [29,32]. All cylindrical configurations have a diameter/height ratio strictly above 4 to avoid finite-size effects [40].

Particle configurations with a range of volume fractions have been obtained with a molecular dynamics code based on the Lubachevsky–Stillinger method [46], which allows us to generate disordered packings of spheres up to a volume fraction of $\sim 64\%$ (corresponding to the theoretical limit for an infinite compression rate [47]). A slightly bidisperse population of spheres was used with a radii ratio of 0.98 to reduce the formation of ordered subdomains in the highly packed configurations. Starting from a set of $P=22$ packing fractions in the range of $0.01,0.04,0.07,\dots ,0.64$, we generate $P(P+1)/2=253$ independent configurations corresponding to every combination of volume ($f_{\mathrm{v}}$) and packing ($f_{\mathrm{p}}\ge f_{\mathrm{v}}$) fractions compatible with nonoverlapping spheres. In other words, each configuration is characterized by a certain degree of correlations (inherited by an original packing of density $f_{\mathrm{p}}$), and a filling fraction obtained by properly rescaling the point pattern until the target $f_{\mathrm{v}}$ is reached. Note that this approach, which conveniently allows us to introduce spatial correlations in a controlled way and independently of the final volume fraction, is equivalent to the packing of core–shell particles with an inner diameter $d$ and an outer diameter $d_{\mathrm{excl}}=d(f_{\mathrm{p}}/f_{\mathrm{v}}{)}^{\mathrm{1/3}}\ge d$, such that the shell has the same refractive index as the surrounding medium acting as a spacer [11,48]. Following this analogy, densities are defined as $f_{\mathrm{v}}=\alpha d^3$ and $f_{\mathrm{p}}=\alpha d_{\mathrm{excl}}^3$ with $\alpha =\mathrm{\pi }\mathit{N}/6\mathit{V}$. A few exemplary configurations with different levels of particle density and spatial correlations are shown in Fig. 1.

 

Fig. 1. Rendering of illustrative configurations relative to values of $f_{\mathrm{v}}$ and $f_{\mathrm{p}}$ equal to 0.01, 0.04, 0.16, and 0.64. Panel (a) shows configurations with an equal particle/area ratio ($\mathit{N}/\mathit{A}=1.125·{10}^2{\mathrm{\mu m}}^{-2}$) and variable thickness ($0.75\mathrm{\mu m}\le L\le 48\mathrm{\mu m}$) comprising up to $3.6·{10}^6$ spheres. Panel (b) shows configurations with equal thickness ($L=3\mathrm{\mu m}$) and variable density ($7{\mathrm{\mu m}}^{-2}\le \mathit{N}/\mathit{A}\le 450{\mathrm{\mu m}}^{-2}$) comprising up to $8.0·{10}^4$ particles per sample. For simplicity, particles beyond a diameter/height ratio of 8 are discarded in the actual $T$-matrix calculations.

Download Full Size | PDF

To obtain a complete picture of how density and correlations affect the overall hiding power of a scattering slab, we consider two distinct sets of sample geometries, allowing us to address the problem from two complementary perspectives. In the first case [Fig. 1(a)], we generate different arrangements of particles while keeping a constant density per unit area. This approach is suitable to describe how the independent scattering approximation breaks down with increasing density, as was also investigated in a recent work with similar numerical simulations on few smaller aggregates [40]. The alleged rationale behind this sample design choice is that Ohm’s law for transmission through a scattering slab does not depend on the sample thickness at a fixed particle density per unit area. However, using Ohm’s law as a reference can lead to imprecise conclusions. First, one would need to consider extremely large samples in order to apply the diffusive approximation, well beyond those considered so far. Second, even in the simple diffusive approximation, additional terms must be considered that are not accounted for in Ohm’s law, resulting in a different scaling with thickness [49]. Finally, transmission in the diffusive regime is also substantially affected by the effective refractive index of the scattering slab, which clearly increases with increasing volume density [45,50,51]. Notably, disregarding the increased refractive index mismatch at the boundaries of a scattering slab is a known source of significant overestimation of the transport mean free path [52,53], highlighting the importance of a rigorous approach that does not rely on the hypotheses of approximated models.

Therefore, in addition to this first set of configurations, we also present a second complementary survey, where we consider slab samples with fixed thickness [Fig. 1(b)], looking for particle configurations providing the highest turbidity at a given wavelength. Besides being a more relevant quantity from an application point of view, absolute maximum integrated reflectance represents also a more rigorous parameter compared to the rather ill-defined transition between independent and dependent scattering regimes [54,55]. More importantly, for both sets of configurations, our approach based on the statistical properties of the particle point patterns allows us to provide an interpretation that does not rely on any approximated transport theory.

B. $T$-Matrix Calculations

Several numerical studies have been published in recent years, describing rigorous solutions of Maxwell’s equations in dense discrete random media. In this respect, the superposition $T$-matrix method provides a convenient approach, as it does not require any underlying mesh or regular point lattice to describe the investigated medium [3]. Even so, to date, most studies have considered only small clusters of particles or individual configurations [31,32,35,40,56–61] due to computational limitations.

In this study, we use a freely available, CUDA-accelerated implementation of the $T$-matrix method named CELES [62] to largely reduce the computational burden. By leveraging massively parallel execution on graphical processing hardware, a block-diagonal preconditioner and a look-up table approach for the evaluation of costly functions, CELES allows us to efficiently address large electrodynamics problems on inexpensive consumer hardware. All calculations are performed at a wavelength of ${\lambda }_0=532\mathrm{nm}$ using a perpendicular Gaussian beam illumination with a beam-waist size at least 8 times smaller than the cylinder diameter. For every configuration, we recorded the overall transmitted and reflected power integrated over the forward and backwards hemisphere, verifying energy conservation to a relative error of $\sim {10}^{-3}$.

Figures 2(a) and 2(b) show the obtained reflectance values for all different configurations, along with a cubic spline adaptation that allows us to smooth the fluctuations of the individual disorder realizations. Additional statistical averaging is inherently provided by the large impinging beam width and size of the samples. As expected, in the case of constant particle density, we find that overall opacity is almost constant at low volume fractions and drops abruptly above a 20% density. Similarly, the fixed-thickness configurations exhibit a maximum opacity around a volume fraction of 25%, which is in agreement with previous results. More interestingly, in both sets of calculations, we find a consistent trend where maximum opacity shifts from highly correlated configurations to maximally uncorrelated ones going from low to high-volume fractions. The net result is shown in the insets of Fig. 2, which considers for simplicity the overall reflectance along the two extreme values of the packing fraction, namely $f_{\mathrm{p}}=f_{\mathrm{v}}$ and $f_{\mathrm{p}}=0.64$. As can be seen, a crossover is observed in both cases, switching between two opposite regimes at a volume density of $\sim 30\%$.

 

Fig. 2. Numerical results and cubic spline adaptation of integrated reflectivity calculations. Panels show data relative to (a) fixed particle area density and (b) fixed-thickness configurations. Edges of the phase space corresponding to $f_{\mathrm{p}}=f_{\mathrm{v}}$ and $f_{\mathrm{p}}=0.64$, representing, respectively, the least and maximum amount of spatial correlations, are plotted in the inset. Pairs of configurations most impacted by spatial correlations are highlighted in panel (b).

Download Full Size | PDF

It is interesting to compare these results to the enhanced scattering strength that few authors claim to be achievable by acting solely on an optimized inter-particle spacing (see, e.g., discussion in Ref. [33] and references therein). In this respect, our controlled calculations cast a rigorous bound on the role of spatial correlations to improve reflectance at the wavelength of interest. While we indeed find that perfect particle segregation (corresponding to the so-called cermet topology) enhances opacity up to moderately high densities, the opposite is true after the crossover. A more quantitative take on this issue can be obtained by considering the fixed-thickness survey of Fig. 2(b). Limiting ourselves to the configurations with the highest ($f_{\mathrm{p}}=0.64$) and least ($f_{\mathrm{p}}=f_{\mathrm{v}}$) degrees of spatial correlations, we seek these pairs of samples exhibiting equal reflectance, while having the largest relative difference of particle volume concentrations. This is a relevant figure of merit especially for paint manufacturers, as it directly connects to the possibility of having two equally performing paint formulations using maximally differing amounts of ${\mathrm{TiO}}_2$ pigments—an expensive component in paint systems. As shown in the inset of Fig. 2(b), we find that an uncorrelated configuration with a particle volume concentration of $f_{\mathrm{v}}=0.0765$ achieves the same reflectance ($R=0.57$) of the maximally correlated configuration with a particle concentration of just 0.0635, corresponding to a 17% ${\mathrm{TiO}}_2$ reduction. Conversely, the opposite result is observed at higher densities, where a total reflectance of $R=0.68$ can be obtained, saving up to 9.3% of particles by completely avoiding spatial correlations, e.g., moving from a correlated configuration at $f_{\mathrm{v}}=0.436$ to a maximally random one at a smaller volume density of $f_{\mathrm{v}}=f_{\mathrm{p}}=0.396$.

Another prominent configuration that we successfully identify is that with the highest absolute reflectance, obtained at $f_{\mathrm{v}}=0.25$ and $f_{\mathrm{p}}=0.49$ for the fixed-thickness survey. The result of this particular calculation, corresponding to a total reflectance $R>77\%$ achieved in just 3 μm of thickness, is particularly remarkable if we consider the relatively low refractive index contrast of 1.78 between the ${\mathrm{TiO}}_2$ nanoparticles and the host matrix, and it represents a good estimate of the ultimate scattering strength that can be achieved in these systems at typical visible frequencies. Figure 3(a) shows a cross-cut along $f_{\mathrm{v}}=0.25$ of the reflectance data plotted in Fig. 2(b), featuring a peak for an intermediate degree of correlations. Based on the optical and geometric parameters of the simulated structures, we solve the radiative transfer equation for these configurations using the Monte Carlo method [63,64] in order to estimate an associated transport mean free path $l^{*}$. To this purpose, we ran several Monte Carlo simulations with ${10}^9$ photons each, considering an average slab thickness of $L+d/2$ with an uncertainty of $\pm d/2$. The effective refractive index of the scattering layer is estimated numerically from the exctinction of the coherent (ballistic) component of a plane-wave excitation for different configurations [65] (see Section 1 of Supplement 1); its exact value and dependence on the degree of correlations are nonetheless found to only marginally affect the inverse solution due to the static nature of the integrated reflectance. Even at fixed particle density, optimizing the degree of correlations allows to decrease the transport mean free path length by more than 10%, reaching an estimated minimum value of 0.43 μm. It is worth noting that this figure, corresponding to a value of $k_{\mathrm{eff}}{l}^{*}>9$, is still far from the expected transition to the strong localization regime. This is due to the moderately low refractive index contrast between the ${\mathrm{TiO}}_2$ nanoparticles and the host medium: much smaller values for $k_{\mathrm{eff}}{l}^{*}$ have been reported in the case of dry aggregates of nanoparticles, even though they are not sufficient to reach the Anderson localization regime [66,67]. Comparable figures for the total reflectance are obtained at different wavelengths in the visible range, as reported in Section 2 of Supplement 1. Figure 3(c) shows a cross-cut of the electric-field distribution in the $xz$ plane of the most opaque configuration illuminated with a Gaussian beam from the negative $z$ direction, illustrating how a large fraction of incoming radiation is back-reflected within a few units of $l^{*}$.

 

Fig. 3. Maximum integrated reflectance for equal-thickness samples. Panel (a) shows a cross-cut of Fig. 2(b), showing that the maximum opacity is reached for an intermediate degree of spatial correlations. Panel (b) shows the corresponding estimated transport mean free path $l^{*}$ obtained by a Monte Carlo fit of the total reflectance. Error bars correspond to a thickness uncertainty between $L\le L+d/2\le L+d$. The absolute value of the total electric field in the $xz$ plane for the configuration with $f_{\mathrm{v}}=0.25$ and $f_{\mathrm{p}}=0.49$ is depicted in panel (c).

Download Full Size | PDF

These observations underline the twofold role of spatial correlations in either enhancing or weakening the overall turbidity of a scattering layer based on the volume fraction occupied by the scattering particles. The onset of the high-density behavior challenges the common view that high-density configurations are associated with a so-called “particle crowding” effect that is detrimental to the overall opacity of the scattering layer. On the contrary, we find that perfect separation between particles does eventually weaken opacity above a certain filling fraction. To the same extent, we also show that the simple interpretation of a near-field assisted regime for light transport that has been recently proposed [39,40] provides an incomplete description for the observed behavior for intermediate and moderate particle densities, as it would predict a monotonic decrease of scattering strength with decreasing inter-particle distance. On the other hand, the typical “particle crowding” picture is recovered at lower densities and is particularly evident for the fixed-area/density configurations, where the presence of strong spatial correlations is capable of preserving high scattering performance up to higher particle concentrations. A consistent picture is also obtained on a broader wavelength range (see Section 2 of Supplement 1), where it is shown that high spatial correlations are indeed not generally associated with optimal scattering performances.

3. SPATIAL POINT ANALYSIS

In order to provide a more quantitative explanation of the observed behavior, we introduce two useful descriptors associated with the different roles of spatial correlations. To this purpose, we use Voronoi tessellation [68,69] as a versatile framework to define and test several functionals of the 3D point pattern statistics of the scattering configurations. A few notable examples have already been reported in the literature, also in connection with the local homogeneity of point configurations and their consequences for effective medium theories in electrodynamics [69,70]. In our case, we exploit Voronoi tessellation as the most computationally efficient tool to define the list of all neighbors to each particle, as well as to calculate the pore-size distribution of the medium with the trial sphere method, based on the Voronoi cells where the trial sphere centers fall [71]. We choose these two figures because of their intuitive connection to the two main mechanisms weakening the turbidity of a scattering layer, namely the near-field contacts between neighboring spheres and the high spatial homogeneity arising in highly correlated configurations.

If we consider parameters normalized to spheres of unit diameter, the mean distance between nearest-neighbor particles for a statistically isotropic 3D packing of hard spheres is defined as [72]

Examples of numerically evaluated probability density functions for ${\delta }_{\mathrm{nn}}-1$ and $2{\delta }_{\mathrm{p}}$ are shown in Figs. 4(a) and 4(b) for a representative set of particle densities. The final values of $⟨{\delta }_{\mathrm{nn}}⟩-1$ and $2⟨{\delta }_{\mathrm{p}}⟩$ have been evaluated numerically on ensembles of $8·{10}^6$ packed spheres and are shown in Fig. 5 for the whole space of possible configurations. In the case of maximally correlated configurations, the distribution of ${\delta }_{\mathrm{nn}}$ exhibits a narrow peak that spans from 1 to 4, i.e., from 200 to 800 nm, which might suggest the presence of sharp spectral features as seen, e.g., in 2D inverse structures [48]. Nonetheless, our numerical results do not exhibit any such feature, either as a function of correlations (Fig. 2) or as a function of wavelength (see Section 2 of Supplement 1). Several reasons explain this apparent difference with the 2D case: first, in 3D, TE and TM polarizations cannot be decoupled into a scalar description. Second, we study the case of direct structures where the high-permittivity inclusions do not form a continuum phase—a condition that is inherently associated with weaker dispersion [74], especially at the low index contrast of our structures. Finally, the use of bidisperse sphere packings further weakens any wavelength dependence. In this respect, our results are consistent with previous surveys that have shown a weak dependence of scattering properties with spatial correlations [42] and wavelength [61].

 

Fig. 4. Illustration of structural descriptors and their typical probability density functions as derived from auxiliary Voronoi tessellations for $f_{\mathrm{v}}=f_{\mathrm{p}}=\{0.01,0.04,0.16,0.64\}$ (solid lines) and for $f_{\mathrm{v}}=0.01$, $f_{\mathrm{p}}=\{0.04,0.16,0.64\}$ (dotted lines). Panel (a) refers to the normalized surface-to-surface distance between nearest-neighbors particles ${\delta }_{\mathrm{nn}}-1$. Panel (b) shows the distribution for the normalized pore diameter $2{\delta }_{\mathrm{p}}$ estimated using the trial sphere method with ${10}^6$ trials. Insets show simplified illustrations of the structural descriptors in a 2D geometry. All distributions have been evaluated on independent configurations containing $8·{10}^6$ particles.

Download Full Size | PDF

 

Fig. 5. Structural descriptors for the ensemble of particles over the space of configurations. Panel (a) shows the average surface-to-surface distance between nearest-neighbor spheres. Panel (b) shows the average size (diameter) of empty regions in the particle ensemble. All quantities are normalized by the particle diameter $d=200\mathrm{nm}$.

Download Full Size | PDF

Apart from the trivial cubic-root scaling along $f_{\mathrm{v}}$, it is interesting to comment on the opposite behaviors exhibited by the nearest-neighbor distance and the average pore size along the direction of increasing spatial correlations. As can be seen, the former always increases with increasing spatial correlations, resulting in a reduced near-field coupling between neighboring particles. Conversely, high spatial correlations—corresponding to perfect particle segregation—are systematically accompanied by a shrinkage of the average pore diameter. This effect is connected to a homogenization of the scattering medium as local density fluctuations are progressively suppressed, which eventually overcompensates for the increased nearest-neighbor spacing to give a decreased scattering strength. In this respect, the dependency of the overall turbidity on spatial correlations emerges from the interplay between these two opposite effects, which are weighted differently depending on the physical parameters of the problem.

To get a clearer overview of this interplay, we build a phase diagram representation to highlight the twofold role of spatial correlations in the space of all configurations tested, as shown in Fig. 6. For this purpose, we rely on the reflectance data calculated for our fixed-thickness survey [Fig. 2(b)], since it allows us to unambiguously identify the structural parameters resulting in the highest opacity. In particular, we define two regions characterized by an opposite behavior based on the sign of $\partial R/\partial f_{\mathrm{p}}$. In other words, we propose a classification of the particle configurations depending on whether increasing the degree of spatial correlations enhances or weakens the overall scattering strength of the random medium. The boundary between these two phases is formed by the loci of the scattering configurations that provide the maximal reflectance for a fixed filling fraction, shifting from highly correlated to maximally uncorrelated arrangements with increasing particle density.

 

Fig. 6. Phase diagram for the role of spatial correlations in terms of the pore size and the nearest-neighbor distance based on the numerical results of Fig. 2(b). Structural parameters are shown normalized either by the diameter $d$ of the scattering spheres (lower and left axes) or by the wavelength in the host medium $\lambda ={\lambda }_0/1.5$ (upper and right axes). Circles indicate the location of the simulated samples. Dotted lines highlight the combination of structural parameters resulting in the highest turbidity.

Download Full Size | PDF

Having units of length, these structural parameters can be directly compared to the wavelength of incident light to get a more direct connection to the physics of the problem. For this reason, we cast our phase diagram in units that are normalized both by the particle diameter $d$ and the wavelength in the host medium of the scattering film $\lambda ={\lambda }_0/1.5$. In this framework, the configuration exhibiting the highest reflectance corresponds to a separation between nearest-neighbor particles of $0.261d\sim \lambda /6.8$ and a pore diameter size of $0.375d\sim \lambda /4.7$. Notably, due to the type of packing and rescaling algorithm that followed to generate these configurations, these optimal parameters can be reproduced, in principle, by packing together core–shell particles with an internal radius of 100 nm for the ${\mathrm{TiO}}_2$ core and a $\sim 125\mathrm{nm}$ outer shell with an index of refraction matched to the host matrix.

4. CONCLUSIONS

In this work, we tried to answer two simple questions regarding multiple scattering in discrete random media, namely: “What is the optimal density for a system of scattering particles?” and “Given the optimal particle volume concentration, do spatial correlations enhance or weaken the overall turbidity?” To investigate these issues, we addressed the problem from a twofold perspective, each pertinent to different applications. In one case, we considered a finite amount of “resources” (e.g., number of ${\mathrm{TiO}}_2$ nanoparticles) per unit area without any constraint on how to arrange them (e.g., their dilution ratio in the resin matrix). Indeed, understanding the role of spacing and dilution of nanoparticles is of great interest to paint manufacturers, as ${\mathrm{TiO}}_2$ is an expensive component of the formulation of paint systems and optimizing its scattering efficiency is a key cost/performance factor. On the other hand, we looked at the same problem using an unconstrained number of particles and a fixed-thickness gap for where to fit them. Arguably, the underlying question in this case—how far can we push density to gain opacity—is even more interesting, as it inherently accounts for how the particle density inevitably affects spatial correlations and the effective permittivity. This approach is in turn relevant for the design of diffuse reflectors for photon-management applications in the illumination and photovoltaic fields that are likewise subject to similar thickness constraints.

As we have shown, answering these simple questions revealed a surprisingly complex role played by spatial correlations and inter-particle separation, whose effect is not simply that of universally enhancing the overall opacity. Indeed, we found that particle configurations can be classified into two groups, depending on whether inter-particle separation enhances or diminishes their overall scattering strength. This allowed us to provide well-defined structural specifications relative to the highest scattering performance at a given wavelength, as well as to those density parameters that are most impacted by the introduction—or the complete removal—of spatial correlations.

By performing a complete survey in the phase space of possible configurations, we were able to set a limit for the shortest possible transport mean free path that can be achieved in a discrete random medium of spherical particles, considering the optical parameters of a typical white paint formulation (size parameter $\mathrm{\pi }d/\lambda =1.77$ and index contrast $n=1.78$), a figure with potential implications both for applied and fundamental aspects of the physics of light transport [66,75]. In this case, we observed that an intermediate amount of correlations is able to induce a reduction of $l^{*}$ of $\sim 13\%$ compared to its value for a random packing of particles.

A lot of work remains to be done beyond the oversimplified model of equal-sized spheres. Relevant aspects remain to be considered, such as arbitrary polidispersity, more realistic clustering distributions, and particle interaction potentials, or the fact that high-volume concentrations are usually associated with air inclusions as the host matrix cannot perfectly wet all particles. Similarly, it could be interesting to further study multiple scattering for optically soft spheres, where evidence exists for yet a different interplay with correlations [19], as well as for core–shell or spheroidal particles, which are also sometimes associated with enhanced scattering strength [76,77]. Beyond the illustrative case of common paint coatings, further exciting opportunities are open to studying more exotic media that, due to their long-range correlations or multi-scale nature, have been prohibitive for studying numerically in representative 3D geometries.