High-density hyperuniform materials can be transparent

We show that materials made of scatterers distributed on a stealth hyperuniform point pattern can be transparent at densities for which an uncorrelated disordered material would be opaque due to multiple scattering. The conditions for transparency are analyzed using numerical simulations, and an explicit criterion is found based on a perturbative theory. The broad applicability of the concept offers perspectives for various applications in photonics, and more generally in wave physics.


INTRODUCTION
The study of light propagation in scattering media has been a very active field in the past decades, stimulated by fundamental questions in mesoscopic physics [1,2] and by the development of innovative imaging techniques [3]. Recently, a new trend has emerged, with the possibility to control electromagnetic wave propagation in disordered media up to the optical frequency range. On the one hand, wavefront shaping techniques offer the possibility to overcome the distortions induced by a scattering material, even in the multiple scattering regime [4][5][6]. On the other hand, the possibility to engineer the disorder itself by controlling the degree of structural correlation opens new perspectives for the design of materials with specific properties (e.g., absorbers or filters for photonics) [7][8][9][10][11]. These materials combine the advantages of disordered materials, in terms of process scalability and robustness to fabrication errors, with the possibility of developing a real engineering of their scattering and transport properties through the control of the degree of correlation in the disorder. For example, it has been shown that correlations can substantially change basic transport properties, such as the mean free path [12], the density of states [13,14], including the appearance of bandgaps [15][16][17], or the Anderson localization length [18].
A specific class of correlated materials has appeared recently, initially referred to as "superhomogeneous materials" [19], and now called "hyperuniform materials" [20]. These materials are made of discrete scatterers distributed on a hyperuniform point pattern, a correlated pattern with a structure factor Sq vanishing in the neighborhood of jqj 0. The geometrical properties of hyperuniform point patterns have been extensively studied, in particular in terms of packing properties [21][22][23][24][25]. Regarding wave propagation, it has been shown that bandgaps could be observed for electromagnetic waves in two-dimensional (2D) disordered hyperuniform materials [26][27][28][29][30]. Although understanding the origin of the bandgaps is still a matter of study [31,32], these results have stimulated the design and fabrication of threedimensional (3D) hyperuniform structures for wave control at optical frequencies [33,34].
In this Letter, we demonstrate that stealth hyperuniform point patterns, a special class of hyperuniform structures for which Sq 0 in a finite domain around jqj 0, offer the possibility to design disordered materials that can be both dense and transparent in a specific and broad range of frequencies and directions of incidence. The analysis is based on full numerical simulations and theoretical modeling. In the single scattering regime, transparency can be explained in simple terms, as a direct consequence of the vanishing of the structure factor. Interestingly, transparency can survive in the multiple scattering regime under a general condition that we establish using a theoretical analysis that applies to a broad range of materials.

HYPERUNIFORM POINT PATTERNS
A distribution of points in space is hyperuniform when the normalized variance of the number of points within a sphere of radius R vanishes when R tends to infinity [20]. In other words, the fluctuations in the number of points in a volume increase slower with the volume than the average number of points, a feature at the origin of the name "hyperuniform." In Fourier space, this is equivalent to stating that the structure factor defined as [35,36] vanishes in the neighborhood of jqj 0 [23], with N as the number of points and r j as their positions. In this work, we study the scattering of light or other electromagnetic waves by coincide with a hyperuniform point pattern. In the present work, we restrict the numerical study to 2D materials, which is convenient for the sake of computer time, while the theory developed in the last section covers both 2D and 3D geometries. Although the design of 3D materials is a major goal, the interest of 2D architectures should not be underestimated. Artificial 2D structures can produce efficient microwave reflectors or filters [27], and promising 2D materials for photonics and optoelectronics can be fabricated using bottom-up approaches [37]. Even natural materials, such as the cornea, can exhibit photonic properties resulting from an underlying 2D microstructure [38]. The first step consists of generating 2D stealth hyperuniform point patterns. We have adapted the algorithm described in Refs. [21,39,40] and based on the minimization of an interaction potential (see Section 1 of Supplement 1). The 2D medium used in the numerical simulations fills a square region with size L and volume V L 2 , with periodic boundary conditions to avoid finite-size effects. We consider hyperuniform point distributions, such that Sq 0 in a square domain Ω of reciprocal space with size K , centered at the origin. The square shape of Ω has been chosen for convenience. Although it induces an anisotropy in the structure factor, this choice does not affect the generality of the transparency property discussed in this work. Due to the periodic boundary conditions, the structure factor Sq is controlled on a finite number of points located on a square lattice inside Ω and separated by 2π∕L. The number of points M K depends on the size K of the domain Ω and satisfies 2M K 1 K 2 L 2 ∕ 4π 2 , where the factor of 2 results from the symmetry property Sq S−q and the "1" correction accounts for the fact that the value of the structure factor at the exact point q 0 cannot be controlled [Sq 0 N ]. For large K , the system is very constrained, creating stronger correlations in the point pattern. The degree of order in hyperuniform structures is usually measured by the parameter χ M K ∕2N . For χ 0, the system is uncorrelated (fully disordered), while for χ 1, the system is a perfect crystal [21]. From the expression of M K , one immediately gets K 2πL −1 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4N χ 1 p , showing that by increasing K , one increases the degree of order in the point pattern. An example of point distribution generated numerically is shown in Fig. 1(a), together with the corresponding structure factor averaged over 20 configurations in Fig. 1(b). In the square area Ω, the structure factor vanishes [except at the origin, where Sq 0 N ]. For the sake of comparison, an uncorrelated pattern with the same number of points is shown in Fig. 1(c), with the corresponding structure factor in Fig. 1(d). The structure factor of the uncorrelated disordered medium is almost everywhere close to unity.

