Topological invariance in whiteness optimisation

Maximizing the scattering of visible light within disordered nano-structured materials is essential for commercial applications such as brighteners, while also testing our fundamental understanding of light-matter interactions. The progress in the research field has been hindered by the lack of understanding how different structural features contribute to the scattering properties. Here we undertake a systematic investigation of light scattering in correlated disordered structures. We demonstrate that the scattering efficiency of disordered systems is mainly determined by topologically invariant features, such as the filling fraction and correlation length, and residual variations are largely accounted by the surface-averaged mean curvature of the systems. Optimal scattering efficiency can thus be obtained from a broad range of disordered structures, especially when structural anisotropy is included as a parameter. These results suggest that any disordered system can be optimised for whiteness and give comparable performance, which has far-reaching consequences for the industrial use of low-index materials for optical scattering.


Introduction
Light diffusion in random structures is a fairly well established phenomenon that has a profound impact in fundamental research as well for applications: from random lasing to light harvesting to imaging and white paints. 1,2 n fact, bright whiteness is typical of materials where light undergoes multiple scattering events over a broad range of wavelengths.Such diffusive reflectance usually requires high refractive index media (with related safety concerns 3 ) or thick scattering layers to ensure sufficient amount of scattering events.Recently, natural examples of disordered materials have drawn much attention due to their ability to express record transport mean free path using low refractive index components.In particular, the white beetles, e.g.Lepioda stigma and Cyphochilus sp.-which exploit an anisotropic chitin network (n c ≈ 1.55) to achieve a bright whiteness -have become the poster children of disordered photonics [4][5][6][7][8][9] .Several investigations have attributed these optical properties to the specific characteristics of the random network structure of the materials.However, due to limited "metrics" used to compare and analyse random structures, the role of the various features of these complex nanostructures in the scattering properties remains still unclear.
Disordered systems are loosely defined as inhomogeneous materials that exhibit small (but not long range) structural correlations at relevant length distances r.The most popular measures used to characterise disordered systems in optics are based on first and second order statistics such as form-and structure factors, 2-point correlation function S 2 (r), filling fraction, V 0 (:= S 2 (0)), and correlation length l c .It has been suggested that unique structural similarities between disordered systems can be inferred from identical reflectance spectra R(λ ) and S 2 (r) 10 .Yet in the case of random network structures inspired by white beetles, many different synthetic models have been reported to match the beetle structures using these functionals [10][11][12][13] , suggesting its non-uniqueness.This is not surprising, given that in the field of disordered studies, which predates the study of disordered photonics, it is well known that only in few exceptions, e.g. in the case Gaussian random fields, also known as Gaussian Processes (GP), the structural properties of stochastic fields are uniquely determined by second order statics alone 14 , and for an arbitrary ones more robust topological measures are needed for unique characterisation.Therefore we propose to use a particularly useful metrics for this quantification: the Minkowski functionals V 0 ,V 1 ,V 2 ,V 3 , which are based on surface integrals, and are proportional to the filling fraction, surface area, integrated mean and total curvature respectively 15,16 .
Using such topological identifiers we tested the scattering efficiency of a series of disordered structures synthesised in silico using both top-down and bottom-up models to cover a wide range of topological features in terms of porosity, connectedness, branching and length scale.We quantified these differences using the Minkowski functionals and measured the reflection properties using finite-difference time domain (FDTD) simulations (cf.Fig. 1).Then, by correlating the two, we were able to understand how each feature contributes to overall scattering properties, demonstrating that disordered systems are sufficiently explained by second order statistics alone.Finally using the the tensorial Minkowski measures 17,18 , we also quantified structural anisotropy and investigated its role in whiteness optimisation.

FDTD Simulations
InSilico Synthesis The optical results are investigated in respect to structural descriptors which include second order statistical measures such as 2-point correlation functions S(r), filling fraction V 0 , and correlation length l c , but also more unique structural measures like Minkowski functionals, V 0 ,V 1 ,V 2 ,V 3 which are the integrated volume, surface area, mean and Gaussian curvatures respectively (colours represent piecewise contributions), and also Minkowski tensors W r,s ν , and the associated eigenvalue ratios β r,s ν which measure the structural anisotropy.

