Emergent scattering regimes in disordered metasurfaces near critical packing

Domain: physics
Published: 2025-12-15
Editorial outcome: minor_revision

Abstract
Disordered metasurfaces provide a versatile platform for harnessing near- and far-field scattered light. Most research has focused on either particulate topologies composed of individual, well-identified metaatoms or, to a lesser extent, semi-continuous aggregate topologies without well identified inclusions. Here, we uncover an intermediate ‘critical packing’ regime characterized by metasurface morphologies in which a significant fraction of metaatoms begin to connect. We experimentally demonstrate that, at this threshold, the properties of the scattered light abruptly change and, via a statistical quasi-normal mode analysis, interpret this change as a marked transition in the statistics of the photon density of states. Unlike percolation in semi-continuous metal films, this transition affects not only the specular but also the diffuse components of the scattered light in a profound way. Our results introduce critical packing topologies as a novel design strategy for manipulating the spectral and angular characteristics of light using ultrathin optical coatings. Emergent functionalities include colour shifts in diffuse light driven by multiple scattering and surface whitening, with potential applications in display technologies—for example, to reduce glare in electronic screens.

Introduction
In solid-state systems, rapid changes in the electronic density of states (DoS) are often associated to spectacular phenomena, such as phase transitions1–3. In random photonic materials, the analogue of the electronic DoS is played by the photon density of states (PDoS) and the states are quasinormal modes (QNMs) with complex energies, owing to non-Hermiticity induced by leakage or absorption4. The PDoS contains crucial information on the optical transport properties of the material5–11, with important implications on light localization12–14, bandgap formation15,16 or lasing in the presence of gain with long-lived quasi-extended modes17. In this work, we report a phase-transition-like phenomenon in the context of dielectric disordered metasurfaces, revealing a close photonic analogue to some behaviour observed in electronic systems.

Optical metasurfaces are artificial surfaces in which the constituent inclusions (or metaatoms) are meticulously engineered from dielectric or metallic subwavelength elements to achieve tailored optical properties and functionalities4,18,19. Among these, disordered metasurfaces hold a prominent place, both from fundamental20,21 and applied22,23 perspectives. They have led to new designs for optical encryption24,25, light extraction26,27, transparent displays28, low-emissivity coatings29, non-iridescent colouration30–33, visual appearance design34, and optical diffuser35, among others.

This work introduces disordered metasurfaces operating at the transition between particulate and aggregate topologies, where metaatoms begin to connect (Fig. 1). In this unexplored regime, referred to as critical packing hereafter, we theoretically predict a sudden transition in the statistical distribution of QNMs when morphological parameters like the metaatom density or spatial correlation are varied. We experimentally verify that this transition results in significant changes in the behaviours of both specular and diffuse light, offering new opportunities for creating ultrathin fade-resistant coatings designed to implement advanced functionalities, unattainable through multilayer thin films22,23,30–34,36,37.

Results
Manifold visual appearances around critical packing

To experimentally investigate the light scattering transition between particulate and aggregate topologies, we designed a series of disordered metasurfaces written in a negative resist using electron beam lithography and then transferred into a 145 nm-thick silicon layer on a glass substrate (Fig. 2a and “Methods” section). The metasurfaces are circular with a 300 μm diameter. They are fabricated on the same sample to reduce fabrication variability. All the electron beam patterns consist of random arrays of square features, whose centres are positioned using a Poisson disk sampling algorithm40. The (x, y) coordinates of the centres are formatted into the graphic data system (GDS) files.

We vary three control knobs of the square features, their density, size and disorder correlation. Structural correlation is controlled by setting a minimum interparticle distance, a, between particle centres, resulting in an effective surface coverage or ‘packing fraction’ defined as p = \rho \pi a^2 / 4.

This process yields a comprehensive set of 27 metasurfaces. The designed patterns feature squares with lateral sizes ‘ = 95, 130 and 170 nm, densities \rho = 2, 5 and 10 μm^-2, and packing fractions p = 0.1, 0.3, and 0.5. These values represent nominal target parameters, but different topological states emerge during fabrication due to spatial overlaps between the square features at certain densities, sizes, or disorder correlations. Scanning electron microscope (SEM) images of all metasurfaces are shown in Fig. 2c, with a few previously displayed in Fig. 1. The metasurfaces can be categorised into three topological states: particulate, aggregate, and critical-packing topologies, where particles begin to touch.

Both the specular and diffuse reflection responses of the metasurfaces are measured using a gonio-spectrometer setup and a slightly focused supercontinuum laser. After normalising the measurements with the diffuse response of a reference diffuser, which offers high diffuse and Lambertian reflectance, we derive the bidirectional reflectance distribution function (BRDF)41. This multidimensional function characterises how the metasurfaces scatter light across all viewing directions, depending on wavelength, incident angle and polarisation. In addition to these quantitative measurements, the light scattered off the metasurfaces is visualised on a centimetre-scale hemispheric diffusive screen.

A halved ping-pong ball serves as the screen, with an elongated band of apertures along the plane of incidence to prevent blocking the incident beams and their specular reflections (Fig. 2b). This simple setup provides a direct visualisation of the colour and scattered intensity produced by the metasurfaces, offering a qualitative representation of their visual appearances in reflection. Photographs of the ping-pong ball, capturing the diffuse light from all 27 metasurfaces, are taken using a Canon EOS 1000D DSLR. Specular reflection is also photographed by aligning the camera with the specular direction. Different camera settings are used for specular and diffuse light images to prevent saturation, but the settings remain consistent across all recordings, allowing for direct intercomparison of the images.

The 27 × 2 photographs capturing both the diffuse components of scattered light at 45° incident angle and the specular components at a 10° incident angle under unpolarized light illumination are displayed in Fig. 2c. They are presented alongside corresponding SEM images of the metasurfaces, with the specular light photographs appearing as rectangular colour boxes positioned just to the left of the SEM images.

Particulate metasurfaces are primarily located in the leftmost column and the upper portions of the central and rightmost columns of Fig. 2c. With the exception of the singular metasurface in the upper right corner, which will be discussed separately, they exhibit two dominant, vivid hues, primarily determined by the resonance frequency of individual metaatoms. For uncorrelated disorder, p = 0.1, the colour and brightness of the diffuse light remain independent of the observation direction. In contrast, for correlated disorder, p = 0.5, brightness is reduced around the specular direction.

This trend is a direct manifestation of the engineered correlations, quantitatively captured by the static structure factor, S(r) (see Supplementary Note 1, Supplementary Fig. 1d and ref. 23). For weakly correlated patterns (p = 0.1), S(r) is nearly constant, resulting in the observed uniform, isotropic scattering. As correlation increases (e.g. p = 0.5), a pronounced dip develops in S(r) at small in-plane wave vectors (q || ~ 0), which directly causes the reduction of diffuse light around the specular direction. Concurrently, a strong peak emerges for q || ~ 2\pi \sqrt{\rho}, creating the bright, ring-shaped halo of enhanced scattering. The dispersive nature of this peak, where its angular position depends on wavelength, is what gives the halo its distinct colour. Furthermore, as the particle density \rho increases from 2 to 5 μm^-2 (moving from the left to the middle column in Fig. 2c), the primary peak of the structure factor shifts to shorter wavelengths (Supplementary Fig. 1b). This shift directly explains why the coloured halo appears to expand and change hue, providing another knob for tuning the visual appearance34. These observations are further analysed in Supplementary Note 1, which presents additional variations in vivid hues obtained by adjusting the metaatom size parameter.

Metasurfaces operating near critical packing appear in the middle and lower sections of the second and third columns. They display pastel hues with low saturation, incorporating a significant amount of white light. Some of these metasurfaces, particularly at packing fractions of 0.5 and 0.3, exhibit an angle-dependent brightness, with a noticeable reduction of diffuse light around the specular direction. This confirms that correlation in the initial electron-beam pattern continues to affect scattering properties, even as the metaatom overlap is significant.

Aggregate metasurfaces, which form at high density and low packing fraction, are primarily found in the lower part of the third column. They are characterised by high brightness, close to white, with soft hues that slightly vary between the specular and backward directions. Unlike particulate metasurfaces, this angular-dependent colouring does not stem from correlated disorder, as aggregate metasurfaces emerge from a weakly correlated configuration and do not exhibit a reduction in diffuse light intensity around the specular direction. A closer inspection reveals pronounced shifts in hue as the system approaches and crosses critical packing, including both a progressive blue-shift and a general whitening of the diffuse light. These spectral trends are not immediately obvious in Fig. 2 alone, but will be quantitatively analysed in detail in the following sections.

Specular colours tend to be relatively dull for aggregate metasurfaces, as well as for particulate metasurfaces with low densities, \rho = 2 and 5 μm^-2, where blue-green or brown hues dominate. However, bright colours emerge near critical packing at a density of 10 μm^-2, and these colours differ significantly from those observed in diffuse light. This distinction arises from the fundamentally different physical processes governing the specular and diffuse components of the scattered light. Specular reflection is formally associated with the average field in the material, which propagates in an effective homogeneous medium and is therefore subject to classical thin-film interference effects. This can lead to significant variations of the specular colour with small changes of metasurface parameters and illumination conditions. Diffuse reflection, on the other hand, is associated with the fluctuations of the scattered field around the average. Diffuse colours are mostly driven by the individual scattering elements dressed by the interaction with their neighbours, namely, for particulate metasurfaces, by dressed Mie resonances and particle pair correlations. Colour variations with the incident or scattered angles are often milder for the diffuse component than for the specular one.

Overall, the photographs in Fig. 2c and Supplementary Fig. 1 illustrate a striking diversity of visual appearances, including angular-dependent intensity, variations in lightness, whitening effects, distinct diffuse and specular colours, and changes in chroma as density and correlation vary. What makes this diversity particularly remarkable is the transition from vivid hues to softer tints near critical packing. Equally noteworthy is that this broad range of optical effects is achieved using simple, manufacturable samples on a neutral silica substrate. The fact that all samples share identical raw materials, uniform square-shaped e-beam patterns, and a fixed etching depth underscores the profound role of the nanoscale topology itself in generating such diverse appearances.

Statistical quasinormal mode analysis around critical packing

Let us now examine the cause of the sudden changes in the specular and diffuse reflectance arising as metasurfaces transition around critical packing. By analogy with solid-state systems, we consider the PDoS. Because of leakage in the air clads and absorption in silicon, the QNM frequencies are complex valued, and we thus define the PDoS as the normalised number of photon states available at a particular energy and given quality factor.