SCATTERING FROM STEALTH HYPERUNIFORM MATERIALS
A hyperuniform scattering material is built by dressing each point with an electric polarizability αω describing the electrodynamics' response of the individual scatterers at frequency ω (in this study, the scatterers are assumed to be much smaller than the wavelength and are treated in the electric-dipole limit). For 2D waves with an electric field parallel to the invariance axis of the scatterers, the polarizability is αω −4γc 2 ∕ω 0 ω 2 − ω 2 0 iγω 2 ∕ω 0 , where ω 0 is the resonant frequency and γ is the linewidth, this expression being consistent with the optical theorem (energy conservation). The scattering medium is illuminated by a Gaussian beam at normal incidence (direction defined by wavevector k i ) and focused in the middle of the medium, as shown in Fig. 2(a). The beam waist w is chosen large compared to the wavelength to mimic a weakly focused beam similar to a plane wave. The electric field is calculated numerically by solving Maxwell's equations using the coupled dipoles method [41]. A description of the method in a similar geometry can be found in Ref. [42]. It allows us to compute the scattered electric field E sca r at any position inside or outside the medium. In wave scattering by disordered media, a measurable quantity is the average diffuse intensity, proportional to the average power radiated in a given direction in the far field (defined by wavevector k s ). When k 0 r ≫ 1, r being the observation distance and k 0 ω∕c 2π∕λ, with c as the speed of light and λ as the wavelength in vacuum, the scattered field takes the asymptotic form where θ is the observation angle (angle between k i and k s ) indicated in Fig. 2(a). By definition, the far-field average diffuse intensity is where h…i denotes an average over the configurations of disorder (positions of scatterers). We first consider the single scattering regime, defined by b B ω < 1, where b B ω L∕l B ω is the optical thickness in the independent scattering (or Boltzmann) approximation, for which the scattering mean free path is l B ω ρσ s ω −1 and the condition k 0 l B ω ≫ 1 is satisfied. Here, ρ N ∕V is the number density of scatterers, and σ s ω k 3 0 ∕4jαωj 2 is their scattering cross section. Under plane-wave illumination with amplitude E 0 , it is well known that the diffuse intensity can be directly written in terms of the average structure factor. The full expression takes the following form (see Section 2 of Supplement 1):