In silico synthesis of disordered structures
As the size of the parameter space for disordered structures is enormous, we selected 10 different model systems (GP1-2, FC1-5, and SD1-3, cf. Figure 2a) using three different stochastic approaches and mapped the reflectivity landscape using line searches by varying i) the correlation length l c , ii) filling fraction V 0 and anisotropy of these systems.We believe that our selection covers a significant wide range of different topological features relevant to the topic of scattering optimisations.
In the first approach we used Gaussian Processes (GPs), which are completely defined by their mean m(r) and S 2 (r) (the latter is also known as covariance kernel in machine learning literature), to test the dependency of photonic properties of disordered structures on second order statistic alone.We selected two different correlation functions, one with sinc-type and other with squared exponential where l is a scaling parameter, to investigate the effect of a fixed correlation length, (GP1), vs. one with distribution of correlation lengths (GP2) (cf. Figure 2a), on the optical response.
In the second case, we used a popular phase-field approach, the Cahn-Hilliard (CH) model 19 to generate more complex structures.The CH model successfully describes the time evolution of demixing process (known as spinodal decomposition) of homogeneous mixture into separate domains and was chosen as it has recently been suggested as the formation mechanism behind the beetle scale structures 10 .We used three spinodal decomposition (SD) models, which we initialised using 50 %, (SD1), 30 %, (SD2), and 70 %, (SD3) filling fractions.
Although the CH model has been successfully used to simulate disordered structures in many different fields, it's mainly limited to bicontinuous networks and isolated micelle morphologies, that exclude the possibility to investigate e.g.tubular and cellular morphologies.A particularly interesting approach, that was developed to gain control over these features, is the Functionalised Cahn-Hilliard (FCH) model [20][21][22] which takes into account hydrophobic and mixing entropy effects, allow to control the surface curvature.Therefore in the third approach we adopted a known FCH protocol with 20 % initial filling fraction 20 , to simulate a variety of structures from colloidal, (FC1), and tubular, (FC2), to branched and cellular (FC3-FC5) ones).These structures are also interesting as they offer flexibility in visual and quantitative (using the V n ) matching of synthetic disordered systems with biological ones.

FDTD simulations
FDTD simulations were carried out using a broadband plane wave source and a reflectance monitor above the sample.The reflectance spectra were then integrated over the monitor area, and spectrally averaged to obtain total reflectance R tot .Because the parameter space, spanned by correlation length l c and V 0 , is rather large, and the FDTD simulations are computationally relatively costly, we carried out series of line scans by keeping the l c fixed and varying V 0 , and vice versa, for each model.Starting from l c = 300 nm and V 0 = 30% one would have expected (under the assumption that optical properties are not solely determined by 2nd order statistics) to see clear difference in reflectivity between different structures.Surprisingly not only were the average reflectances between different structures, very similar, from 0.49 to 0.53, except for the cellular structure FC5 with R ≈ 0.45, but they all behaved in a similar unimodal manner as function of V 0 = [10, 20, 30, 40, 50, 60]%, as shown in Figure 3a, with highest reflectivity reached between V 0 =30-50 % in agreement with earlier studies 23 .Next, we therefore continue to scan over correlation lengths between l c = [100 nm, 200 nm, . . ., 900 nm] while keeping the V 0 = 30 % fixed, cf. Figure 3b.Here it could be seen that the value l c = 300 nm was a convergence point for most structures.When moving to higher values, the differences became more pronounced and in favour of colloidal systems, like FC1 and SD2, and interestingly at the expense of inverse/cellular structures, such as FC5 and SD3.We then line scanned V 0 = [10, 20, 30, 40, 50, 60]% only to discover similar unimodal behaviour and an optimal region of V 0 =30-40 %, with the exception of more pronounced separation in averaged reflection between the different structures.

