Enlightening complexity: route to strong localization

By employing Random Matrix Theory (RMT) and first-principle calculations, we investigated the behavior of Anderson localization in 1D, 2D and 3D systems characterized by a varying disorder. In particular, we considered random binary layer sequences in 1D and structurally disordered photonic crystals in two and three dimensions. We demonstrated the existence of a unique optimal degree of disorder that yields the strongest localization possible. In this regime, localized modes are constituted by defect states, which can show subwavelength confinement properties. These results suggest that disorder offers a new avenue for subwavelength light localization in purely dielectric media.


Introduction
One of the most fascinating phenomena occurring in disordered systems is Anderson localization [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19]. This dynamics relies on the restoration of interference effects due to coherent scattering, which alters the eigenmodes of the system from extended to localized and traps the electromagnetic radiation into long-living resonant states. The transition to Anderson localization is normally manifested as the breaking of light diffusion, with the emergence of localized states whose properties are described by a characteristic length -namely the localization length l c -that measures the spatial degree of localization achieved through Anderson modes [20]. Diffusion and localization can share an intermediate stage known as weak localization, which acts as a precursor for Anderson localization to settle in [20,21]. In the research panorama initiated by Anderson, an open problem is represented by the regime of strong localization of light, where interference effects are expected to be more pronounced [20]. In this situation, a full vectorial solution of Maxwell's equations [3] can be applied to rigorously study the problem.
Among the various open issues, a fundamental question concerns the role of disorder in sustaining strongly localized states, and how to effectively vary the randomness in order to control localized modes and achieve strong localization effects. This problem is essential not only to capture new insights on the dynamics of localization, but also to identify new avenues for the application of disordered systems to the manipulation of light properties. In this Article, we combine Random Matrix Theory (RMT) with massively parallel Finite-Difference Time-Domain (FDTD) simulations, and study the behavior of light localization versus disorder in 1D, 2D and 3D systems. More specifically, we begin with one dimensional media and employ RMT [22] to demonstrate the existence of a unique minimum in l c , as well as the existence of an optimal degree of disorder for localizing light. We then considered 2D and 3D geometries, and evaluated through parallel FDTD simulations the behavior of l c in the spectral region where the highest density of localized states is expected to be observed. In particular, we considered different realizations of disordered photonic crystals and analyzed the properties of localized states near the band edge of the photonic bandgap, which is expected to support a high density of localized states [2]. As in the 1D case, we found that the variation of l c is characterized by a bell-shaped behavior showing a unique minimum associated to the strongest localization possible. For average disorder values different from this optimal condition, strong localization of light is completely lost. This contributes, in part, to explain the elusive character of this phenomenon. Quite remarkably, our analysis reveals that in the condition of maximal localization, electromagnetic modes are constituted by defect states, whose localization length is smaller than one cell of the periodic crystal and has a subwavelength spatial extension, which in both 2D and 3D cases is approximatively 60% smaller than a single internal wavelength in the high index refractive material. This demonstrates that a proper engineering of disorder can be effectively employed for subwavelength control of light within purely dielectrics.