Research Article
where q k s − k i is the scattered wavevector. On the right-hand side, Θq is the Fourier transform of a window function equal to unity inside the volume V occupied by the medium and vanishing outside. This function describes the diffraction pattern of the finite size volume V , and in practice, it is sharply peaked around q 0. In particular, for q 0 corresponding to forward scattering, the two terms on the right-hand side in Eq. (4) exactly compensate each other.
In the single scattering regime, the vanishing of Sq in a finite domain Ω suppresses scattering for a given frequency and/or angular range, a property at the origin of the denomination of stealth hyperuniform materials [39]. To illustrate this interesting behavior, we choose the parameters such that b B 0.5, and we change the incident wavelength, or equivalently, k 0 , in order to tune the range of scattered wavevector q k s − k i k 0ks − k 0ki that is involved in the scattering process (here,ˆdenotes a unit vector). The first situation is chosen with k 0 K ∕8, so that the circle described by the scattered wavevector q in Fig. 2(b) lies entirely within the domain Ω. In this case, we expect a substantial reduction of scattering for a hyperuniform material [with Sq ≃ 0 in Ω], compared to a fully disordered material with the same density of scatterers. This is observed in Fig. 2(c), where we plot the average diffuse intensity versus the observation angle θ calculated from the full numerical simulation without approximation (the averaging is performed over 20 configurations). Also note that this regime is observed as long as the circle described by the scattered wavevector q lies entirely within the domain Ω, independently of the exact shape of Ω. For the uncorrelated disordered medium (red dashed line), a diffuse intensity pattern is observed, while for the hyperuniform structure made with the same scatterers and the same density, the diffuse intensity (blue line) is negligible [reduced by a factor of about 800, so that it is hardly visible in Fig. 2(c)].
For smaller wavelengths (larger k 0 ), another behavior can be observed. Choosing k 0 K ∕3.5 as an illustrative example, the circle described by the scattered wavevector q in Fig. 2(b) is now only partially included inside the domain Ω. Therefore, a reduced level of scattering, due to a vanishing structure factor, is expected only for the observation angles θ ∈ −θ lim ; θ lim such that q falls within Ω. Theoretically, we find θ lim arccos1 − K ∕2k 0 138°. The result of the full numerical simulation of the diffuse intensity pattern is shown in Fig. 2(d). The change between the pattern produced by an uncorrelated disordered medium (red dashed line) and the hyperuniform material (blue line) is clearly visible. For the latter, scattering is suppressed for a wide angular range, with the value θ lim 134°in good agreement with the estimate from single scattering theory (for b B 0.5, the higher-order scattering is not fully negligible). Finally, let us note that in terms of frequency bandwidth for a predefined angular range −θ lim ; θ lim , scattering would be suppressed for all wavelengths satisfying λ > λ lim , with λ lim 4π∕K 1 − cos θ lim .
The examples above have illustrated the interest of hyperuniform point patterns in the design of materials with a tunable level of scattering. An interesting question is to see under which conditions these properties could survive in the multiple scattering regime. The transport of the average intensity in this case is well described by the radiative transfer equation [43], which is derived from the Bethe-Salpeter equation in the limit k 0 l B ≫ 1 [2]. In this framework, the angular dependence of the diffuse intensity is driven by the phase function pu; u 0 , where u 0 and u are unit vectors denoting the incoming and outgoing directions. For subwavelength scatterers, the phase function is connected to the average structure factor through the relation (see Section 3 of Supplement 1) where the term involving Θ appears as in Eq. (4), but with a different weighting factor. This difference results from the connection between the phase function and the vertex of the Bethe-Salpeter equation that contains irreducible terms only. Note that for an uncorrelated disordered medium, the phase function is constant [i.e., pu; u 0 1]. In the case of a hyperuniform structure, under the condition k 0 < K ∕4, the scattered wavevector k 0 u − u 0 falls inside the domain Ω [see the construction in Fig. 2(b)], and the phase function vanishes for all directions, except in the forward direction u u 0 . Indeed, the second term involving Θ in Eq. (5) does not exactly compensate the value of the structure factor at the origin. This means that scattering occurs, but only in the forward direction. For an infinite medium, this would lead to an effective scattering mean free path l tending to infinity. In practice, for a medium with size L, the width of Θk is on the order of 2π∕L, keeping l finite but permitting it in principle to reach arbitrary large values. This is particularly interesting, since getting l ≫ L leads to transparency. To check this idea, we  Fig. 3, where we clearly observe that the average diffuse intensity for the hyperuniform medium (blue line) is extremely small compared to that obtained for a fully disordered medium (red dashed line), except in the forward direction, as expected.
Our analysis demonstrates that, in conditions (in terms of type and density of scatterers) under which a fully disordered material would induce strong multiple scattering, a hyperuniform material can chiefly generate forward scattering, leading to an effective scattering mean free path l ≫ L, which is the condition for transparency.

