Bandgap Simulations in Randomized 3D Photonic Crystal Supercells

1. Introduction

Photonic crystals (PCs), or photonic bandgap (PBG) crystals, are periodic dielectric materials that are designed to control how electromagnetic waves propagate [1]. The crystal is designed to prevent or inhibit the propagation of light at certain frequencies, either through photonic bandgaps or by significantly modifying the density of electromagnetic states [2,3,4]. PCs are an area of active research, with strong emphases on both two-dimensional and three-dimensional PCs for their many practical uses, such as communication [5], computing [6,7] and materials science [8,9]. Photonics, more widely, is the control of photons in free space or matter, in analogy to the control of electron-charge seen in electronics [10]. The most widespread use of photonics at present is for telecommunications through the optical fiber waveguide [11,12], which allows for data transfer at much higher bandwidth, higher efficiency, and lower cost than equivalent electronic communications methods. Optical computing, also called photonic computing, is another future application of photonics that proposes using light rather than electrons as the medium for computers for the similar reasons [13]. PCs constructed with bandgaps to control the propagation of certain frequencies would be the foundation to such photonic computers [14,15].

Recent research has sparked interest into the effects of structural disorder within Photonic crystals. Studies focused on biological and biomimetic materials [16,17,18], reveal that such disorder can exert both positive and negative influence on the efficacy of light propagation in these structures. Whether light localizes or delocalizes hinges on the specific disorder introduced. The extrinsic introduction of disorder into photonics, and the engineering of topologies associated with disordered structures [19], has revealed new ways for manipulating the light spectrum, transport, and wavefront properties [20,21,22]. These variations can also lead to unique optical properties and functionalities not present in perfectly ordered PCs. For instance, the disorder can facilitate the phenomenon of Anderson localization of light, where light waves are trapped and localized due to scattering from random defects in the crystal structure, or irregularities in the crystal lattice [23,24]. Notably, the utilization of materials lying along the spectrum between the regimes of perfect order and disorder shows promise in yielding innovative optical phenomena [25,26]. The simulation of disorder in photonic crystals aids in comprehending its impact on light propagation and in the design of materials with tailored optical properties. Using self-stabilized colloid deposition, scalable methods to fabricate tailored disorder have been presented [27], this allowed control over lateral and vertical statistics of the patterns and accurate prediction of particle patterns is investigated using structure factor and optical measurements.

In this paper, the effect of randomization disorder on photonic crystals was simulated using the Massachusetts Institute of Technology (MIT) Electromagnetic Equation Propagation (MEEP) and Photonic Bands (MPB) open source freewares aimed at simulating Maxwell’s equations, where the evolution of electromagnetic waves over time is modelled in materials [28,29]. The output from MEEP and MPB was confirmed to match expected results using photonic crystal unit cells (fcc lattices of packed spheres) with known band structures from literature. Supercell structures were created, with significant variations in the radius size added as an element of disorder.

Three-dimensional inverse opal structures (see Figure 1) of a face-centered cubic (fcc) lattice with air spheres in a dielectric, are simulated both in the ordered state and with a randomization of the sphere radii. A fcc unit cell known to have a bandgap of 3.5% (of normalized frequency range) in the ordered configuration, with sphere radius 0.34a where a is the lattice constant, was used as the reference. The sphere radii are then randomized according to a Gaussian distribution, with a mean radius μ /a = 0.34. Multiple standard deviations (σ) were used in simulations, ranging from σ/a = 0.0025 (~0.7% of 0.34a) to σ/a = 0.020 (~5.9%). In order, to generate a statistically significant body of outputs, ~500 independent simulations runs were executed for each value of σ studied.

A reduced prevalence and width of the bandgap in 3D systems with a small variance disorder was shown and is expected [30] and shows that the introduction of random disorder has a crucial influence on the maintenance and degradation of PBGs. However, some PBGs were shown to be larger than the ordered configuration. This ability to use engineered disorder to create tailored optical bandgap properties, following the outputs of simulation, is of particular interest and warranting of further investigation and future optimization.