1D systems: a unique minimum in the localization length
We begin by considering the (X, Z) propagation of light into a disordered one-dimensional system characterized by N different layers, each possessing a refractive index n i and thickness h i , which are randomly varying on the X axis. Maxwell's equations are solved by the transfer matrix approach, by looking for solutions of the type E y (x, z) = E 0 2 U (x)e iωt−iβ z + c.c. and H z (x, z) = E 0 2icµ 0 V (x)e iωt−iβ z + c.c., where β is the propagating constant of the mode with frequency ω, x = X/k 0 , z = Z/k 0 dimensionless coordinates (k 0 = 2π/λ ) and U , V dimensionless fields with E 0 is an arbitrary electric field constant. Input U 0 = U (0), V 0 = V (0) and output U l = U (L), V l = V (L) (L is the sample length along x) are related by the transfer matrix [23] with γ j = n 2 j − n 2 β and n β = β /k 0 . The matrix M j belongs to the Lie group SL(2,R) and provides a non orthogonal rotation to the input fields U 0 and V 0 . According to Furstenberg theorem [22], when the elements of M j are random variables, the following limit holds: lim N→∞ log Tr(M) = δ > 0 [ ... denotes an ensemble average] and all modes get exponentially localized as U , V ∝ e −δ |x| , where δ is the positive Lyapunov exponent of the random matrix M [22]. The localization length l c can be directly computed from the Lyapunov exponent and reads l c = δ −1 . We focused our analysis in the important physical situation where the disordered material is composed by a random binary sequence of two different media, each characterized by a diverse index n i and/or thickness h i , with i ∈ [1,2]. These define two different transfer matrices M 1,2 , which we labeled A and B for notation simplicity. Random deviates are chosen from a binary probability density characterized by two discrete values: p and q = 1 − p, which represent the probability to extract n 1 (h 1 ) or n 2 (h 2 ) values, respectively. A binary superlattice is not only physically interesting per se [24,25], but it is also fundamental to understand the behavior of the localization length l c in higher dimensions. In 2D and 3D systems, in fact, the majority of the experiments is conducted in samples showing the random aggregation of two different materials such as, e.g., scatterers hosted in a liquid matrix (a photonic colloid) or disordered solids where the second medium is composed by air holes [4,6,7,10,26]. In all these cases, a 1D representation of the dynamics along a specific direction of propagation, or parallel to a light ray that is randomly scattering inside the sample, would be equivalent to that of a random sequence of two different layers. We analyzed the dynamics of l c versus disorder by applying a fully analytic theory based on the Microcanonical method. The method is well documented in [22,27] and in this paper we show only the main results achieved with this approach. We focused on the following general case: with M 0 a constant matrix and ∆M a random term. In order to evaluate the behavior of l c , we calculated the generalized Lyapunov exponent L (1) = lim N→∞ log Tr(M) =δ that, apart from negligible fluctuations, coincides with δ [22,27]. The ensemble average can be expressed as follows: with Tr(M 0 ) a constant value and Tr(∆M) expressed by the following Cauchy integral [22]: Equation (4) can be evaluated with the saddle point method and, after straightforward algebra, we obtain: with h a complex value that satisfies the saddle point equation (5), which can be employed to study the behavior of the localization length versus disorder. The latter can be measured by the probability p of extracting the matrix A from the binary ensemble. Values of p = 0 or p = 1, in particular, describe an ordered system composed by a constant configuration of A or B matrices, while p ∈ (0, 1) yields different disordered sequences of A and B. By applying a derivative on p, we get ∂ δ /∂ p = 2 tan −1 (1 − 2p) − lnh, which shows a single zero at p * : in the material. We verified this theoretical prediction by performing a series of numerical simulations. Without any loss of generality, we set the arbitrary term M 0 to M 0 = 0. We then considered two matrices A and B with γ a = π/4, γ b = γ a + 0.2 and h a = h b = 1 (since h j appears always in the factor γ j h j , a random variation of γ j is enough to have a binary random sequence) and compared the localization length l c with Eqs. (5)- (6). An excellent level of agreement is found between theory and numerical results [ Fig. 1(b)], with the localization length showing a single minimum as predicted by Eq. (6).