Modal analyses of waves in disordered systems often rely on point scatterers with electric dipole resonance to linearise the eigenvalue problem42. This restriction is completely relaxed hereafter by using full-field electromagnetic simulations and a supercell approach with Bloch boundary conditions in the lateral x and y-directions and perfectly matched layers (PMLs) in the z-direction (see “Methods” section). To ensure the accuracy of our predictions, we systematically test them by progressively increasing the supercell size.

Our smaller supercells consist of approximately 30 silicon metaatoms and have a lateral size of 1.2 × 1.2 μm^2. The metaatoms' side length is ‘ = 100 nm, and their height is 145 nm. All metasurface patterns are generated using a Poisson disc sampling algorithm40 for a fixed filling fraction of 0.2. The QNMs are computed and normalised using the freeware MAN (Modal Analysis of Nanoresonators)43. For computational efficiency, we assume a fixed silicon permittivity of 20 - 0.5i, an average value derived from literature data44 for the visible spectrum (\lambda = 400 to 750 nm); the negative imaginary part signifies optical loss under the e^i\omega t time-harmonic convention used by our solver. Additionally, we approximate the electromagnetically inert silica substrate as air, introducing a planar symmetry at z = 0, which reduces the computational time. Based on this assumption, we classify modes as even or odd according to the symmetry or asymmetry of the x-component of the QNM electric field. The PDoS maps are computed from a total of ≈ 7600 QNMs over 40 independent realizations, providing a statistically stable dataset. Additional QNM maps can be found in Supplementary Note 2.2.

To unify the varied structural configurations that arise from tuning the metaatom size, density, and correlation, we introduce an effective parameter, the metaatom overlap rate O_e defined as the fraction of metaatoms which are touching at least one neighbour. This fraction is calculated directly from the random pattern generated for the electron beam writer with the Poisson disc sampling algorithm. The overlap rate provides a direct measure of the structural connectivity and serves as a fundamental variable to track the system’s progression from the particulate (O_e = 0) to the aggregate regime.

Figure 3a summarises our key computational results for the PDoS maps for various topologies. To provide a comprehensive picture, each map uses two distinct colour scales: the main heatmap visualises the PDoS itself (right colour bar, logarithmic scale), while the horizontal colour bands above the maps represent the averaged QNM excitation probability (left colour bar, linear scale). The four panels show the progression through the different regimes: (i) particulate O_e = 0%, (ii) and (iii) critical packing O_e = 10% and 27% and (iv) aggregate O_e = 49%. The corresponding packing fractions are (i) p = 0.5, (ii) p = 0.2, (iii) p = 0.16 and (iv) p = 0.1.

The spectral range covered by the maps corresponds to the visible region, within which the silicon subwavelength metaatoms exhibit several distinct low-loss resonances. These include a quadrupolar-electric (EQ) mode, a degenerate dipolar-electric (EDx and EDy) mode, and a dipolar-magnetic (MDz) mode. These modes are highlighted with bright cyan squares on the maps; the x-component of their electric fields in the xy-plane are shown in the black inset in Fig. 3b. The maps in Fig. 3 reveal a clear redistribution of the QNMs in the complex frequency plane as the metasurfaces transition through critical coupling, from isolated clouds in (i) to a single, regular cloud in (iv).

PDoS diversity. For particulate metasurfaces (i), the PDoS map appears as a patchwork of distinct QNM clouds, each centred predominantly around the individual QNM frequencies. One prominent example is the high-Q (~ 80) cloud at ‘ = \lambda ~ 0.16, which primarily results from the hybridisation of MDz modes that mostly radiate in-plane45. A representative QNM from this group is Mode 1 in Fig. 3b (i). This mode exhibits a highly ordered field distribution localised in air, which is reminiscent of leaky guided resonances at k = 0 in grating waveguides46. Additional QNM clouds emerge between the real frequencies of the individual modes due to hybridisation between different mode types—for instance, the cloud spanning 0.17 < ‘ = \lambda < 0.2 is driven by the interaction between EDx and MDz modes.

A distinguishing feature of the particulate PDoS map is the emergence of pronounced pseudogaps, with a notable example at ‘ = \lambda ~ 0.21. These gaps arise from destructive interference among waves scattered by the strongly correlated metaatoms, analogous to the formation of photonic band gaps in crystals. The analogy is reinforced by the mode profiles: QNMs just below the gap localise their fields primarily within silicon (as in dielectric bands), whereas those just above the gap concentrate their fields in surrounding air (as in air band)45.

At the other end, for aggregate topologies (iv), metaatoms organise into extended chains—comparable in size to the incident wavelength—that merge into irregular clusters. Due to substantial polydispersity in cluster size and morphology, the resulting QNMs are more sparsely and almost uniformly distributed across the visible spectrum. This leads to a single, broad QNM cloud with smaller Q-factors, in agreement with the fact that larger metaatoms are better radiators. This broad cloud contrasts sharply with the fragmented PDoS of particulate topologies and arises from the complete closure of the pseudogaps—a hallmark of the profound reorganisation of the electromagnetic mode structure at critical packing. Typical QNM field profiles are shown in Fig. 3b (iv) as Modes 4 and 5. Their electric-field distributions differ markedly from those of isolated Si nanopost resonances: instead of being confined within the posts, the fields form ‘hot spots’ predominantly along the boundaries of the irregular Si clusters.

At critical packing, (ii) and (iii), the PDoS exhibits characteristics intermediate between the particulate and aggregate regimes. Even within the particulate topology, hybridised modes can shift by several linewidths from the resonance frequencies of isolated Si metaatoms, signalling the onset of strong inter-metaatom coupling. As the metaatoms begin to overlap, this coupling becomes the dominant factor shaping the PDoS, as seen in panels (ii) and (iii). In (ii), additional QNMs appear at lower frequencies, reflecting the formation of larger effective scatterers. In (iii), new hybridised modes emerge to fill the pseudogaps, marking the breakdown of the destructive interference that defines the particulate state.

Collectively, these features identify critical packing as the transitional regime where the PDoS undergoes a fundamental reorganisation. The nature of this transitional regime is exemplified by Mode 3 in Fig. 3b (ii), which exhibits a hybrid character: on isolated metaatoms (left side of the unit cell), the internal fields show a high degree of phase coherence, with each particle exhibiting a similar field profile. In contrast, when the metaatoms begin to connect, this local order breaks down, and the field profile becomes non-uniform, signalling the onset of strong, disorder-induced hybridisation.

QNM excitation probability. While the appearance of the PDoS map is informative, it is equally important to identify which modes are preferentially excited. The QNM excitation probability densities are represented by the colour bands above the PDoS maps in Fig. 3a and as colour dots in Supplementary Figs. 2.2–2.4, ranging from blue (weak excitation) to red (strong excitation). These probability densities are computed by averaging the QNM excitation coefficients over narrow spectral intervals (see “Methods” section for the detailed procedure), thereby highlighting the spectral regions of the PDoS that are most effectively excited.

A striking feature of the particulate topology (i) is the strongly bimodal nature of the excitation distribution: only QNMs within two narrow frequency bands, near the EDx and EQ resonances, have significant excitation. At critical packing (ii) and (iii), this bimodality persists, though the peaks broaden and the maximum excitation strength diminishes. In the aggregate regime (iv), the excitation distribution becomes more heterogeneous, forming a multicoloured pattern that reflects generally weak-to-intermediate excitations across a broader frequency range. This trend aligns with the equipartition theorem47,48, which suggests that, in highly disordered or chaotic systems, energy becomes uniformly distributed among all accessible degrees of freedom.

The near-field distributions in Fig. 3b provide further insight into these excitation trends. High cooperativity, where many metaatoms exhibit induced polarisations with similar phases, lead to the strong, coherent scattering required for efficient excitation by the incident plane wave. This is clearly seen in the highly excited Mode 2. Conversely, the loss of long-range coherence in the aggregate regime, exemplified by the disordered field patterns of Modes 4 and 5, results in weaker, less efficient excitation. The hybrid nature of the critical packing modes, such as Mode 3, naturally corresponds to the intermediate excitation strengths observed during the transition.

Supplementary Note 2.3 examines the odd-symmetric modes, where E_x(z) = -E_x(-z), and similarly identifies the critical packing transition as a pivotal point in the evolution of mode structure and excitation behaviour.

Phase-transition diagram. The appearance of the metasurfaces in Fig. 2 and the PDoS in Fig. 3 both suggest that the transition occurs when the metaatoms start to connect, irrespective of the specific trajectory chosen to transition, whether by varying \rho, p, or ‘. Consequently, the metaatom overlap rate O_e is expected to be a universal, trajectory-independent parameter for tracking the phase change. Figure 3c confirms this hypothesis. It provides a phase diagram that delineates the boundaries between particulate (blue) and aggregate (red) topologies across the full \rho - p - ‘ morphological parameter space. The critical-packing regime is highlighted by two constant-O_e contours, corresponding to 10% (dashed) and 30% (solid) overlap rates. This phase diagram underscores that the critical-packing regime is trajectory independent and emerges as a general feature governed by the combined effects of all the morphological parameters.

To further reveal the impact of the transition in the PDoS maps, we further introduce the order parameter, \Phi, defined as the largest connected–component fraction. This metric, widely used in the statistical analysis of complex networks49, measures the degree of fragmentation: \Phi ranges from 0 to 1, with small values indicating highly fragmented networks.

Applied to probability density functions such as the PDoS, \Phi is obtained by binarising the PDoS maps and computing the ratio of the largest cloud area to the total area occupied by all clouds. In other words, \Phi represents the fraction of “active” PDoS pixels that belong to the largest connected cloud49. Figure 3d plots \Phi as a function of the overlap rate O_e, showing a sigmoidal evolution toward unity, with a sharp transition near O_e = 20%. This evolution marks a shift from fragmented QNM clouds separated by band gaps to a single, spanning QNM cloud that extends across the entire visible spectrum. The coloured insets along the curve highlight this evolution, illustrating the largest connected region in each map.

Intuitively, as adjacent metaatoms begin to connect, the likelihood increases that additional metaatoms will also link up, leading to the rapid formation of larger clusters. This, in turn, results in a significant shift in PDoS connectivity. Thus, although Fig. 3a, d follow one specific trajectory in parameter space (with only p varied), the abrupt transition in PDoS connectivity—and the associated pseudogap closure or optical response—should be considered a universal feature. This transition is expected to occur for any path that passes through the critical-packing regime. This universality is likely to hold when the optical wavelength is comparable to the metaatom size (\lambda ~ n‘, where n is the refractive index of the metaatoms), although the precise visual appearance (e.g. hue or saturation) may still depend on the specific trajectory taken.