EXPLICIT CRITERION FOR TRANSPARENCY
It is possible to establish a criterion for the existence of transparency in a dense (stealth) hyperuniform disordered material, based on a perturbative theory valid for both 2D and 3D materials. To proceed, we use a diagrammatic expansion of the phase function in the form [44] (6) where the circles denote scatterers, the horizontal solid lines represent free-space Green functions G 0 , the vertical solid lines stand for identical scatterers, and the dashed lines connect different but correlated scatterers. The first two terms in Eq. (6) correspond to Eq. (5), and according to the analysis above, they compensate for a stealth hyperuniform medium, provided that k 0 < K ∕4. The value of the effective scattering mean free path is obtained by an angular integration (over direction u) of the phase function [44]. For a "standard" short-range ordered material, the first-order correction is obtained by integrating the first two terms, leading to an expression equivalent to that introduced initially in Ref. [12] to describe multiple scattering in interacting (but non-hyperuniform) suspensions. In the specific case of stealth hyperuniform disorder, the first two terms do not contribute, and the effective scattering mean free path is obtained by an angular integration of the next four diagrams. For subwavelength scatterers with polarizability αω, this leads to RαωG 0 rhrF k 0 rr d −1 dr; where h is the pair correlation function [20] and d is the dimension of space. Note that for large scatterers with sizes comparable with the wavelength (such as Mie spheres), the dependence of the scattering amplitude of an individual scatterer on the scattered wavevector q would change the analytical expressions. We can expect, in this case, the results to change qualitatively. For scalar waves in 2D, F x J 0 x and G 0 r i∕4H 1 0 k 0 r, with J 0 (H 1 0 ) as the zero-order Bessel (Hankel) function of the first kind. For scalar waves in 3D, F x 2 sincx and G 0 r expik 0 r∕4πr. An estimate of the right-hand side can be obtained by computing the self-energy Σ. To first order in k 0 l B −1 , one has [2] (8) The integral term in Eq. (9) is on the order of k 0 l B −1 , which can be used to estimate the right-hand side in Eq. (7), leading to l ∼ l B × k 0 l B . Transparency is observed provided that l ≫ L, which leads to the criterion b B ≪ k 0 l B : (10) Under this condition, a stealth hyperuniform material becomes transparent for frequencies satisfying k 0 < K ∕4. In particular, this means that even for b B ≫ 1 (corresponding to the multiple scattering regime for a fully disordered material), transparency can be reached, provided that Eq. (10) is satisfied. Therefore, our analysis establishes a criterion for transparency of a hyperuniform material even at high density, i.e., far beyond the single scattering regime.

CONCLUSION
In summary, we have shown that materials made of nonabsorbing subwavelength scatterers distributed on a stealth hyperuniform point pattern can be transparent, even in conditions under which, for the same density of scatterers, an uncorrelated disordered material would be opaque due to multiple scattering. This occurs under very general conditions that we have established theoretically for 2D and 3D materials. In this first study, the numerical examples have been restricted to 2D materials (which by themselves are of practical interest) for the sake of computer time, but there is no fundamental limitation for the extension of the numerical simulations to 3D. The transparency property of dense hyperuniform media opens new perspectives in the engineering of disordered materials combining specific photonic properties with a high robustness to fabrication errors. More generally, the design of dense and transparent hyperuniform materials is in principle achievable for any kind of wave, and the applications of the concept cover many fields of wave physics.

METHOD TO GENERATE HYPERUNIFORM POINT PATTERNS
This section is dedicated to the presentation of the method used to generate 2D hyperuniform point patterns. The algorithm is adapted from Refs. [1][2][3]. Numerically we can only manipulate a finite number of points and to mimic an infinite medium, we divide the space into identical square cells of size L, each containing N points, as shown in Fig. S1. The system is supposed to be L-periodic and we only consider points in one unit cell. The structure factor for the full system is thus given by exp[iq · (n x Lu x + n y Lu y + r j )] 2 (S1) where u x and u y are unit vectors along the x and y axes respectively. This leads to the following expression where q = (q x , q y ) and is the structure factor of the unit cell. This expression implies that S tot (q) = 0 for all q such that q x and q y are not multiple of 2π/L. Thus to get S tot (q) = 0 for all q ∈ Ω, we only have to require a vanishing of the last term for any q = (2πn/L, 2πm/L) ∈ Ω, that is where Ω ′ is the ensemble corresponding to all q ∈ Ω such that q = (2πn/L, 2πm/L), as shown in Fig. S2. B(q) is always positive or zero. As a result, to make the structure factor vanishing in Ω ′ , B(q) must reach its global minimum (i.e. zero) for all q ∈ Ω ′ . The simplest way to implement this is to define a potential φ by and to minimize this potential. In the case of a square domain Ω with size K, this leads to the following expression for φ: Ensemble Ω ′ of points where the structure factor of the elementary cell should vanish to get a hyperuniform point pattern.