Role of anisotropy
As structural anisotropy is known to be an important feature affecting the reflectance properties of disordered systems 24 , we decided to use the Minkowski tensors, which are based on strong mathematical foundation of integral and convex geometry 17 , to quantitatively correlate the structural anisotropy to the optical response.In fact, to the best of our knowledge, a robust way to quantify the anisotropy of an arbitrary structure 24 is still lacking in the photonics community.Popular methods include the simple comparison of ratios of correlation lengths along different directions, but question how such values should be interpreted remains open, and can be rather inaccurate for latent anisotropy, as we will later demonstrate.
In 3D there are several different linearly-independent tensors, and a particularly suitable one for two-phase structures is which measures the distribution of surface normals n 18,25 .The degree of anisotropy can then be expressed as the ratio of minimal and maximal eigenvalue of the tensor  .which we will now on referred as β .An isotropic structure will have β ≈ 1 and lower values signify increased anisotropy.By replacing the scale parameters and diffusion coefficients in equations ( 1), (6), and ( 9) with ones that are independent along different axes we can extending our simulations to model also anisotropic structures 26 and using the Minkowski tensors, we can both induce and quantify structural anisotropy, and therefore systematically map its effect on the photonic response.Thus, for all the model types, we created a series of structures with anisotropy ranging between β ∈ [0, 1], while keeping the other parameters, V 0 = 0.3 and l c = 300 nm fixed.The values of β were roughly evenly spaced since a structure with an arbitrary anisotropy value has to be iteratively searched.
The result are shown Figure 4.In all cases the reflectance shows improvement with increasing anisotropy usually peaking to R ≈ 0.6 between β ∈ [0.4,0.6].It is interesting to note that, thanks to release of the X-ray tomography data sets of the beetle scales to public domain by Burg and coworkers 27 , we were also able to precisely quantify that that disordered structures inside both Lepioda stigma and Cyphochilus sp.scales are also anisotropic, β = 0.7 and β = 0.6 respectively.In the latter case there has been some speculation whether the structure is isotropic or not, but the quantitative analysis confirms that the original anisotropy claim 5,28 is valid.What is most surprising, for the simulated structures, is that the system that shows worst reflectivity in the isotropic case, FC5, reaches the highest reflectivity R ≈ 0.63 at β ≈ 0.15 among all the models, suggesting One should note that all the structures, at high anisotropy, start to resemble 1D multilayer-like systems, and therefore we would expect them also show high reflectivity due to consistency with literature reports 29 .In the supporting information we demonstrate that correlation lengths are merely shifted and by reoptimisation the higher reflectance levels can be recovered (cf. Figure S2).
Moreover, by comparing mean free paths as function of anisotropy, as shown in Figures 5 (cf.ESI for transport mean free path l t ) we observed the spectrally averaged mean free paths are rather similar between the different structures, and differences are mostly related to spectral variances.Furthermore these variances are decreased with increasing anisotropy with bottle neck like optima between β ∼ 0.2-0.7 suggesting anisotropy plays an important role in whiteness optimisation.

Feature analysis
So far we have mostly related topologically invariant features such filling fraction, correlation length, and anisotropy to reflectance properties.Ultimately one aims to arrive at an analytical model for predicting photonic properties of an arbitrary disorder structure via examining above mentioned structural properties.Deriving such model is very difficult task, and instead a 5/12 black box approach is needed.Machine Learning (ML) methods have become very popular in such problems.In particular deep learning based (and similar) ML methods are becoming popular for disorder structures investigations in photonics and material science 30,31 , and allow the use of raw 3D structures in the learning process.While those approaches are in principle very powerful in predicting the output (photonic) from the inputs (3D structure), given large enough training dataset, the challenge of interpreting the structure-property relationship between the two remains.Given that we have the various structural descriptors, such as Minkowski functions and anisotropy values, using them instead of the raw 3D structures for ML regression analysis, we can quantitatively estimate the importance of each features, and since each of those features physical interpretation also the inference will be more interpretable.
Thus we simply collect all the relevant structural features (a.k.a labels), for each simulated structure and reflectance spectra R(λ ), and use regression analysis with the following mapping where f min is the correlation strength (cf.ESI), β r,s ν are the anisotropy factors 18 , and the Betti numbers b 0 , b 1 , b 2 measure the number of particles, loops, and cavities respectively 32 .A practical way to assess the importance of each feature is to divide the data to training and test sets, and compare how accurately a ML model, trained with former dataset, can predict a particular feature from the latter.Features that are easy to predict correctly can be interpreted to have a higher significance for the output, (R(λ )), and thus imply a stronger structure-property relationship.
For instance using the popular Random Forest (RF) ML regression model for equation (4), we observed in Figure 6a), that most easily predictable features, from reflectance spectra, are the surface area V 1 , integrated mean curvature V 2 , and correlation length l c .Furthermore, using quantitative importance analysis, we can give percentual estimates of importance of each 15 features to reflectivity, with the five most important being V 1 (37 %), V 2 (13 %), l c (13 %), V 0 (11 %), β (8 %).These values should however be taken as tentative, since the importance analysis is rather sensitive parameter limits.E.g. structures with very low filling fractions, V 0 < 10 %, will have very low reflectance, and including them feature analysis will result in higher weighting of V 0 in importance, and the same is true for l c < 100nm.Therefore we performed the feature analysis above these limits, where we still have reasonably amount of reflectance.The suggestion that the surface area V 1 and the mean curvature V 2 are the most important features is a consequence of the fact that the colloidal systems (FC1 & SD2) seem to be most reflective.