Blue shift of diffuse light in dense particulate metasurfaces

Strong modal coupling between localised resonances is a key factor driving the widespread use of subwavelength metaatoms46, particularly when their optical near fields strongly overlap50,51. This principle is well established and widely exploited in ordered systems—such as small clusters of nanoresonators or periodic arrays—to finely tune the structural colours arising from collective resonances52,53.

In contrast, the behaviour of disordered ensembles of nanoresonators remains less well understood. Due to the random spatial arrangement, the polarisations induced by multiple scattering fluctuate and tend to statistically average out23,54, weakening coherent interactions. Consequently, strong coupling effects in disordered metasurfaces often manifest more subtly54,55.

Figure 4 highlights a striking and unexpectedly blue shift that stands out as particularly remarkable in comparison. The shift is displayed for an incidence angle of 30°, with three ping-pong ball images obtained for p = 0.5 and increasing densities. At a low density of \rho = 2 μm^-2, multiple scattering is minimal. The diffuse light appears cadmium green, with an angle-dependent hue marked by a suppressed diffusion near the specular reflection direction and a bright halo marked with a dashed arc line. These features—specular suppression and halo enhancement—are not central to the present discussion and are addressed in Supplementary Note 1.

As the density increases to 5 μm^-2, the halo shifts to larger angles. The shift is accompanied by a modest increase in brightness and the onset of a perceptible bluing effect. At a high density of 10 μm^-2 just below the critical packing threshold, a vivid sapphire blue emerges. Further increases in density would lead to a whitening of the diffuse light due to metaatom overlapping, as discussed in the next section.

Given that the shape and size of the metaatoms are precisely controlled during fabrication, the blue shift is attributed to direct electromagnetic interaction between metaatoms, particularly through evanescent and radiative couplings that are enhanced at higher densities. In most studies of cooperative scattering in particulate systems, particles are treated as point scatterers with resonant electric polarizabilities, where multiple scattering is modelled through dipole-dipole interactions, as seen in research on cold atomic clouds56,57 and micro- and nanoscalediscrete disordered media58,59.

Instead, our high-index, 145 nm-tall metaatoms support two electric dipolar modes, two magnetic dipolar modes, and one electric quadrupolar mode within the visible spectral range, see Fig. 3b. While it is possible to incorporate these multipoles using a closed set of equations60, deriving a closed-form expression for diffuse light intensity that accounts for both multiple scattering and multipolar interactions is notoriously challenging, as it requires solving for the field-field correlation function21.

To model the blue-shift effect, we employ a recent model that incorporates multiple scattering effects in a mean-field manner34,36, without requiring explicit consideration of the multipolar response. In this model, the diffuse contribution to the BRDF is expressed as:

d\sigma_{diff} / d\Omega_s = \rho (d\sigma_s / d\Omega_s) S(r) C(k_s, e_s, k_i, e_i), (1)

where the first three terms—the density, the form factor d\sigma_s / d\Omega_s which describes how a single isolated metaatom scatters light and the structure factor S(r) —neglect the electromagnetic interaction between metaatoms. Only the final term, the correction factor C, accounts for this interaction and is therefore responsible for the observed blue shift. In Eq. (1), k_i and k_s denote the incident and scattered wave vectors, and e_i and e_s represent the corresponding polarisations.

Without delving into the details provided in Supplementary Note 3, the correction factor is theoretically inferred in the model from a mean-field estimation of the field driving every scatterer in the presence of multiple scattering. Through energy conservation arguments, it is possible to show that the sum of the specular transmittance (T_s) and reflectance (R_s) plays a key role in inferring the mean field. In ref. 34, these specular coefficients were inferred theoretically using advanced homogenization theories61.

Here, we instead use measured values of T_s and R_s, improving resilience to topological and material variations and minimising potential inaccuracies in the homogenization theories. The spectral and angular maps of (1 - T_s - R_s), used to infer the correction factor, are shown in Supplementary Fig. 3a. They are recorded with a 2-degree step variation in the incident angle and averaged over TE and TM incident light. Note that the maximum of (1 - T_s - R_s) is blue-shifted for the denser metasurface for \rho = 10 μm^-2, compared to the lower densities. Consequently, the correction factor peaks near this wavelength (Supplementary Note 3), suggesting that the averaged local field exciting the metaatoms is maximised due to electromagnetic interaction.

To quantitatively assess the model accuracy, we measure the BRDF of the metasurfaces and compare the experimental data with theoretical predictions from Eq. (1). The intermediate steps for calculating the correction factor and subsequently deducing the BRDF are detailed in Supplementary Fig. 3, which additionally includes colour maps of C(k_s, e_s, k_i, e_i) as functions of wavelength and scattering angle.

The agreement between theory and experiment is notably good. The model effectively captures the BRDF increase and the resonant blue shift at high density. Remaining discrepancies likely stem from inaccuracies in modelling the metaatom shape and size, as well as moderate variations in these parameters, which significantly influence metaatom resonances.

The pronounced blue shift observed here is notably different compared to previously reported coupling effects induced by multiple scattering in high-density disordered metasurfaces54,55. This shift points to a novel mechanism for harnessing structural colouration in disordered metasurfaces, paving the way for innovative design strategies that move beyond conventional periodic architectures.

Whitening effect across the critical packing transition

The previous effect arises for metasurfaces with tightly packed non-overlapping metaatoms. As the fabrication parameters are further tuned, the critical packing threshold is reached, and the diffuse light shows important hue variations.

These variations are illustrated in Fig. 5a for a fixed metaatom size ‘ = 95 nm. Two transitions from particulate to aggregate topologies are highlighted, either by decreasing the packing fraction from p = 0.5, to 0.3 and then 0.1 (transition (i)) or by increasing the density from \rho = 2, 5 and then 10 μm^-2 (transition (ii)), causing the non-iridescent blue and green diffuse colours of the ping-pong ball images to converge into a white hue. The whitening is confirmed by in-plane BRDF measurements showing the spectral and angular dependence of the diffuse light for an incident angle of 30° (Fig. 5a). As we cross critical packing, the initial spectral distributions of the particulate topologies initially centred in the blue-green spectral region suddenly broaden, covering almost uniformly the entire spectral range from 420 nm to 700 nm.

To quantitatively assess the purity of the white hue of the spectrally broad BRDF, we map the ping-pong ball images onto the CIE 1931 chromaticity diagram. This diagram is widely used to correlate reflection or transmission spectra with the colours perceived by human vision. The two transitions are indicated by the trajectories (i) and (ii) in Fig. 5b, where we additionally mark the white point with coordinates (0.3127, 0.3290) with a blue circle. This point represents a standard reference for many colour applications with the D65 illuminant. Since our supercontinuum illuminant is different from the average daylighting, we also mark the white point with the D50 illuminant used in graphic arts, making it easier to evaluate the purity of the observed white, which is very close to the two marks. This trend is further quantified in Supplementary Note 6, where the 3D CIELAB representation confirms the trajectories’ convergence toward the white point in chromaticity, while also explicitly showing the marked increase in lightness (L*) that accompanies the visual effect.

The whitening phenomenon can be intuitively understood through the QNM analysis in Fig. 2. It can be attributed to an inhomogeneous broadening of aggregate topologies, where particles overlap and even cluster. Then, numerous modes, which typically cover all visible frequencies, are significantly excited by collimated beams, scattering light at all visible frequencies in all directions.

A key performance metric for white paint or high haze surfaces is the brightness, which quantifies the fraction of the light that is effectively diffused. High brightness values are typically achieved with materials that exhibit strong scattering, often involving complex multilayer structures or thick particulate coatings21,62. Figure 5c presents the spectral dependence of the diffuse light—actually, 1 - T_s - R_s includes ~ 15% Si absorption in the blue—for the two metasurfaces that generate the most pronounced whitening. Both metasurfaces have minimal correlation (p = 0.1) and maximum density (\rho = 10 μm^-2), with the only difference being the metaatom sizes, ‘ = 95 and 170 nm. The metasurface with the larger metaatom size forms larger aggregates, resulting in a brightness approaching 80% across the entire visible spectrum.

While this value is lower than the brightness achieved by certain animals63 or ultra-white paints64, the 145-nm-thick metasurfaces strike a relevant balance between brightness and thickness. This makes them a promising choice for compact, photonic CMOS-compatible devices in applications like materials science and lighting, particularly for reducing glare in electronic displays.

Broadband mirror effect at critical packing

The two preceding effects are observed for diffuse light. Despite their stark contrast in nature, both diffuse and specular lights are influenced by the PDOS and the QNM excitation probabilities. Thus, much like diffuse light, specular light is expected to exhibit rapid variations at critical packing.

This expectation is confirmed in Fig. 6, where we analyse the specular reflection spectra of five metasurfaces transitioning from particulate to aggregate topologies by increasing the silicon post size and disorder correlation. As expected, we first observe a colour transition, from forest green to a dark shade, passing through a distinct scarlet red at critical packing.

Further insight is gained from the reflectivity spectra in Fig. 6b (i). For the particulate metasurface with the smallest metaatoms, the specular reflectivity R_s remains weak across the visible spectrum, peaking below 10%. This metasurface primarily transmits light coherently, with a specular transmission T_s > 80% for \lambda > 500 nm. Similarly, the aggregate metasurface with minimal correlation exhibits weak specular reflectivity. Here, the polydispersity of cluster sizes and shapes, combined with their weakly correlated positions, suppresses coherence, leading to dominant diffuse scattering across all wavelengths (Fig. 5).

This dual behaviour explains why aggregates can appear bright white in diffuse reflection while remaining dark in the specular direction. The broad continuum of resonant states ensures efficient scattering of all visible wavelengths into diffuse channels, giving rise to whiteness. At the same time, the randomised phases of these scattered waves suppress specular reflection. The effect is robust across particle sizes large enough to scatter efficiently in the visible, as illustrated in the last three rows of the third column of Fig. 2c.

Between these extremes, the intermediate metasurfaces exhibit significant changes in specular reflectivity as they approach critical packing. The peak reflectivity first increases, reaching 50% in the red, before declining. Notably, the ratio R_s / T_s at critical packing is particularly high, reaching ~ 15 (Fig. 6b (ii)).