2. Materials and Methods

The MEEP software package for Python (v1.29.0), installed on Linux Ubuntu 22.04, was used to generate all the reported simulations. NumPy was used extensively throughout, as well as Matplotlib for general data plots and VisPy for 3D plots. An AMD Ryzen 7 processor with 32GB RAM was used for development, while simulations were performed on an Intel i7 processor and with the use of Supercomputing Wales. Both systems produced numerical results identical to within the pertinent number of significant figures.

2.1. Simulation Software

MEEP and MPB are a bundled set of free open-source software for simulating electromagnetics in materials. MEEP is specialized for time domain simulations using the finite-difference time domain (FDTD) method, where it simulates how electromagnetic waves evolve over time in a material [28,31]. FDTD is a numerical analysis method for solving computational electrodynamics; Maxwell’s equations are replaced by a set of finite difference equations, such as Δf = f(x + 1) − f(x), and then these are solved for given boundary conditions and appropriate field points [32].

MEEP and MPB do not use standard SI units [31,33]. Instead, they use a set of scale-invariant dimensionless units where c = 1. The frequency f is given in normalized units of c/a (and angular frequency ω in 2πc/a). MPB outputs the frequencies from 0 to c/a in the Brillouin zone. Hence, supposing if MPB showed a bandgap at f = 0.5 c/a, and the lattice constant was known to be a = 250 nm, the bandgap would be at f = 600 THz or wavelength λ = 500 nm.

MPB specializes in frequency solving, especially for band structures, in PCs [29,33]. Central to both MEEP and MPB is the ’geometry’ variable, which predictably specifies the physical geometry simulated. In MEEP, it primarily affects how electromagnetic waves propagate through the cell, while in MPB, combined with the ’geometry lattice’ variable, it specifies the periodicity of photonic crystals. The Brillouin zone of the crystal is described through the ’k-points’ variable in values of 2π/a.

2.2. Three-Dimensional Systems

The system chosen for this investigation was fcc lattices of packed spheres in general and inverse opal structures [34] in particular. A silicon-based inverse opal structure with air sphere radius r/a = 0.34 was used as the primary structure to determine how randomized variance in radius influenced the bandgap. Simulations of three-dimensional supercells were highly constrained by processing power. The time taken to compute a single band diagram is strongly dependent on the number of bands calculated, and the number of bands that need to be calculated significantly increases based on the size of the supercell. The silicon-based inverse opal supercell used was 2 x 2 x 2, for a total of 8 air spheres in the silicon material, enabling a larger number of simulations and this a larger sample size, within reasonable constraints of computing time. However, of note, the small size of the supercell (2 x 2 x 2) is nonetheless a good representation of all the nearest neighbor lattice interactions. An initial calibration cell was calculated with constant air sphere radius r/a = 0.34 to determine the number of bands needed to be calculated and was found to be about 150 bands.

A Gaussian distribution was used to allocate the radii of the air spheres within each independent simulation run. In order to maintain physically realistic outputs, the upper limit was set to radius r/a = 1/√8 ≈ 0.354 so the spheres cannot spatially overlap, with a lower distribution limit of r/a = 0.30. The median radius was selected to be μ /a = 0.34 with a standard deviation σ. These parameters were used to create an array of eight values, which was then used to overwrite the sphere radii in the original supercell, with the randomly assigned radii thus updated for each successive simulation run. Hence, outputs are collected in a “Monte Carlo” fashion from many multiple independent simulations.

Multiple standard deviation values were investigated to compare how the deviation of sphere radii, representing the physical disordering, affects the PBG. The target number of simulations for each value of σ was ~500 to produce statistically significant results. The largest deviation was chosen to be σ/a = 0.020, approximately 5.9% of the mean 0.34a. The smallest deviation was σ/a = 0.0025 (0.7% of the mean). A key range of interest for values in between was chosen based on the results produced, with a focus on deviations between σ/a = 0.0075, ~2.4% of the mean, and σ/a = 0.010, ~2.9% of the mean.

3. Results