RELATIONSHIP BETWEEN THE AVERAGE DIFFUSE INTENSITY AND THE STRUCTURE FACTOR IN THE SINGLE-SCATTERING REGIME
In this section, we derive the relationship between the average diffuse intensity and the structure factor in the single scattering regime. We consider a 3D geometry, but the derivation is essentially the same in 2D (and can be easily adapted). The average diffuse intensity is given by the difference between the average intensity and the intensity of the average field (ballistic component) through the relation where q = k s − k i . The first term corresponds to the average structure factor and is given by The second term is more complex. To have a more explicit expression, we have to consider the probability P(r) = V −1 of having a scatterer at position r which gives We now define Θ as Thus Θ is a window function, the Fourier transform of which gives the second term in the average diffuse intensity. Finally, we obtain which corresponds to Eq. (4) of the main text.

RELATIONSHIP BETWEEN THE PHASE FUNCTION AND THE STRUCTURE FACTOR
In this section, we derive the relationship between the phase function that enters the radiative transfer equation (RTE) [4] and the structure factor in a 3D geometry (the derivation can be easily adapted to 2D). We recall that the phase function corresponds to the radiation pattern of a single scattering event in a macroscopic sense (after averaging over disorder), and the RTE is a transport equation for the specific intensity (averaged over disorder as well). The RTE can be derived from first principles (Maxwell equations) [5,6]. More precisely, the starting point of the derivation is the Bethe-Salpeter equation [7][8][9] which, in statistically homogeneous medium, reads as where . . . denotes the ensemble average over the configurations of disorder, E(r, ω)E * (ρ, ω) being the spatial correlation of the field and G(r − r ′ , ω) the average Green function. Γ is called the intensity vertex operator. It is a very complex object that contains all multiple scattering events that cannot be factorized when averaging over disorder. In particular, Γ contains the information on the spatial correlations of the disrodered medium. From Eq. (S12), and considering a dilute system with k 0 ℓ B ≫ 1, one can derive the RTE, and in particular the expression of the phase function that reads p(u, u ′ ) = A Γ(k r u, k r u ′ , k r u, k r u ′ , ω) (S13) where the constant factor A is such that p(u, u ′ )du = 4π.
is the effective wave vector (i.e. the wave vector appearing in the wave equation for the average field). In a dilute system, we have k r ∼ k 0 . Γ is the reduced vertex operator in Fourier space, given by Γ(k, k ′ , κ, κ ′ , ω) = (2π) 3 δ(k − k ′ − κ + κ ′ ) × Γ(k, k ′ , κ, κ ′ , ω). (S14) In a dilute system, it is convenient to use a Taylor expansion of Γ with respect to the small parameter [k 0 ℓ B ] −1 . In the case of a correlated medium (such as a hyperuniform material), the first two terms are needed, and the expansion is where the first term corresponds to a scattering process for the field and its complex conjugate involving the same scatterer (solid line between both circles representing the scatterers, the upper row being related to the field and the lower row to its conjugate). The second term describes a scattering process involving two different but correlated scatterers (dashed line). Analytically, for subwavelength satterers with polarizability α(ω), we obtain Γ(r ′ , r ′′ , ρ ′ , ρ ′′ , ω) = N V k 4 0 |α(ω)| 2 δ(r ′′ − r ′ )δ(ρ ′′ − ρ ′ ) where h is the pair correlation function [10]. This leads to the following expression for the phase function: To get the relationship with the structure factor, we now explicitely write S in terms of the pair correlation function. We have × P(r j )P(r k )[1 + h(r j , r k )]d 3 r j d 3 r k × [1 + h(r 1 − r 2 )]d 3 r 1 d 3 r 2 From Eqs. (17) and (18), we finally obtain Eq. (5) of the main text: The appearence of the correction term in the right-hand side results from the definition of the phase function, that does not involve the same diagrams as the structure factor. In the phase function, or equivalently in the intensity vertex, we only take into account irreductible diagrams (i.e. diagrams that can not be factorized after the average process). In the structure factor, we have an additionnal diagram corresponding to two different scatterers that are not correlated, so that the expression of the averaged structure factor reads S = + + . (S20)