2D and 3D media: defects states and subwavelength effects
We studied the behavior of higher dimensional systems by resorting to FDTD simulations, and performed a massively parallel simulation campaign (about 10M of single-cpu hours) on a variety of experimentally feasible 2D and 3D structures. As a measure of l c , we employed the Inverse Participation Ratio (IPR): withl c the IPR of a single realization of disorder, V the system volume and E = 1 2 (E · D + H · B) the electromagnetic energy density. The operator ... represents an average over realizations. Random samples are constructed by disordering a Photonic Crystal (PhC) of air holes on a Silicon (Si) substrate. In particular, we considered positional disorder of the type x = x 0 + ∆, with ∆ being a two/three element random vector whose values ∆ i are given by a normal distribution with ∆ i = 0 and standard deviation ∆ 2 i = χ. In order to prevent folding effects, we constrained the global shift of each hole to be lower than the lattice constant of the PhC, and chose the standard deviation χ of disorder accordingly. In our FDTD simulations, in order to get accurate results, we adopted a spatial discretization with at least 20 internal points per wavelengths in the high-index (Si) material. Two dimensional samples are constructed from a triangular lattice, 66µm ×66µm wide along x and y, characterized by air holes of radius r = 0.3µm. Figure 2(a) shows the refractive index distribution of a single realization for χ = 0.15. In our simulation campaign, we considered on average 20 disorder realizations for each value of χ, which we have verified is sufficiently high to yield a good convergence of the averaged observables. In order to maximize the probability of observing localization effects, we employed a wideband gaussian source with FWHM=250nm centered at λ 0 = 4.8µm [ Fig. 2(b) dashed line], whose wavelength is close to one of the gap edges of the photonic crystal, as shown by the Local Density Of States (LDOS) of the system for χ = 0 [ Fig. 2(b) solid line]. In each simulation, we launched the gaussian source and then monitored the evolution ofl c (t) over time t [Fig. 2(c)]. After an initial transient of ≈ 4 · 10 5 timesteps, the source propagated over the system and the localization lengthl c reached a steady state, characterized by the electromagnetic modes excited in the system after the interaction with the input wave [ Fig. 2(c)]. The behavior of l c , after averaging, is qualitatively identical to the one dimensional case, with a single minimum attained in the dynamics [ Fig. 2(d)]. In order to investigate the properties of light in the regime of strong localization, we begin by calculating the LDOS for the strongest localization, which is attained for χ = 0.15 [ Fig. 2(e)]. As seen in Fig. 2, electromagnetic states are composed by a few eigenfrequencies located inside the badgap of the photonic crystal, and they effectively behave as single-defect states. As a result, their spatial localization is nor originated by Anderson interference effects, neither by multiple impurity scattering, but conversely by the localization properties of the electromagnetic bandgap. The consequences of this dynamics are quite interesting, as shown in Fig. 2(f), which displays the spatial distribution of the electromagnetic energy E for one realization at As seen, the electromagnetic radiation is able to reach a subwavelength confinement.
The FWHM of the energy spot, in fact, is calculated to be FW HM = 0.58µm, which is 60% smaller than one internal wavelength λ = λ 0 /n = 1.38µm, where n is the refractive index of Si. It is worthwhile remarking that such a subwavelength confinement can be only observed in the narrow range of χ that yields a strong localization regime, i.e., for χ ∈ [0.13, 0.18]. As seen in Fig. 2(d), in fact, a variation of a few percent of average disorder from its optimal condition yields an increase of the localization length, with the vanishing of subwavelength confinement. We completed our theoretical analysis by the study of 3D systems. FDTD simulations have been carried out by disordering a 18µm × 18µm × 18µm fcc lattice of air holes (with radius r = 0.35µm) in silicon, and by employing a Gaussian input source with λ 0 = 1.5µm and FWHM = 100nm. As in the 2D case, we have centered the source at one edge of the photonic crystal gap in order to maximize the probability of observing localized states [ Fig. 3(a) dashed line]. Figure 3(b)-(d) summarizes the results arising from numerical simulations. The behavior of the localization length l c is qualitatively identical to the 2D case [ Fig. 2(d)], with a unique minimum observed as χ is varied. Quantitatively, the strongest localization is attained at χ = 0.05 [ Fig. 3(b)], which corresponds to a small variation of disorder of only 5%. In this regime, as for the two-dimensional case, strongly-localized states are characterized by defectstates originated in the electromagnetic bandgap, as observed by the LDOS reported in Fig.  3(c). In order to further investigate the degree of spatial confinement achieved in the strong localization regime, we considered a single realization for χ = 0.05 and evaluated the FWHM of the hot energy spot [ Fig. 4(a)-(b)], obtaining a value FW HM = 0.19µm, which is ≈ 56% smaller than one internal wavelength λ 0 /n = 0.43µm in silicon. The three dimensional system therefore shows the same qualitatively scenario of the two dimensional case, with comparable subwavelength confinement properties -i.e., 60% in 2D vs 56% in 3D-although it offers a much narrower disorder window for achieving strong localization effects, as observed by comparing Fig. 3(b) with Fig. 2(d).