6/12
This indicated that disordered assembly of particles would give the highest reflectance.To test this hypothesis we performed series of molecular dynamics particle simulations with the same correlation lengths l c =100-900 nm.To also check if the reflectance is dependent on the particle shape, we carried out the simulations with both spheres and tetrahedrons, as the latter has a higher integrated mean curvature, due to sharp edges, and increased surface area compared to the former.
The results (cf.ESI, Figure S1) indeed show that the spheres outperform the other systems in reflectivity, whereas tetrahedrons show only moderate reflectance levels, suggesting that the minimal surface area with high mean curvature, under the constrain of filling fraction and correlation length, is a relevant factor for reflectance.Finally, whether spheres are, in general, the optimal shape for particle system in respect to reflectance is interesting but out of the scope of this work.

Conclusion
Using different stochastic models we synthesised, in silico, a diverse set of disordered structures varying in shape, connectiveness, curvature, surface properties, volume, percolation, number of particles and anisotropy.We then measured their reflectance properties, quantified their structural differences, and analysed the structure-property relationships.Our extensively research revealed that all the investigated disordered systems exhibit similar unimodal behaviour in respect to filling fraction and correlations lengths.By utilising the Minkowski measures we were able for the first time quantitatively relate structural features, including anisotropy, to optical properties suggesting that with proper tuning of the second order statistical features (l c and V 0 ), and anisotropy brilliant whiteness can be achieved in any system.While colloidal systems have the highest reflectance due their to small surface area V 1 and high mean curvature V 2 , due to the lack physical support at optimal filling fraction V 0 =30-40 %, they are not a feasible solution in the context of material fabrication.Therefore whiteness must be realised with continuous structures.However in this case, we observed that choice of particular morphology becomes irrelevant as we have demonstrated.From an industrial and design point of view this is very fortunate, as one is not limited to trying to realise a particular morphology for optimal whiteness, but instead the efforts can be focused on optimising the before mentioned invariant parameters.
The results also suggest that one has to be careful when inferring the about underlying formation mechanism of biological disordered systems using optical response and S 2 (r) alone, due to non-uniqueness of these characteristics.Conversely this would also explain why there are many different disorder system in nature that exhibit whiteness [33][34][35] , as the results suggest that any disordered system generated by a tuneable mechanism can be optimised for brilliant whiteness.
In conclusion our work suggests that there is no unique route to brilliant whiteness, and instead it can be realised from many different starting positions.

Methods
All the in silico syntheses were carried out in a N × N × N cubic grid with periodic boundary conditions, and set to correspond to a L 3 =5 µm×5 µm×5 µm box.A value of N = 100 was used for all structures with l c ≥ 200 nm, and N = 200 otherwise.

Phase-field simulations
The Cahn-Hilliard model used for generating the spinodal decomposition structures SD1-3 is given by where M and ε are the mobility and the diffusion constants, and W ( f ) is the mixing energy between the two equilibrium phases, and was numerically solved using semi-explicit spectral method 40, 41 where ),with time step ∆t = 1 × 10 −3 , A = 1, M = 1, and the coefficient ε was chosen from [0.04,1.00] to reach the desired length scale.The FC1-5 models were created with Functionalised Cahn-Hilliard model 20,21 ∂ using the CUDA code of Jones 22 with η h = 5 and η m = [−7.25,−0.5, 3, 6, 10] for FC1-5 respectively.