3.1. Uniform Cells

As described, the calculations were validated using a unit cell with a known bandgap [1] for a given radius r/a = 0.34 and comparing to the uniform 2 x 2 x 2 supercell. The matrix material has dielectric constant ε = 13 (silicon) and the spheres ε = 1 (air). These two baseline simulations produced identical results to within three significant figures, a bandgap of 3.50% in frequency. Variation beyond these significant figures is attributable to the inherent uncertainty in the calculations. Identical simulations with no parameters changed showed variation beyond the three significant figures. Figure 2 shows the band structures, with Figure 2a showing the unit cell and Figure 2b showing the supercell. The supercell has significant band folding, resulting in a very dense band diagram compared to the unit cell.

3.2. Radii and Bandgap Distributions

All radii values produced for each standard deviation of radii σ were plotted in Figure 3 in a fifty-bin histogram, within the physically limited range of 0.30 < r/a ≤ 0.354 to show the overall distribution for each deviation. Each simulation produced eight radii values, and the overall shapes align with the expected Gaussian distributions. The narrower distributions are notably sharper than the broader distributions, resulting in a non-uniform y-axis scaling.

The bandgap magnitudes (by bandgap percentage) were likewise sorted into fifty-bin histograms, plotted in Figure 4. These were plotted against the frequency of simulation occurrences. A dotted black line marks the size of the uniform cell bandgap of 3.5%. The frequency of the zero bin (bandgap = 0%) increases sharply as σ increases, resulting in very “squashed”-appearing plots. A Gaussian line-fit, shown as solid blue lines, was applied to the histograms, excluding the zero bin, for all plots that showed a clear Gaussian shape. Figure 4h, σ/a = 0.015, and Figure 4i, σ/a = 0.020, did not fit such a Gaussian curve, and the mean bandgap percentage, Μ, is notionally ≤ 0. The parameters of the Gaussian fits (Μ and the standard deviation Σ) are collated in Table 1.

The MPB documentation [33] states that bandgaps of less than 1.0% are considered within the bounds of uncertainty and can represent false gaps produced by uncertainty in the calculation or the crossing of two frequency bands. This was not considered in plotting the histograms or in Table 1, although the general reliability of Gaussian fitting in Figure 4 suggest this is not a significant problem. The MPB documentation does suggest ways of interrogating suspected false gaps by targeting the specific frequency bands individually, but this was not done due to the large number of simulations produced.

In terms of the general trends observed for these inverse opal supercell lattices, the mean bandgap magnitude by percentage (Μ) increased with decreasing disorder from sphere polydispersity (σ) while the corresponding standard deviation (Σ) decreased with decreasing σ, shown in Table 1 and Figure 5. The errors in Μ, Σ were calculated using a standard error, where the mean and standard deviation of the Gaussian fit were divided by the square root of the number of simulations. A simple linear extrapolation of the plots in Figure 5a gives an inferred Μ = 0 intercept between σ/a = 0.013 and 0.015, and corresponding Σ value exceeding 1.5%. This is strongly corroborated in Figure 5c, where the cumulative fraction of simulations with non-zero bandgap (Ω) can be intuitively fitted to a sigmoidal function, with Ω tending to 1 as σ/a approaches zero and tending to 0 as σ/a becomes large. The half-height (Μ = 0.5) point of the trendline is found to be at σ/a = 0.014, determining a point at which there is quite a rapid pseudo-transition from a regime whereby bandgaps are mostly present in simulations to one where they are predominantly not. Moving into the distribution tail, the Μ = 0.1 (10%) point is at σ/a ≈ 0.021, and the M = 0.01 (1%) point at σ/a ≈ 0.028.

To complete this analysis, in Figure 5b, the fraction of individual simulations (Λ) with bandgap magnitude greater than the 3.5% crystalline baseline is plotted versus increasing sphere polydispersity (σ/a). For σ/a below ~0.005, these cases are widespread, but occurrences decline rapidly, thus conveniently fitting an exponential trendline, with very few simulations above σ/a = 0.015 exceeding the 3.5% benchmark.