Discussion and conclusions
We investigated the role of disorder in sustaining strong localization of light in 1D, 2D and three-dimensional systems. In one-dimensional media, we applied a Random matrix analysis and demonstrated the existence of a unique minimum that yields the strongest localization possible. A set of numerical simulation has then been employed to verify this result. In 2D and 3D media, conversely, we tackled the problem by an extensive campaign of FDTD simula-tions, which kept into account the vectorial nature of Maxwell's equations. More specifically, we considered a series of samples characterized by structurally disordered photonic crystals. The behavior of the localization length l c in these systems, as obtained by FDTD simulations, appears qualitatively identical to the one dimensional case, with the manifestation of a single minimum that sustains the strongest localization of light possible. From a more quantitatively perspective, in both 2D and 3D samples, FDTD simulations predicts the possibility of observing strong localization effects for a relatively small average disorder : χ = 15% in two dimensions and χ = 5% in the 3D case. In addition, the numerical analysis shows that a small variation of disorder from the optimal condition leads to a complete vanishing of the strong character of localized modes, with the localization length sharply increasing to the initial value pertaining to the unperturbed crystal. The average disorder range for the observation of strong localization effects is quite narrow: approximatively 5% in 2D and only 1% in three dimensions. These results contribute, to a certain extent, to also explain the elusive nature of this phenomenon. An interpretation for the qualitatively similar behavior of l c versus disorder in 1D, 2D and 3D systems can be drawn by generalizing the one dimensional RMT analysis to higher manifolds. In 2D and 3D cases, in fact, we can fold the photons dynamics over the scattering trajectory of light rays. A single ray, in the most common configuration where only two materials are employed (i.e., a suspension of particles in a host matrix or air holes on a disordered substrate [26]), undergoes a chaotic motion inside the sample and practically experiences a random succession of two layers. The presence of a higher dimensional manifold allows a richer dynamics of the evolution of light paths if compared to the one dimensional case, which might lead to a different scenario for light localization (e.g., mobility edge, marginal localization, coexistence of localized and extended states, ...). However, among all the various paths available, we can intuitively expect that the variation of disorder can lead to the formation of a minimum in l c when, on average, light rays scatter along trajectories whose disorder is optimal for their one dimensional evolution. Regardless of the presence of a higher manifold for the dynamics of light, we can therefore qualitatively expect the appearance of a single minimum, if the average disorder is able to optimally localize the dynamics of light rays. As a final step of our theoretical analysis, we investigated the nature of localized modes in the strongest localization regime. Calculation of the local density of states revealed that in this condition light localization is sustained by defect states originated by a few eigenfrequencies inside the electromagnetic band gap. This results in a very strong spatial confinement for light, due to the excitation of such disordered-induced impurity modes. From a physical perspective, such a strong confinement can be explained by considering that a defect mode results from the interference of a large number of plane waves that, being localized in the spectral bandgap, cannot couple to any extended state of the structure and are forced to interfere at every wavevector k. The ideal case of an infinite sum of plane waves, which coherently interfere at every angle, yields the smallest spot achievable, i.e., a dirac-delta δ function. Conversely, the opposite case of an incoherent summation of plane waves originates the well-known Rayleigh distribution (see, e.g., [28]) that does not show any hot-spot with subwavelength confinement. Our analysis shows that suitable values the disorder support a coherent combination of plane waves that lives in between these two extrema -i.e., a full coherent combination and a totally incoherent interference-, which results into a subwavelength spatial confinement for light. A proper engineering of disorder can therefore be employed as an avenue for subwavelength light localization, where extremely narrow hot energy spot can be generated to highly enhance absorption and/or emission properties of purely dielectric materials. This is expected to open new applications in the field of high power lasers, sensing and energy harvesting devices.