A large coherent reflection of 50% is a singular observation. It is neither explained by existing models nor intuitively anticipated by considering the morphology of the critical-packing metasurface in Fig. 6a. To contextualise this observation, we compare the metasurface reflectance with that of a benchmark highly coherent system: a Si film on glass with the same 145 nm thickness. The Si-film reflectance computed using tabulated refractive index data of silicon44 are shown with dashed curves in Fig. 6b (ii–iii). It exhibits different characteristics, such as an absence of a marked peak for R_s / T_s and a reflectance nearly independent of the incident angle. Interestingly, the maximum reflectance of 60% is only slightly higher than that of the metasurface, highlighting the singular capability of the critical-packing metasurface to coherently reflect light despite its strong heterogeneities.

No robust theoretical framework currently exists for modelling metasurfaces near critical packing21,38,46,59,65,66. However, to achieve a qualitative insight into the singular properties of the critical packing metasurface, we attempted to infer effective parameters using effective medium theory46. For the two particulate metasurfaces, predictions based on an Extended–Maxwell–Garnett model67 align well with specular reflectance measurements (Supplementary Fig. 4), by using the electric and magnetic dipole Mie coefficients of the individual metaatoms. As expected, the model is less quantitative at critical packing, where it significantly underestimates the reflection peak at 670 nm (Supplementary Fig. 4). While limited confidence can be placed in the model within this regime—where it exceeds its range of validity—we note an increasing role of the magnetic dipole resonance as the metaatom sizes increase to reach critical packing, where we find that the reconstructed permittivities and permeabilities are both remarkably close to zero at the peak wavelength.

The QNM analysis furthers support the significance of magnetic effects. At critical packing, the excitation probability distribution for even-symmetry QNMs in Fig. 3a (ii) peaks around the electric dipolar modes of individual metaatoms, while the odd-symmetry QNM PDOS in Supplementary Fig. 2.5 peaks from magnetic dipolar modes. Since both figures indicate a single dominant QNM with comparable excitation efficiency—such as mode 3 for the even-symmetry QNM in Fig. 3a (ii)—we expect that the metasurface at critical packing exhibits both significant effective permeability and permittivity.

Methods
Sample fabrication

The metasurfaces were fabricated by etching structures in a 150 nm-thick polycrystalline silicon (Si-poly) layer. First, Si-poly was deposited on both sides of a 4-inch fused silica (FS) wafer by LPCVD at 605 °C. Then, Si-poly was removed from one face of the FS wafer using fluorine-based based reactive ion etching (F-RIE). The sample was then covered by a 160 nm-thick layer of a negative resist (maN 2405—Micro Resist Technology) followed by a 40 nm-thick conductive layer (ELECTRA 92—All Resist) to avoid charging effects during subsequent electron beam lithography.

Samples were exposed by EBL with a beam energy of 20 keV at a current of 135 pA with a step size of 4 nm and an nominal dose of 150 μC/cm^2. Resist was developed using MF-CD-26 (Micropposit) at 20 °C for 50 s, rinsed with deionised water for 1 min and dried with nitrogen. Negative maN resist patterns were transferred in the Si-poly layer using F-RIE, and the remaining resist mask was removed by oxygen plasma. Depending on density and applied dose, one can vary the size of the squares in the resist after development using proximity effects. Increasing the dose can then lead to fusion of squares as it can be seen in Figs. 1, 2, 5 and 6.

Visual observation with the ping-pong ball

A 200 μm supercontinuum laser beam (Leukos, Rock 400) is slightly focused onto the metasurface samples. A semi-translucent half-ping-pong ball is positioned at the centre of the metasurface, serving as a hemispherical screen for the diffusely reflected light. The images of the ping-pong ball are recorded with a Canon EOS 1000D camera mounted on an optical table with fixed exposure times for each series. Full details on the characterisation of diffuse and specular light, including the diffusive transmission properties of the ping-pong ball, are provided in Supplementary Note 5.

BRDF characterization

BRDF measurements are conducted using a custom-built goniospectrophotometric setup equipped with a supercontinuum laser and a 1 mm-diameter optical fibre connected to a spectrometer (Ocean Insight, HDX). The incidence (\theta_i) and scattering \theta_s, \phi_s angles are precisely controlled by three stepper motor rotation stages (Newport URS75, URS150, and SR50CC).

To ensure measurements in the visible range, the unpolarized supercontinuum laser beam is filtered using a short-pass filter (Schott KG-1). The incident laser radiant flux is measured within the same setup, with the fibre detector aligned to capture the focused laser beam.

Electromagnetic computational results

The QNMs are computed for a null in-plane wavevector with the QNMEig solver of the freeware MAN43. This solver automatically identifies a large number of modes near a user-defined complex frequency. The computed eigenvectors include both QNMs and numerical PML modes. To eliminate PML modes, we implement a systematic procedure detailed in Supplementary Note 2.1. After filtering, we normalise the QNMs43,73 and rank them based on their excitation coefficient43

|\alpha_m|^2 = |\int \int \int (\epsilon_{Si} - \epsilon_b) Q_m E_b \cdot E_m d^3r|^2, (2)

for a normally incident plane wave at a frequency matching the real part of the QNM frequency. Here, \epsilon_b = 1 and \epsilon_{Si} = 20 - 0.5i denote the background and silicon relative permittivities, respectively. The negative sign of the imaginary part reflects COMSOL’s use of the e^i\omega t time-harmonic convention, where a negative imaginary component corresponds to material loss. The spatial overlap integral between the incident electric field E_b and the normalised QNM electric field E_m is evaluated within the silicon posts of the supercell.

Figure 3a using even-symmetric QNMs collected from 40 independent metasurface realisations. Twenty of these 40 realisations are displayed in Supplementary Figs. 2.2−2.4, where each dot represents a QNM complex frequency, with the colour encoding the excitation coefficients (blue for low values, red for high). In each realisation, only a small subset of QNMs is strongly excited. On average, approximately 190 QNMs are computed per realisation, resulting in a total of around 7600 QNMs used to estimate the PDoS maps. Additional realisations for odd-symmetric modes at critical packing density are presented in Supplementary Fig. 2.5, with the corresponding PDoS map shown in Supplementary Fig. 2.6.

Importantly, note that consistent results are observed for supercells four times larger (Supplementary Figs. 2.7 and 2.8), reinforcing our confidence that the statistical data from the 1.2 × 1.2 μm^2 supercells are sufficiently converged to support our conclusions on the transition of the PDoS across particulate, critical-packing, and aggregate topologies.

Data availability

The data that support the findings of this study are available within the article and its Supplementary Information. All data generated from this study are available from the corresponding authors upon request.

Code availability