Synthesis of anisotropic structures
Anisotropic structures where synthesis using a the method of Essery 26 where the mobility coefficient ε in equations ( 6), and ( 9) was replaced by tensor and a similar method for the GP models by anisotropic scaling (5).

Conversion to binary fields
The stochastic fields f (x) : R 3 → R synthesised with the different methods were converted to the two phase disordered structures I(x) : R 3 → {0, 1} using a simple thresholding scheme where I(x) is an indicator function , and 0 and 1 represent the empty and the solid phase, respectively, and p 0 is the threshold value.Thus the final filling fraction V 0 = I(x) of structures is determined by the choice of p 0 .While such systems might be physically difficult to realise, we are primarily interested in understanding the optical properties of disorder systems, and question about chemical synthesis of potential structures are beyond the scope of this paper.
For the FDTD simulation, the structures f (x) were imported to USCF chimera 42 and converted to STL stereolitography files using a similar level set scheme.

Particle simulations
The additional particle simulations (see ESI) were conducted with HOOMD-blue 43,44 package using molecular dynamics simulation with Langevin integrator and Weeks-Chandler-Andersen potential for spheres, and hard particle Monte Carlo simulation for the tetrahedrons.
where T b and Ttot are the ballistic and total transmission respectively, and the coefficient R = 0.23 was calculated using Maxwell-Garnett approximation.The ballistic transmission was measured using mode expansion monitor in Lumerical.
The mean free paths where calculated using the least squares slope solution [2]   l where are T tot , y = L for ls and l t respectively.

Figure 1 .
Figure 1.Workflow model of the disorder optics investigations.a A disordered structure, f (x), is synthesised using stochastic model in silico, and b imported to finite difference time domain (FDTD) optics simulator with periodic boundary conditions in x-and y-directions and a plane wave propagating along z to obtain broad band reflectance spectra, R(λ ).cThe optical results are investigated in respect to structural descriptors which include second order statistical measures such as 2-point correlation functions S(r), filling fraction V 0 , and correlation length l c , but also more unique structural measures like Minkowski functionals, V 0 ,V 1 ,V 2 ,V 3 which are the integrated volume, surface area, mean and Gaussian curvatures respectively (colours represent piecewise contributions), and also Minkowski tensors W r,s ν , and the associated eigenvalue ratios β r,s ν which measure the structural anisotropy.

Figure 2 .
Figure 2. a Investigated simulated structures morphologies including Functionalised Cahn-Hilliard (FC*), Gaussian Processes (GP*), and spinodal decomposition (SD*) using the Cahn-Hilliard model.b The second order correlation functions, S 2 (r), of these structures demonstrating that very different topologies are often indistinguishable using second order statistics alone.
Correlation length l c (nm)

Figure 3 .
Figure 3. Line searches along different filling fractions V 0 and correlation lengths l c for different morphologies.a,c Varying filling fraction with fixed l c = 300 nm and l c = 600 nm.b Varying correlation length with fixed V 0 = 30 %. d) Extrapolated surface for the most reflective system FC1..

FC5Figure 4 .Figure 5 .
Figure 4. Effect of anisotropy β on average reflectance.A value of β = 1 indicates complete structural isotropy and β = 0 full structural anisotropy.The lines present Gaussian Process fits to measured (FDTD simulated) reflectances.The data points have been omitted for clarity, and can be found ESI, Figure S3.Snapshots of FC1 structures are shown in the bottom to illustrate the effect of the anisotropy on a morphology.

Figure 6 .
Figure 6.Feature importance analysis using Random Forest regression.a) Plots of predicted and actual values of the various features (Only the 9 most relevant are shown for space considerations).A good correlation between predicted and actual value (diagonal distribution of points) is interpreted to signify high relevance to reflectance.b) Quantitative analysis of the feature importance.Value of 1 signify high and 0 low importance.

Figure S5 :Figure S6 :
FigureS5: The transport mean free path, l t , as function of wavelength λ, and anisotropy β. a) 3D plots of l t vs β, λ and b) marginal plots where the shaded area presents variation between spectral extremes (red curve for λ = 800 nm, and blue for λ = 300 nm) and the markers indicate spectral averages.The plots demonstrate how, especially on large wavelengths, the scattering can be increased with anisotropy optimisation