3.3. Band Diagram Types

Each simulation run has a corresponding stored set of 8 sphere radii within the supercell, enabling a deterministic reconstruction of specific outcomes, relating to the photonic bandgap and dispersion relation properties.

As an illustration of this build-up of a “library” of designs, four major groupings of band diagrams (and associated PBG magnitudes) were produced across all simulation runs; the original ordered cell produced a bandgap (${\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{u}\mathrm{d}\mathrm{e}\mathsf{\Delta }}_0$) of 3.5%. Figure 6a is an example from an individual simulation output, of where the diagram shows no bandgap present ($\mathsf{\Delta }=0)$.. Figure 6b is an example of a small bandgap, with a bandgap size percentage of 1-3% ($\mathsf{\Delta }$< ${\mathsf{\Delta }}_0$). Figure 6c shows a comparably sized bandgap, ($\mathsf{\Delta }$≈ ${\mathsf{\Delta }}_0$,size of 3-4%). Figure 6d shows an example of a diagram with a bandgap size larger than the ordered cell ($\mathsf{\Delta }$ ≥ ${\mathsf{\Delta }}_0$) of 4-5+%.

3.4. Simulation Scaling

MPB scales significantly in time for larger lattice sizes, more bands, and higher resolution. Parallelization for both is limited to distributed memory systems such as CPU clusters through OpenMPI. OpenMPI showed speed increases over the non-OpenMPI simulation runs, especially for this system where multiple runs of the same simulation are used. The speed increase is described in Appendix B.

4. Discussion

In our supercell lattices, the distribution of random radii chosen was relatively narrow to determine how even small degrees of polydispersity affects the bandgap. It is demonstrable that randomization can also lead to the creation of bandgaps greater in magnitude than the uniform ordered system and is an interesting path for future research into, and combinatorial engineering of, three-dimensional photonic crystal systems.

In this pertinent exemplar system of fcc inverse opal lattices in silicon-air, Monte-Carlo based simulations of disorder from sphere radius polydispersity (σ) yield insightful qualitative and quantitative analyses. The imposition of a Gaussian distribution variation in the sphere radii within supercells gives a corresponding Gaussian distribution in the resultant photonic bandgaps, with the mean bandgap magnitude decreasing from the 3.5% baseline value for an ordered supercell with increasing σ, and the statistical distribution of bandgaps broadening such that a fraction of simulations no longer exhibits a complete bandgap in dispersion relations. Above a polydispersity of σ/a ≈ 0.014, a pseudo-transition characterized by a sigmoidal trend occurs into a regime whereby only a rapidly diminishing minority of simulations still returns a finite bandgap. In this high σ regime, there is an exponential tailing of occurrences with a bandgap greater than the 3.5% benchmark. Some isolated configurations, with a high degree of uniqueness, can exhibit enhanced bandgap properties despite considerable random positional disordering.

Future research and simulations into these photonic crystals should focus on scaling the lattice size. This will allow simulations to iterate towards genuinely randomly disordered structures, not constrained by supercell size or periodic boundary conditions. Additionally, we might consider the effects of random deviations from ordering on other photonic crystal symmetries and custom lattices through packing geometry. Systems with inherent disorder, such as those with random packing arising from the lack of next-nearest-neighbor interactions, may be of particular interest [35] for simulation.

5. Conclusions

The simulation tools MEEP and MPB were used to study how randomization of the structural order in photonic crystals has a significant impact on photonic properties and bandgaps. By understanding the interplay between lattice disorder and photonic bandgap properties, researchers can develop strategies to mitigate or implement disorder-induced effects and enhance both the functionality and the performance of photonic devices. Overall, the study contributes to the understanding of disorder effects on PBGs and illustrates the potential use of engineered randomness in supercell systems to create targeted photonic crystal properties and functionality.

Recent advances in computing and methods such as machine learning will conceivably allow for more in-depth FDTD simulations and modelling of engineered disorder in photonic crystals, offering avenues for further prescient research into designing and optimizing photonic structures for practical applications.