The COMSOL model for computing the QNMs will be provided in the freeware MAN (Modal Analysis of Nanoresonators) V9 and the following versions.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript presents an interesting study on color generation using disordered metasurfaces. To the best of my knowledge, the concept of “critical packing” is proposed here for the first time and is shown to play an important role in controlling the color output. I find the work novel and impactful, and well-aligned with the high standards of Nat. Comm. However, before I can recommend acceptance, I would like the authors to address the following concerns:\nMajor Comments\nA. Material Dependence and Generality\nThe study focuses on dielectric particles without localized gap modes.\nI am curious about whether similar phenomena could be observed using metallic nanoparticles with the same topology.\nAdditionally, since the authors employed high-index silicon, what would happen if lower-index materials (e.g., silica) were\nused?\nIn other words, how sensitive is the observed effect to the choice of material?\nB. Shape and Topology Variations\nNanocubes were used as the meta-atoms in this work.\nWhat would be the expected impact if nanodisks or nanospheres were used instead?\nI am particularly curious about the role of both particle morphology (shape) and connectivity (topology) in the observed\neffects.\nC. Unified Parameterization\nThe authors use three parameters—density, size, and disorder correlation—to characterize the topology.\nIs it possible to define a single effective parameter (e.g., T) that captures the essential physics of the transition?\nIt would be valuable to make a quantitative demonstration of the transition, similar to a phase diagram. Is that possible?\nD. Density of States in Fig. 3a(i)\nNotably, this redistribution eliminates the localized bandgap-like regions previously observed in the particulate case at\nleading to a continuous PDoS across the entire visible range.\nIn Fig. 3a(i), there is a black gap around 0.21. It is unclear why the DoS is particularly low at that frequency. A clarification\nwould help.\nE. Color Range in Aggregates\nIt is quite interesting that both black and white colors can be achieved in the aggregated cases.\nI recommend a more detailed discussion on this point.\nParticle size plays a key role?\nAchieving both black and white from the same material (Si) is hard in conventional systems, as far as I know.\nMinor Comments\n1.Fig. 2 – Color Transition\nThe \"sudden changes\" in color mentioned in the text are not immediately obvious from Fig. 2.\nHowever, Figs. 4 and 5 provide clearer insights. A forward reference in the text noting that this will be explained in more\ndetail later would help guide the reader.\n2.Fig. 2c – Black Panels\nWhy are there three black panels in row 2, column 1 of Fig. 2c?\nThis anomaly is not explained in the figure caption or main text.\n3. Excitation Probability Densities\nThe method used to calculate excitation probability densities is not described.\nA brief explanation or a reference would improve reproducibility.\n4. Fig. 5b – Intensity Information\nThe CIE color map in Fig. 5b lacks intensity information. Since the configuration impacts both color and brightness, I suggest\nadding a separate intensity plot or using a 3D color space to more comprehensively show the effect.\n5. ref recommendation\nIn solid-state systems, rapid changes in the electronic density of states (DoS) are often associated to spectacular\nphenomena, such as phase transitions for instance.\nConsider the broad readership of Nat Comm, some reference is recommended.\n6. Radiative Cooling Discussion\nThe manuscript compares the brightness of the metasurfaces to that of white paints used for radiative cooling.\nHowever, radiative cooling relies not only on high visible/NIR reflectivity but also on high emissivity in the mid-infrared (~10\nµm).\nI suggest briefly discussing this aspect, and whether Si is suitable for such applications.\n7. Permittivity of Silicon\nIn Fig. 3, the real part of the permittivity for silicon is set as 20.\nHowever, in Phys. Rev. B 27, 985–1009 (1983), it is reported to be 16.7.\nIs the increase due to quantum confinement effects?\nA reference or short explanation would be helpful.\n8. Additional References on Structural Coloration\nThe manuscript mentions recent developments in coloration technologies. I recommend including the following relevant\nreferences to strengthen the background and highlight related progress:\nMao, Peng, et al. \"Manipulating disordered plasmonic systems by external cavity with transition from broadband absorption\nto reconfigurable reflection.\" Nature Communications 11.1 (2020): 1538.\nCencillo-Abad, Pablo, et al. \"Ultralight plasmonic structural color paint.\" Science Advances 9.10 (2023): eadf7207.", "reviewer_2": "This work analyses the light scattering properties of disordered metasurfaces made of silicon nanoparticles (metaatoms) dispersed on a glass substrate. By controlling the particle size, the surface density and the packing fraction (i.e. the range of effective “hard-sphere” correlations), the authors measure the BRDF and the specular reflection in different regimes. They focus on changes in the visual appearance in regimes ranging from isolated particles to aggregates, and identify the intermediate regime in which metaatoms begin to connect as “critical packing”. They put forward the peculiarities of this regime by different theoretical analyses, and show its potential in the design of thin coatings with substantial impact on visual appearance.\nThe work is original, and of potential high interest in the context of metasurfaces for the control of light waves in the visible spectrum. I have some concerns, that are summarized below (major points and minor points). They need to be addressed by the authors before any recommendation can be made. At this stage, I cannot recommend publication in Nature Communications.\nMajor points\n1) A huge literature has emerged in the 1980s on metal (gold or silver) semicontinuous films, with studies on colors and scattering in the far field and in the near field. The authors seem to consider from the outset that their systems are very different, and that there is no need to draw on the very many previous studies on metal particle film. I find this a bit hasty. In other words, the context should be revisited.\nStatements like “Unlike percolation metal films, this transition affects not only the specular but also the diffuse components of the scattered light in a profound way” or “Unlike percolating metal films, the transition observed in our system influences not only absorption, but also both specular and diffuse optical responses” or “in our dielectric metasurfaces, we do not observe the localized gap modes that play a key role in semicontinuous metal films” seem to me too superficial. Is it true that diffuse scattering of semicontinuous films has not been studied ? Is it clear that their properties are influenced only by absorption ? Is it clear that the localized gap modes are responsible for most of their properties ? I would expect a more careful analysis of similarities/differences since comparing both systems might be natural to many readers.\n2) The spatial correlations between particles induced by the hard-sphere type potential (exclusion region) are mentioned as an important degree of freedom. I would have expected to see a characterization of the degree of correlations other than the packing fraction p. It could be useful to show the structure factors of some key configurations. The structure factor appears in Eq. 1, and is briefly mentioned in the SI, but its role is largely overlooked.\n3) The word “percolation” hardly appears the manuscript. The authors may have a good reason for that, but this should be briefly explained. A (probably naïve) approach could be to see the critical packing regime as an intermediate regime between a highly correlated system with isolated particles and an aggregate system that could percolate at large density. Whether a point of view based on 2D percolation (another widely explored field) is relevant or not should be explained.\n4) The analysis based on the photonic density of states (PDoS) seem to show that the statistics is important on revealing the different regimes. I am not sure to understand why the authors do not study the fluctuations (variance of distributions) of the PDoS both spatially and spectrally. This could be an interesting observable (see for example Krachmalnicoff et al., PRL 105, 183901, 2010) and would make the statistical analysis more substantial.\nMinor points\n5) Introduction, line 3: The authors write “the analogue of the electronic DoS is played by the statistics of the photon density of states”. Why “statistics” at this stage ?\n6) Page 3, paragraph starting with “The 27x2 photographs”: It seems that the incidence angle for specular and diffuse measurements are not consistent with the caption of Fig. 2. This needs to be adjusted.\n7) Page 5, paragraph starting with “specular colours tend”: I am not sure to understand the statement\n« This distinction arises from the fundamentally different physical processes governing diffuse and specular reflection. In diffuse reflection, light experiences the granular nature of the materials and is scattered in multiple directions. The material microscopic resonances—such as the properties of individual scatterers—play a central role in this process, while thin-film interference effects are relatively minor. In contrast, specular reflection is driven by collective phenomena at a mesoscale ».\nThe usual understanding is that specular reflection is associated with the average field, while diffuse reflection with the field fluctuations. In this setting, I do not see how the processes can be discriminated between “individual” and “collective”. The authors should clarify their point.\n8) The choice of references is sometimes surprising. For example, Page 8, bottom line, Ref. 47 is not about atomic clouds, Ref. 51 is not about discrete media…. There are many inconsistencies in references all along the manuscript. The authors should systematically check the relevance of their references, and use either the first reference in a field, or a review paper, or a textbook… and not references that may seem to have been chosen at random."}, "author_response": "Response to the referees\nWe sincerely thank the referees for the time and effort they dedicated to reviewing our manuscript and for\ntheir valuable comments, which have greatly improved our work. All points raised have been carefully\naddressed and discussed below, along with the corresponding modifications made to the text.\nReferee #1\nReferee report\nThe manuscript presents an interesting study on color generation using disordered metasurfaces. To the\nbest of my knowledge, the concept of “critical packing” is proposed here for the first time and is shown to\nplay an important role in controlling the color output. I find the work novel and impactful, and well-aligned\nwith the high standards of Nat. Comm. However, before I can recommend acceptance, I would like the\nauthors to address the following concerns:\nOur response\nWe thank the reviewer for summarizing our work and highlighting its novelty and importance.\nSpecific comments and questions and our step-by-step answer\n1. Material Dependence and Generality\nThe study focuses on dielectric particles without localized gap modes. I am curious about whether similar\nphenomena could be observed using metallic nanoparticles with the same topology. Additionally, since\nthe authors employed high-index silicon, what would happen if lower-index materials (e.g., silica) were\nused? In other words, how sensitive is the observed effect to the choice of material?\nReply:\nHow our results are affected by material choice is indeed an important question we have asked ourselves.\nIn short, we think that they are very sensitive to material choice. It is difficult to quantitatively answer\nwithout a comprehensive study. However, we can reasonably say that:\nMetallic Metasurfaces (e.g., gold or silver): As noted by the reviewer, dielectric system does not support\nthe localized gap-plasmon modes that are characteristic of closely-packed metallic nanoparticles. Close to\ncritical packing, the electromagnetic fields of metal topologies would support a bunch of additional QNMs\nwith intensely localized fields in the nanogaps between particles. These modes are notoriously lossy and\nbroadband. On top of that, metallic particles are also more absorbent than dielectric ones. We can\nreasonably anticipate that the transition at critical packing is likely to be dominantly impacted by absorption,\nreducing the visibility of effects observed close to critical packing with our dielectric metasurfaces, e.g. blue-\nshift and whitening effects.\nLow-Index Dielectric Metasurfaces (e.g., SiO₂ or polymers): The use of high-index silicon is critical for many\n(all?) effects we report. For low-index contrasts, the resonances would be significantly weaker and much\nbroader (i.e., have a very low Q-factor), confinement within the particles would be poor, and colours would\nbe less vivid, and the sudden transition we report would be largely absent. The sharp spectral features that\nunderpin the multiple-scattering-induced blue shift and the dramatic changes at critical packing would be\nwashed out. While some variation in scattering would occur with changing density, the sudden and striking\ntransitions we report would be largely absent.\nIn summary, the phenomena of a multiple-scattering-induced blue-shift, a distinct whitening effect, and a\nbroadband mirror at critical packing are intrinsically linked to the use of a low-loss, high-index dielectric\nmaterial. This material choice supports strong, spectrally-distinct scattering resonances without the\noverpowering absorptive losses of metals.\nAction taken:\nTo make these important points clearer in the manuscript, we have added a paragraph in\nthe Outlook section addressing the expected behaviours of metallic and low-index dielectric systems. We\nrepeat the new paragraph here for clarity:\n“The use of high-index-contrast Si subwavelength metaatoms with a few distinct low-loss resonances is central in the\npresent work. By comparison, low-index-contrast metaatoms, such as those fabricated from polymers, are expected to\nexhibit much weaker optical effects. Metallic metaatoms, on the other hand, behave quite differently: at critical packing,\nthey support numerous additional resonances with intense fields localized in the nanogaps between nearly touching\nparticles. This produces a photonic density of states (PDoS) landscape distinct from that shown in Fig. 3a. Because these\nextra modes are lossy and span a broad spectral range, absorption strongly influences the critical-packing transition,\ndiminishing the visibility of certain effects—such as the blue shift—that we observe in dielectric metasurfaces. However,\nthese reductions may be offset by the striking emergence of metallic cluster.”\n2. Shape and Topology Variations\nNanocubes were used as the meta-atoms in this work.\nWhat would be the expected impact if nanodisks or nanospheres were used instead?\nI am particularly curious about the role of both particle morphology (shape) and connectivity (topology) in\nthe observed effects.\nReply:\nThe shape of the individual metaatoms directly determines the resonant frequencies and thus the colour.\nHowever, we believe that the precise shape will alter much the properties of the particulate and aggregate\ntopologies. In particular, the large dispersity of appearance of Fig. 2 is likely to be resilient to morphology\nvariation. Similarly, the whitening effect that stems from the dense continuum of modes created by\ndisparity between clusters in the aggregate regime will be observed with other metaatom shapes and sizes.\nThe special case of elongated particles deserves more attention, especially under the condition that their\norientation is not random. Then anisotropy will be preserved in the metasurfaces and new polarization\neffects may be anticipated for both the diffuse and specular light.\nDespite undeniable differences concerning colours, we thus hypothesize that all our observations—a blue-\nshift from multiple scattering, a whitening effect, and a specular mirror transition—would still remain for\nother morphologies. The sudden emergence of a new set of collective resonant modes as the system\ntransitions from a collection of isolated scatterers to a network of connected ones is likely to be resilient to\nmorphology variations.\nAction taken:\nIn the Outlook section of the manuscript, we have added a short paragraph which is discussing shape\nmorphology and material. The paragraph reads as follows.\n“Similarly, variations in the morphology of the metaatoms—such as changes in size or shape—do influence their\nresonance and thus the resulting colour. However, these variations do not hinder the observation of the transition\nphenomena. In particular, the pronounced dispersity of appearances in Fig. 2 is expected to remain robust against\nmorphological fluctuations. A noteworthy special case arises with elongated particles, especially when their orientation\nis not random. In such situations, anisotropy is preserved within the metasurfaces, and new polarization-dependent\neffects may emerge in both diffuse and specular light.”\nWe expect that our responses to the referee questions on morphology and material variations broaden the\nscope of our work.\n3. Unified Parameterization\nThe authors use three parameters—density, size, and disorder correlation—to characterize the topology.\nIs it possible to define a single effective parameter (e.g., T) that captures the essential physics of the\ntransition?\nIt would be valuable to make a quantitative demonstration of the transition, similar to a phase diagram.\nIs that possible?\nReply:\nWe again acknowledge an interesting and insightful remark, which has significantly strengthened the\nmanuscript. To address the reviewer's points, we have introduced two key quantities: a unified physical\nparameter, 𝑂̃, and a corresponding order parameter, Φ, and we present the associated analysis in two new,\ndedicated panels of Fig. 3.\nInspired by studies on semicontinuous metal films—where a phase transition arises as the connecting\nlength diverges at percolation—we define the metaatom overlap rate 𝑂̃, the fraction of interconnected\nmetaatoms. This parameter provides a single, unifying descriptor that captures the system’s progression\nfrom a particulate to an aggregate topology, independent of the specific trajectory in the 𝜌−𝑝−ℓ\nparameter space.\nUsing 𝑂̃, we construct a phase diagram that spans the full 𝜌−𝑝−ℓ parameter space. The diagram\ndelineates three distinct topological regimes—particulate, critical, and aggregate—through contours of\nconstant overlap. This demonstrates that the critical packing transition is not an artifact of particular\nparameter choices, but rather a general and universal feature of the system, accessible through multiple\ncombinations of structural parameters.\nTo further establish the transition’s character, we introduce an order parameter, Φ, defined as the largest\nconnected–component fraction of the PDoS. In the new analysis, Φ exhibits a sharp, sigmoidal increase as\na function of 𝑂̃. This provides robust quantitative evidence for the abrupt, phase-transition-like nature of\nthe phenomenon and highlights the central role of 𝑂̃ as the governing parameter.\nAction taken:\nWe have revised Fig. 3 to directly address this point:\n1． Panel (a) now shows a sequence of PDoS maps illustrating the evolution of the modal structure\nwith increasing 𝑂̃.\n2． Panel (c) presents the phase diagram of the 𝜌−𝑝−ℓ parameter space, with contours of constant\noverlap delineating the particulate, critical, and aggregate regimes.\n3． Panel (d) provides the quantitative demonstration of the transition, plotting the order parameter\nΦ against 𝑂̃. Insets illustrate the evolution of the largest connected component of the PDoS map.\nThe section Statistical quasinormal mode analysis around critical packing has been substantially rewritten.\nIt now includes a new, dedicated subsection titled \"Phase-transition diagram\" that formally introduces and\ndiscusses the unified parameter 𝑂̃, the phase diagram in Fig. 3c, the order parameter Φ, and the\nquantitative results in Fig. 3d.\nFor the reviewer’s convenience, the revised Fig. 3 is reproduced below as Fig. R1.\nFigure R1. Evolution of the photon density of states (PDoS) when transitioning through critical-\npacking topologies. All metasurface patterns are generated using a Poisson disc sampling\nalgorithm for Si metaatoms with lateral size ℓ=100 nm and height of 145 nm. The filling\nfraction is 0.2. a PDoS maps of even electromagnetic QNMs for increasing overlap rate 𝑂̃: (i)\nparticulate (𝑂̃ =0% ), (ii–iii) critical packing (𝑂̃ =10% and 27% ), and (iv) aggregate (𝑂̃ =\n49%). The corresponding packing fractions are 𝑝 =0.5, 0.2, 0.16 and 0.1, respectively. Colour\nbands above each map represent the averaged QNM excitation coefficients (see Methods),\nwhile cyan squares mark the dominant single-metaatom modes—MDz, EDx, and EQ. b\nRepresentative field maps showing the 𝑥-component of the QNM electric field (see Fig. S2.8 for\ndetails). The black inset shows the MDz, EDx and EQ QNM fields in the 𝑥𝑦-plane. c 𝜌−𝑝−\nℓ phase diagram distinguishing particulate (blue) and aggregate (red) topologies. d Quantitative\nanalysis of the PDoS evolution shown in a, from fragmented clouds at small 𝑂̃ values to a unique\nspanning cloud at Φ≈1. Insets show the largest connected component for every PDoS map.\n4. Density of States in Fig. 3a(i)\nNotably, this redistribution eliminates the localized bandgap-like regions previously observed in the\nparticulate case at leading to a continuous PDoS across the entire visible range.\nIn Fig. 3a(i), there is a black gap around 0.21. It is unclear why the DoS is particularly low at that frequency.\nA clarification would help.\nReply:\nWe thank the reviewer for this comment on the gap around ℓ/𝜆 ≈ 0.21 in Fig. 3a(i).\nThis dip originates from destructive interference in the particulate regime: individual nanoposts act as\nstrong scatterers with well-defined resonances (notably the EDx mode just below this frequency). When\narranged with strong short-range correlations (p=0.5), their scattered waves interfere destructively,\nsuppressing the formation of quasinormal modes in this frequency range. The result is a pseudogap,\nanalogous to a photonic bandgap in a crystal but arising here in a disordered system.\nAs shown in Fig. R3, the analogy is reinforced by the field profile of modes: just below the gap the field\nresides mainly inside the silicon posts, while just above the gap it shifts into the surrounding air — closely\nmirroring the dielectric/air localization seen across photonic crystal gaps. In the aggregate regime, this\npseudogap disappears because the posts merge into irregular clusters of varying shapes and sizes, creating\na continuum of hybridized modes that fill in the DOS across the spectrum.\nTherefore, the presence of this gap in the particulate topology and its disappearance at the transition to\nthe aggregate topology is another key piece of evidence for the dramatic change in the underlying\nelectromagnetic mode structure that we define as critical packing.\nAction taken: To ensure this important point is clear to all readers, we have amended the manuscript to\nexplicitly explain the origin of this pseudogap and its significance for the transition.\n“A distinguishing feature of the particulate PDoS map is the emergence of pronounced pseudogaps, with a notable\nexample at ℓ/𝜆≈0.21. These gaps arise from destructive interference among waves scattered by the strongly\ncorrelated metaatoms, analogous to the formation of photonic bandgaps in crystals. The analogy is reinforced by the\nmode profiles: QNMs just below the gap localize their fields primarily within the silicon (as in dielectric bands), whereas\nthose just above the gap concentrate their fields in the surrounding air (as in air band) [46].”\nFigure R3. Field profile (x-component of the electric field) of selected QNMs in the xy-plane. Left panel: the\nfield profiles of QNMs just below the bandgap at ℓ/𝜆≈0.21. Right panel: the field profiles of QNMs just\nabove the bandgap.\nThis figure, which displays typical modes above and below the gap for the present discussion, is not added\nin the revised draft.\n5. Color Range in Aggregates\nIt is quite interesting that both black and white colors can be achieved in the aggregated cases.\nI recommend a more detailed discussion on this point.\nParticle size plays a key role?\nAchieving both black and white from the same material (Si) is hard in conventional systems, as far as I\nknow.\nblack and white are both band response, but sitting at two different sides. One is with complete absorption\nwhile the other without any absorption. I found it is interesting to achieve it in a fixed structure.\nAs mentioned by the authors, a sharp transition occurs from the specular to diffusive ones.\nIntuitively, I expect a gradual change.\nReply:\nThe reviewer correctly identifies that our aggregate structures can appear \"white\" in diffuse light while\nbeing \"dark\" in specular reflection. This behaviour stems from the fact that the aggregate topology acts as\na near-ideal incoherent scatterer, which simultaneously maximizes diffuse scattering and suppresses\ncoherent reflection. The profound structural disorder, with its wide variety of cluster shapes and sizes,\ncreates a dense continuum of resonant states across the visible spectrum. This efficiently scatters light of\nall wavelengths into all directions, producing the bright \"white\" appearance. At the same time, this very\nrandomness completely randomizes the phase of the scattered light. In the specular direction, these phase-\nrandomized waves interfere destructively, effectively cancelling out coherent signal. This is why the surface\nis not truly \"black\" from absorption, but rather \"specularly dark\" due to the suppression of specular\nreflection.\nIf the particles are large enough to scatter efficiently in the visible range, the coexistence of specularly black\nand diffusely white will be robust. This is clearly illustrated in the last three rows of the third column of Fig.\n2c, where increasing the particle size from 95 nm to 130 nm and 170 nm preserves the same dual behaviour.\nAction taken:\nWe have added a new paragraph, a more detailed discussion on this point in section Broadband mirror\neffect at critical packing:\n“Similarly, the aggregate metasurface with minimal correlation exhibits weak specular reflectivity. Here, the\npolydispersity of cluster sizes and shapes, combined with their weakly correlated positions, suppresses coherence,\nleading to dominant diffuse scattering across all wavelengths (Fig. 5).\nThis dual behaviour explains why aggregates can appear bright white in diffuse reflection while remaining dark in the\nspecular direction. The broad continuum of resonant states ensures efficient scattering of all visible wavelengths into\ndiffuse channels, giving rise to whiteness. At the same time, the randomized phases of these scattered waves suppress\nspecular reflection. The effect is robust across particle sizes large enough to scatter efficiently in the visible, as illustrated\nin the last three rows of the third column of Fig. 2c.“\n6. Fig.2 – Color Transition\nThe \"sudden changes\" in color mentioned in the text are not immediately obvious from Fig. 2.\nHowever, Figs. 4 and 5 provide clearer insights. A forward reference in the text noting that this will be\nexplained in more detail later would help guide the reader.\nReply: We agree that a forward reference of Fig. 4 and 5 is important to explain the “sudden changes” in\ncolour when describing Fig. 2.\nAction taken:\nAs suggested, we have added the following sentence in the section Manifold visual appearances around\ncritical packing:\n“A closer inspection reveals pronounced shifts in hue as the system approaches and crosses critical packing, including\nboth a progressive blue-shift and a general whitening of the diffuse light. These spectral trends are not immediately\nobvious in Fig. 2 alone, but will be quantitatively analysed in detail in Figs. 4 and 5.”\n7. Fig. 2c – Black Panels\nWhy are there three black panels in row 2, column 1 of Fig. 2c?\nThis anomaly is not explained in the figure caption or main text.\nReply: This anomaly is now explained in the caption in the last sentence.\n8. Excitation Probability Densities\nThe method used to calculate excitation probability densities is not described.\nA brief explanation or a reference would improve reproducibility.\nReply: We thank the reviewer for pointing this out. The method used to calculate the excitation probability\ndensities is described in the Methods section, and we now explicitly reference this in the main text.\n9. Fig. 5b – Intensity Information\nThe CIE color map in Fig. 5b lacks intensity information. Since the configuration impacts both color and\nbrightness, I suggest adding a separate intensity plot or using a 3D color space to more comprehensively\nshow the effect.\nReply:\nWe thank the reviewer for pointing out the importance of including intensity information.\nAction taken:\nAs suggested, we have added a new figure in the Supplementary Information (Suppl. Section S6), where the\nfive diffuse colours from Fig. 5b are represented in the CIELAB colour space. This representation explicitly\nshows their brightness (L*) in addition to chromaticity. We also provide details of the conversion from\nCIEXYZ to CIELAB space for reproducibility.\nIn the main text, we now explicitly refer to this figure: in the caption of Fig. 5 and in the discussion of Fig.\n5b, where we note that the CIELAB representation both confirms the convergence toward the white point\nand highlights the accompanying increase in lightness (L*).\nFigure R4. Diffuse colours from Fig. 5 represented in the CIE 1976 L*a*b* colour space. This 3D representation\nexplicitly includes brightness (L*) in addition to chromaticity (a*, b*), thereby complementing the 2D\nchromaticity plot shown in Fig. 5b. The conversion from CIE 1931 XYZ coordinates was performed using the\nD65 reference white point.\n10. ref recommendation\nIn solid-state systems, rapid changes in the electronic density of states (DoS) are often associated to\nspectacular phenomena, such as phase transitions for instance.\nConsider the broad readership of Nat Comm, some reference is recommended.\nReply: We fully agree.\nAction taken:\nThree new references have been added in the first sentence of the introduction.\n[1] M. Imada, A. Fujimori, Y. Tokura, \"Metal-insulator transitions\", Rev. Mod. Phys. 70, 1039 (1998).\n[2] E. Yablonovitch, \"Inhibited spontaneous emission in solid-state physics and electronics\", Phys. Rev. Lett.\n58, 2059 (1987).\n[3] E. A. Cornell, C. E. Wieman, \"Nobel Lecture: Bose-Einstein condensation in a dilute gas, the first 70 years\nand some recent experiments\", Rev. Mod. Phys. 74, 875 (2002).\n11. Radiative Cooling Discussion\nThe manuscript compares the brightness of the metasurfaces to that of white paints used for radiative\ncooling.\nHowever, radiative cooling relies not only on high visible/NIR reflectivity but also on high emissivity in the\nmid-infrared (~10 µm).\nI suggest briefly discussing this aspect, and whether Si is suitable for such applications.\nReply:\nRadiative cooling is an important topic nowadays. However, we prefer not to engage in that discussion here,\nsince our metasurfaces absorb predominantly in the visible range.\nAction taken:\nThe part related to radiative cooling has been deleted.\n12. Permittivity of Silicon\nIn Fig. 3, the real part of the permittivity for silicon is set as 20.\nHowever, in Phys. Rev. B 27, 985–1009 (1983), it is reported to be 16.7.\nIs the increase due to quantum confinement effects?\nA reference or short explanation would be helpful.\nReply:\nWe thank the reviewer for pointing out our omission. We have now given a short explanation and a\nreference.\nAction taken:\nWe clarify the reason for choosing 20-0.5𝑖 as the permittivity for Silicon. The sentence is repeated here:\n“For computational efficiency, we assume a fixed silicon permittivity of 20-0.5i, an average value derived\nfrom literature data [45] for the visible spectrum (λ=400 to 750 nm);”\n[45] M. A. Green, M. J. Keevers, \"Optical properties of intrinsic silicon at 300 K\", Prog. Photovolt.: Res. Appl.\n3, 189-192 (1995).\n13. Additional Reference on Structural Coloration\nThe manuscript mentions recent developments in coloration technologies. I recommend including the\nfollowing relevant references to strengthen the background and highlight related progress:\nMao, Peng, et al. \"Manipulating disordered plasmonic systems by external cavity with transition from\nbroadband absorption to reconfigurable reflection.\" Nature Communications 11.1 (2020): 1538.\nCencillo-Abad, Pablo, et al. \"Ultralight plasmonic structural color paint.\" Science Advances 9.10 (2023):\neadf7207.\nReply: We are pleased to cite the first reference, shown as Ref. [33] in the revised manuscript; the second\none is already cited as Ref. [25] in the original manuscript, now Ref. [30].\nReferee #2\nReferee report\nThis work analyses the light scattering properties of disordered metasurfaces made of silicon nanoparticles\n(metaatoms) dispersed on a glass substrate. By controlling the particle size, the surface density and the\npacking fraction (i.e. the range of effective “hard-sphere” correlations), the authors measure the BRDF and\nthe specular reflection in different regimes. They focus on changes in the visual appearance in regimes\nranging from isolated particles to aggregates, and identify the intermediate regime in which metaatoms\nbegin to connect as “critical packing”. They put forward the peculiarities of this regime by different\ntheoretical analyses, and show its potential in the design of thin coatings with substantial impact on visual\nappearance.\nThe work is original, and of potential high interest in the context of metasurfaces for the control of light\nwaves in the visible spectrum. I have some concerns, that are summarized below (major points and minor\npoints). They need to be addressed by the authors before any recommendation can be made. At this stage,\nI cannot recommend publication in Nature Communications.\nOur response\nWe appreciate the general comments of the reviewer highlighting the novelty of the work and its main\nfindings.\nSpecific comments and questions and our step-by-step answer\n1. A huge literature has emerged in the 1980s on metal (gold or silver) semicontinuous films, with studies\non colors and scattering in the far field and in the near field. The authors seem to consider from the\noutset that their systems are very different, and that there is no need to draw on the very many\nprevious studies on metal particle film. I find this a bit hasty. In other words, the context should be\nrevisited. Statements like “Unlike percolation metal films, this transition affects not only the specular\nbut also the diffuse components of the scattered light in a profound way” or “Unlike percolating metal\nfilms, the transition observed in our system influences not only absorption, but also both specular and\ndiffuse optical responses” or “in our dielectric metasurfaces, we do not observe the localized gap\nmodes that play a key role in semicontinuous metal films” seem to me too superficial. Is it true that\ndiffuse scattering of semicontinuous films has not been studied ? Is it clear that their properties are\ninfluenced only by absorption ? Is it clear that the localized gap modes are responsible for most of\ntheir properties ? I would expect a more careful analysis of similarities/differences since comparing\nboth systems might be natural to many readers.\nReply:\nThis comment aligns with the first comment from reviewer #1. We agree that our initial manuscript did not\nprovide a sufficiently detailed comparison with semicontinuous metal films and that this comparison is\nnatural and necessary to properly contextualize our work.\nAction taken:\nIn response, we have revised the introduction to provide a more careful and nuanced analysis of the\nsimilarities and, more importantly, the fundamental differences between the two systems. In the new\nparagraph, we highlight the differences. We repeat here the revised paragraph:\n“The phase transition in our disordered metasurfaces shares similarities with the well-known system of semicontinuous\nmetal films [39,40], yet it also reveals significant differences. Semicontinuous metal films consist of deeply\nsubwavelength (typically <20 nm) nanoislands, which exhibit much larger absorption cross-sections than scattering\ncross-sections and are arranged with limited short-range order. In contrast, our metasurfaces exhibit a different\nhierarchy: scattering dominates over absorption, and we achieve highly precise short-range spatial correlations. These\ndifferences help explain why semicontinuous metal films are predominantly studied for their absorption properties and\nspecular responses [23], while in our metasurfaces, the spectral and angular behaviour of diffuse light plays a central\nrole. Furthermore, the optical properties of semicontinuous metal films are often analysed through a geometric phase\ntransition, where small, disconnected clusters merge into larger, connected clusters—ultimately reaching a percolation\nthreshold that spans the entire system. In our metasurfaces, we also observe a geometric phase transition; however,\nhere, the transition occurs when the metaatoms connect locally, forming small clusters well before a percolation\nthreshold with spanning clusters is reached.”\nWe have also added a new paragraph in the Outlook (see our response to reviewer 1), which highlights\ninteresting possibilities offered by metallic metaatoms significantly larger than those of semicontinuous\nthin films.\n2. The spatial correlations between particles induced by the hard-sphere type potential (exclusion region)\nare mentioned as an important degree of freedom. I would have expected to see a characterization of\nthe degree of correlations other than the packing fraction p. It could be useful to show the structure\nfactors of some key configurations. The structure factor appears in Eq. 1, and is briefly mentioned in\nthe SI, but its role is largely overlooked.\nReply:\nWe agree with the reviewer that explicitly showing the structure factors of hard-sphere type disorder\nprovides a much more direct and intuitive characterization of spatial correlations than the packing fraction\nalone.\nAction taken:\nIn response to this comment, two new figures (shown below as Fig. R5) are added to the Supplementary\nInformation (Suppl. Section S1 and Fig. S1), where we explicitly calculate and display the structure factors\nfor our designs. We have also revised the main text to refer to this new analysis, demonstrating how the\nfeatures of the structure factor directly explain the observed visual phenomena. Specifically, we show that\nthe suppression of 𝑆 at small in-plane wavevectors explains the quenching of diffuse scattering in the\n𝑟\nforward direction, while the emergence and wavelength-dependent dispersion of the main 𝑆 peak account\n𝑟\nfor the formation, colour, and tunability of the halo. This addition emphasizes the central role of the\nstructure factor as a design parameter, complementing the packing fraction description.\nFigure R5. Structure factor characterization. a Static structure factor 𝑆 as a function of the reduced in-plane\n𝑟\nwavevector 𝑞 𝑎 for three different packing fractions 𝑝, where 𝑎 = 1/√𝜌 is the mean interparticle distance.\n||\nb 𝑆 plotted as a function of inverse wavelength for a fixed packing fraction of 𝑝=0.5 and for three different\n𝑟\nparticle densities 𝜌.\nParagraph to be added in main text:\n“This trend is a direct manifestation of the engineered correlations, quantitatively captured by the static structure\nfactor, 𝑆 (see Suppl. Section S1 and Fig. S1a). For weakly correlated patterns (𝑝=0.1), 𝑆 is nearly constant, resulting\n𝑟 𝑟\nin the observed uniform, isotropic scattering. As correlation increases (e.g., 𝑝=0.5), a pronounced dip develops in 𝑆\n𝑟\nat small in-plane wavevectors (𝑞 ≈0), which directly causes the quenching of diffuse light around the specular\n||\ndirection. Concurrently, a strong peak emerges for 𝑞 ≈2𝜋√𝜌, creating the bright, ring-shaped halo of enhanced\n||\nscattering. The dispersive nature of this peak, where its angular position depends on wavelength, is what gives the halo\nits distinct colour. Furthermore, as the particle density 𝜌 increases from 2 to 5 µm⁻² (moving from the left to the middle\ncolumn in Fig. 2c), the primary peak of the structure factor shifts to shorter wavelengths (Suppl. Fig. S1b). This shift\ndirectly explains why the coloured halo appears to expand and change hue, providing another knob for tuning the visual\nappearance.”\n3. The word “percolation” hardly appears in the manuscript. The authors may have a good reason for\nthat, but this should be briefly explained. A (probably naïve) approach could be to see the critical\npacking regime as an intermediate regime between a highly correlated system with isolated particles\nand an aggregate system that could percolate at large density. Whether a point of view based on 2D\npercolation (another widely explored field) is relevant or not should be explained.\nReply:\nWe thank the reviewer for raising this critical point. We agree that the word “percolation” hardly appears\nin the manuscript; this is a deliberate choice. The reason is that percolation and critical packing topologies\nare largely different. Percolation was introduced in the context of semicontinuous metallic films, because\nof the importance of manufacturing continuous films with minimum thickness (characterized by a reduction\nof R ). Then in more fundamental studies, percolation refers to almost continuous metal films supporting\n\nmany random hot spots which are especially important for nonlinear optics.\nOur transitioning is different. We never reach continuous film morphologies in our work. The transition is\nobserved before percolation is reached, between particulate and aggregate topologies. As mentioned by\nthe reviewer in question #1, then if we further for instance increase the density or size, we will transition\nto continuous film morphologies. This could lead to an interesting study.\nHowever, we agree that our original manuscript failed to justify this.\nAction taken:\nWe have added new paragraphs in the Outlook section. Clarifying the difference between the proposed\ncritical packing regime and percolation.\n“It is important to distinguish between critical packing and percolation in semi-continuous metal films, both of which\ninvolve geometrical phase transitions. In disordered metal metasurfaces, percolation refers to the formation of\nspanning clusters [68] that connect the entire sample, resulting in a sharp decrease in sheet resistance. These nearly\ncontinuous systems behave similarly to homogeneous films, dominated by extended networks of hot spots, which play\na key role in absorption and nonlinear optics [39,40].\nIn contrast, the critical-packing phase transition examined here affects local connectivity, where only a small number\nof metaatoms are interconnected, forming small clusters. This occurs much earlier than percolation, long before a\ncontinuous film is established. The morphology remains particulate or aggregate, and the abrupt shift in optical\nbehaviour arises from scattering, mode hybridization, and diffusion rather than percolative transport. Thus, while this\ntransition is reminiscent of percolation due to its geometrical phase transition, critical packing represents a distinct\nregime that broadens the design space for disordered metasurfaces.”\n4. The analysis based on the photonic density of states (PDoS) seem to show that the statistics is\nimportant on revealing the different regimes. I am not sure to understand why the authors do not\nstudy the fluctuations (variance of distributions) of the PDoS both spatially and spectrally. This could\nbe an interesting observable (see for example Krachmalnicoff et al., PRL 105, 183901, 2010) and would\nmake the statistical analysis more substantial.\nReply:\nWe thank the reviewer for this insightful suggestion. Valentina Krachmalnicoff and collaborators study the\nLocal Density of States (LDOS), which is strongly influenced by the underlying structure, particularly by hot\nspots in semicontinuous films. In this context, it is indeed natural to investigate the fluctuations (i.e., the\nvariance of the distributions) of the LDOS, both spatially and spectrally.\nIn contrast, our report focuses on the spectrally resolved global Density of States (DoS), which corresponds\nto the spatially averaged LDOS. This quantity is not only more challenging to compute but also more\nfundamental in the case of distributed emitters, rather than single emitters. Our main emphasis is on the\nspectral dependence of the PDoS during the transition, where we find that the statistical mean of the PDoS\nvaries significantly around the critical packing (Fig. 3). One may also compute the average spectral spacing\nbetween modes (normalized by the average DoS, yielding the “dimensionless conductance” as discussed in\nthe 2018 review on QNMs).\nThat being said, we must admit that we do not have a clear expectation regarding the spectral fluctuations.\nGiven the substantial computational challenges, we have unfortunately not attempted such calculations in\nthe present work.\nAction taken:\nWe have added a comment in the Outlook section to highlight that a deeper statistical analysis of the QNMs\nrepresents an exciting direction for future work. The comment is repeat here for convenience:\n“A deeper statistical characterization of the QNMs, including spectral and spatial fluctuations of the PDoS, would\nrepresent an exciting avenue for future investigation. Such an analysis could reveal additional signatures of the\ntransition and provide complementary insights beyond the mean statistics considered here [69].”\n5. Introduction, line 3: The authors write “the analogue of the electronic DoS is played by the statistics\nof the photon density of states”. Why “statistics” at this stage?\nReply: We fully agree.\nAction taken:\nThe word statistics is removed in the revised draft.\n6. Page 3, paragraph starting with “The 27x2 photographs”: It seems that the incidence angle for\nspecular and diffuse measurements are not consistent with the caption of Fig. 2. This needs to be\nadjusted.\nReply:\nWe thank the reviewer for pointing out the mistake in the caption of Fig. 2.\nAction taken:\nThe incidence angle for specular measurement is corrected in the caption as: “Specular colours recorded\nat 10° incidence for the same metasurfaces are displayed as uniformly coloured rectangles.”\n7. Page 5, paragraph starting with “specular colours tend”: I am not sure to understand the statement «\nThis distinction arises from the fundamentally different physical processes governing diffuse and\nspecular reflection. In diffuse reflection, light experiences the granular nature of the materials and is\nscattered in multiple directions. The material microscopic resonances—such as the properties of\nindividual scatterers—play a central role in this process, while thin-film interference effects are\nrelatively minor. In contrast, specular reflection is driven by collective phenomena at a mesoscale ».\nThe usual understanding is that specular reflection is associated with the average field, while diffuse\nreflection with the field fluctuations. In this setting, I do not see how the processes can be\ndiscriminated between “individual” and “collective”. The authors should clarify their point.\nReply:\nWe agree: the specular reflection is associated to the average field, and the diffuse reflection to field\nfluctuations. It is nevertheless unclear (at least to us) why the resonance of the NPs, that clearly show up in\nthe diffuse light, does not show up in the specular light. Our explanation is that an increase of the effective\npermittivity through a resonance does not directly result in the specular color; our intuitive explanation is\nthat an increase of epsilon may result in an increase or decrease of the transmission through Fabry-Pérot\nresonances built on the coherent field.\nAction taken: We have tried to further clarify the discussion in the revised draft. The revised paragraph is\nrepeated here for convenience:\n“This distinction arises from the fundamentally different physical processes governing the specular and diffuse\ncomponents of the scattered light. Specular reflection is formally associated to the average field in the material, which\npropagates in an effective homogeneous medium and is therefore subject to classical thin-film interference effects.\nThis can lead to striking variations of the specular colour with small changes of metasurface parameters and illumination\nconditions. Diffuse reflection, on the other hand, is associated to the fluctuations of the scattered field around the\naverage. Diffuse colours are mostly driven by the individual scattering elements dressed by the interaction with their\nneighbours, namely, for particulate metasurfaces, by dressed Mie resonances and particle pair correlations. Colour\nvariations with the incident or scattered angles are often milder for the diffuse component than for the specular one.”\n8. The choice of references is sometimes surprising. For example, Page 8, bottom line, Ref. 47 is not about\natomic clouds, Ref. 51 is not about discrete media…. There are many inconsistencies in references all\nalong the manuscript. The authors should systematically check the relevance of their references, and\nuse either the first reference in a field, or a review paper, or a textbook… and not references that may\nseem to have been chosen at random.\nReply:\nWe sincerely thank the reviewer for pointing out the citation error in our manuscript. Indeed, in the\nsubmitted version, Ref. 47 should have been Ref. 48, and Ref. 51 should have been Ref. 54.\nAction taken:\nWe have corrected these mistakes and systematically re-checked all references throughout the manuscript.\nWhere appropriate, we replaced or complemented citations with the most relevant and authoritative\nsources (including key original works, review articles, and textbooks), to ensure accuracy and consistency\nin the references.\nResponse to the referees\nReferee #1\nOverall Comment\nI greatly appreciate the authors' efforts in revising the manuscript, which have significantly enhanced its\nquality.\nI feel all comments have been addressed, and I am pleased to recommend its acceptance.\nOur response\nWe thank again the reviewer for the very positive assessment of our work. We are glad that we could\nanswer all the reviewers’ comments and questions.\nReferee #2\nOverall Comment\nI have reviewed the authors' responses to the questions I had raised, as well as the changes made to the\nmanuscript. The answers provided are accurate, and the changes are appropriate and have clarified points\nthat I felt were essential.\nI also find the work to be highly original, with a broad scope covering the optics of disordered surfaces and\napplications to color rendering. In particular, the concept of critical packing and its essential role in the\noptical properties of complex surfaces is new to my knowledge. I believe it will be of interest to a broad\naudience, with a potentially strong impact. I therefore recommend publication in Nature Communications.\nOur response\nWe thank the reviewer for his/her positive comments. We are glad that we have addressed the reviewer’s\nconcern."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I greatly appreciate the authors' efforts in revising the manuscript, which have significantly enhanced its quality.\nI feel all comments have been addressed, and I am pleased to recommend its acceptance.", "reviewer_2": "I have reviewed the authors' responses to the questions I had raised, as well as the changes made to the manuscript. The answers provided are accurate, and the changes are appropriate and have clarified points that I felt were essential.\nI also find the work to be highly original, with a broad scope covering the optics of disordered surfaces and applications to color rendering. In particular, the concept of critical packing and its essential role in the optical properties of complex surfaces is new to my knowledge. I believe it will be of interest to a broad audience, with a potentially strong impact. I therefore recommend publication in Nature Communications."}, "author_response": "Response to the referees\nReferee #1\nOverall Comment\nI greatly appreciate the authors' efforts in revising the manuscript, which have significantly enhanced its\nquality.\nI feel all comments have been addressed, and I am pleased to recommend its acceptance.\nOur response\nWe thank again the reviewer for the very positive assessment of our work. We are glad that we could\nanswer all the reviewers’ comments and questions.\nReferee #2\nOverall Comment\nI have reviewed the authors' responses to the questions I had raised, as well as the changes made to the\nmanuscript. The answers provided are accurate, and the changes are appropriate and have clarified points\nthat I felt were essential.\nI also find the work to be highly original, with a broad scope covering the optics of disordered surfaces and\napplications to color rendering. In particular, the concept of critical packing and its essential role in the\noptical properties of complex surfaces is new to my knowledge. I believe it will be of interest to a broad\naudience, with a potentially strong impact. I therefore recommend publication in Nature Communications.\nOur response\nWe thank the reviewer for his/her positive comments. We are glad that we have addressed the reviewer’s\nconcern."}]

Decision letter