Photonic crystals (PCs) are materials that alter the dispersion properties of light by means of a periodic modulation in their dielectric constant [1, 2]. It is well established that PCs based on certain structures, like diamond [3, 4], can suppress the omnidirectional propagation of light over a range of frequencies. This range of frequencies for which the photon density of states drops to zero is called a photonic band gap (PBG) and has been the driving force behind a large effort to design periodic dielectric structures that act like semiconductors for photons. Such materials will likely serve as platforms for the integration and miniaturization of novel optical devices like lossless waveguides [5], all-optical switches [6], and zero-threshold lasers [1]. Various techniques to fabricate these materials have been developed, such as self-assembly [7], layer-by-layer lithography [8], two-photon laser writing [9], and holographic lithography (HL) [10, 11].

Holographic lithography, introduced by Berger et al. [10] and extended to three-dimensional (3D) structures by Campbell et al. [11], has emerged as a leading candidate for the versatile fabrication of PBG materials due to its ability to generate an unparalleled variety of thick, defect-free structures that span large areas. Using the interference profile generated by the overlap of four or more non-coplanar laser beams, HL produces a 3D structure within a film of photoresist using just a single exposure. The resultant structure may then act as a template for the infiltration of high-k dielectrics using techniques such as chemical vapor deposition [12–14] and/or atomic layer deposition [15, 16]. These infiltration methods can be modified so that either the polymer template or the vacant background is replaced with a high dielectric constant material. The final structure, having an increased dielectric contrast, will then possess a PBG provided the underlying structural morphology is favorable.

Thus, the design of holographic beam setups that yield structures with large PBGs has been the focus of much research [17–19]. Predominantly, such efforts have relied on analytical approaches based on generating a desired set of symmetries within a resultant structure. While symmetry considerations have played an important role in the design of PBG structures, the mere presence of a particular symmetry group does not strictly guarantee the existence of a PBG. For example, a face-centered cubic (fcc) arrangement of dielectric spheres does not possess a PBG for any dielectric contrast, yet its inverse structure, having the same symmetry group, does. In this sense, a full account of the structure’s morphology matters more than just the symmetry it possesses. This notion is further supported by the various approximations to large-gap structures that possess sizable gaps of there own, despite not retaining the full set of symmetries of the “target” structure [4]. Most notably, the layer-by- layer “woodpile” structure [20], which garners a gap comparable to that of diamond, mimics the morphology of the diamond structure without possessing the same symmetry [4]. With this in mind, analytical design strategies based solely on obtaining a desired set of symmetries patently overlook robust approximations to large PBG structures.

Recently, computational approaches have emerged as an alternative to analytical design strategies. However, their scope has been limited to the optimization of intensity contrast [21] or the simple scanning of small subsets of beam parameters to optimize for gap size [22]. Nonetheless, computational approaches have potential to more thoroughly explore the entire parameter-space, in contrast to analytical strategies which limit themselves to solutions of prescribed symmetry. By lifting this constraint, an effective computational approach may potentially unearth holographic structures that possess little or no symmetry, but are nonetheless close morphological approximates to desirable target structures, particularly those having large PBGs.

In this paper, we describe a versatile computational approach that utilizes genetic algorithms (GAs) to design holographic structures having large 3D band gaps. GAs are inherently parallel search techniques based on the principles of evolution [23, 24]. They utilize the biologically analogous operators of selection, crossover, and mutation to operate on abstract representations of trial solutions. Over time, these operators drive the trial solutions to converge toward a desired objective. For a wide variety of complex optimization problems GAs excel at finding robust approximations because they adeptly avoid convergence to local optima and do not rely on initial starting guesses. They have been successfully implemented in various problems in photonics, such as photonic crystal design [25, 26] and waveguide demultiplexing [27], indicating their viability as an inverse design method for HL.

In this work, a GA was used to search for holographically definable structures that approximate the real-space architecture of the rod-connected diamond (RCD) [28], shown in Fig. 1(a). Because of the structure’s exceptionally large gap of 30% - the largest reported for a silicon-air structure (ε_silicon=12.96) - the design of holographic architectures with complete diamond symmetry has been a focal point [18, 19]. However, the only setup known to produce such a structure calls for the beams to be launched from opposite hemispheres, requiring the use of a transparent substrate. From a fabrication standpoint, this is unfortunate because many of the promising applications for integrated optical circuitry require PBG materials on opaque or absorbing substrates like silicon. Here, we report a holographic design where twelve beams originating from the same hemisphere produce a high-contrast structure that very closely mimics RCD. At a dielectric fill of 17%, the resultant two-component structure shown in Fig. 1(b) has a gap/midgap ratio (Δω/ω) of 28%, the largest yet reported for a holographically definable structure.

 

Fig. 1. (a). Unit cell of the RCD structure having a PBG of 30%. (b). Best holographically definable structure approximating RCD. A PBG of 28% is obtained for a dielectric fill of 17%.

Download Full Size | PDF

Because we require all beams to originate from the same hemisphere, the wave vectors, k_i, are chosen prior to the simulation and universally applied to a “population” of trial solutions. In spherical coordinates they are expressed as k_i=2π/λ (sin θ_i cos φ_i, sin θ_i sin φ_i, cos θ_i), where λ is the wave length, θ_i is the polar angle, and φ_i is the azimuth. To ensure a lattice coincident with that of the desired structure, the differences between wave vectors, written as g_ij=k_i-k_j, must include a set of primitive vectors corresponding to the reciprocal lattice of the target. The remaining g_ij must be linear combinations of these primitive vectors. Here, we choose a twelve-beam configuration that produces an fcc lattice coinciding with that of RCD. For this arrangement, θ_i=cos^-1(7/$\sqrt{65}$) and φ_i=iπ/2 when 0≤i≤3. For the remaining beams, θ_i=cos^-1 (5/$\sqrt{65}$) and φ_i=[tan^-1(-3), tan^-1(-1/3), tan^-1(1/3), tan^-1(3), π+tan^-1(-3), π+tan^-1(-1/3), π+tan^-1(1/3), π+tan^-1(3)] for i=[4, 5, 6, … 11] respectively.

Since the trial solutions and target structure have matching lattices, we need only to approximate the morphology within the primitive cell of RCD. However, due to the convenience of working with an orthogonal basis, we have chosen instead to focus on its unit cell shown it Fig. 1(a). To implement RCD as a “target,” we discretize its unit cell into a 100×100×100 binary array of low (ε_air=1.00) and high (ε_silicon=12.96) dielectric material. We choose the high dielectric material to comprise ≈20% of the volume, roughly corresponding to the fill ratio at which the largest gap occurs. This array is then used to numerically evaluate a given holographic structure on a voxel by voxel basis.

In general, the 3D intensity distribution generated by the interference of N mutually coherent, monochromatic plane waves of amplitude E_i and phase δ_i, is given by

where ∂_ij=δ_i-δ_j, and ê_i denotes the complex unit polarization vector describing elliptically polarized light. We choose to express ê_i in terms of two mutually orthogonal, linearly polarized plane waves out of phase by χ_i and having a ratio in their amplitudes given by E_iy/E_ix=tan γ_i, where γ_i is the angle between the locally defined x-axis and the vector sum between the two orthogonal fields [29]. Explicitly, ê_i is expressed in Eq. (2) as the locally defined Jones vector, J, transformed to global coordinates by the rotation matrix R and is written as,

The structure produced in photoresist is determined by applying an intensity threshold, I _th, to Eq. (1). For negative resists, regions of intensity greater than I _th polymerize into a crosslinked solid phase, while unpolymerized regions less than I _th are dissolved away. The holographic structure is formally expressed by the Heaviside step function,

where “1” typically represents regions where ε=12.96 (silicon) and “0” represents regions where ε=1.00 (air).

In all, there are 6N+2 parameters that determine the holographic structure described by Eq. (3). But since we have predetermined k_i, only 4N+1 parameters remain to be optimized. These parameters are encoded into a binary string, called a “chromosome,” which is the basic entity subject to evolution during the simulation. In general, a parameter, y, lying on the interval [α, β] can be encoded as a binary “gene,” y’, of length l according to the transformation y’=(y-α)/(β-α)×2^l. Here, l, is made just large enough (7–9 bits per parameter) to provide sufficient resolution of the parameter space. Since the layout of the encoded parameters within the chromosome can affect a GA’s performance, we have chosen to cluster important “genetic material” together such that potent binary sequences are less likely to be disrupted during crossover. Figure 2 illustrates our method for distributing “genetic material” within the chromosome.

 

Fig. 2. Bit layout for the binary parameters, y’_ij, within the chromosome. A specific case for which each parameter has a length of l=7 is used as an example. The subscripts i and j correspond to the beam and bit numbers respectively. For clarity, each bit having the same subscript j is labeled in the same shade of grey, where darker shades correspond to larger place values. The j ^th bit of each parameter is placed in the corresponding j ^th substring according to the arrangement shown at the bottom (for j=2). The substrings are then placed into the chromosome (right) using an alternating stacking method so that bits with a higher place value are clustered together.

Download Full Size | PDF

At the start of a simulation, N_c chromosomes are randomly generated to form a population of trial solutions (in this paper N_c=6500). The corresponding holographic structures are calculated according to Eq. (3) using an adjustable I _th numerically chosen to constrain each structure’s fill ratio to ≈20%. Each chromosome is then assigned a fitness, F, equal to the fraction of voxels in one unit cell that match the target. The top 1% are copied directly to the next generation unaltered. This procedure, referred to as elitism, improves a GA’s performance by preserving the “elite” chromosomes, allowing them to reproduce more often than the non-elite.

To produce the subsequent generation, the following 3 steps are repeated until N_c “offspring” are produced: (1) Two chromosomes are selected to crossover with a selection probability, p_s∝F^x (here x=5 provided optimal results). (2) Once selected, the two chromosomes undergo crossover with a probability p_c=0.95. Crossover then occurs by exchanging equivalent segments of the chromosomes about randomly chosen nodes to create two new chromosomes. In our simulations, we employ a simple single-node crossover. (3) Finally, each bit in the resultant offspring is mutated (flipped) with a probability p_m=0.01. Once N_c offspring have been generated, the “parent” generation is replaced by the newly formed offspring. The same procedure is then repeated for all subsequent generations. The values assigned to N_c, p_s, p_c, p_m, and x were either empirically motivated or based on the findings in Ref. [30].

Because GAs cannot guarantee global convergence, we carried out six simulations with different randomly generated starting populations to increase the chances of a successful run. We allowed each simulation to iterate through 300 generations before stopping to save the best chromosome in each. The best fitness is plotted as a function of generation in Fig. 3 for the best run out of the batch of six. Upon completion, the globally optimized solution yielded a 94% match to the target. The parameters encoded in this chromosome were converted into floating point numbers and underwent local optimization using a next ascent hill-climbing (NAHC) algorithm. The NAHC algorithm operates by systematically varying each parameter while keeping those changes resulting in an increased F. When F no longer increases, the size of the variation is decreased by half until a minimum step-size is reached. Because NAHC is highly susceptible to convergence at local optima, multiple runs are performed in series using RCD targets of different fill ratios to prevent local stagnation. In this work, we used various fill ratios ranging between 15% and 21%. The final optimized parameters, listed in Table 1, yield a structure with F=97%. The resultant structure bears a strong resemblance to that of RCD, as seen in Fig. 1.

 

Fig. 3. The best fitness plotted as a function of generation (top). The unit cells (bottom) show the (111) view of the best structures at different generations during the run. Letters a - e match each unit cell to a different point on the graph. The unit cell of the RCD structure, which was used as the target, is shown at the bottom right for comparison.

Download Full Size | PDF

Table 1. List of optimized parameters used to make structure shown in Fig. 1(b). All angles are given in radians.

View Table

From an experimental viewpoint, this design possesses other advantages in addition the unilateral arrangement of beams. In particular, the design yields a large contrast in the resulting intensity profile, I _max/I _min≈125. This is an important consideration for the high-fidelity transfer of the intensity profile into photoresist. Equally important is the bicontinuity of the resultant structure. Bicontinuous topologies enable the dissolution of all underexposed resist and ensure that the resultant structure is self-supporting upon the removal of unsolidified material. Our design possesses a large range of intensity thresholds, I _th, for which a bicontinuous structure is formed. We report a usable range of I _th that spans 59% of the possible intensity range, (I _max-I _min). This corresponds to a range of filling fractions between 6% and 87%. In addition, the field amplitudes listed in Table 1 suggest that several beams can be eliminated to simplify the optical setup. It appears as though the subset of seven beams for which E_i>0.15 will suffice in obtaining a structure close to that of Fig. 1(b). Details of this will be presented elsewhere [31].

Using the MIT Photonic-Bands package [32], we optimized the normalized gap width for the family of structures based on the parameters listed in Table 1. Figure 4(a) shows the normalized gap width as a function of filling fraction (which is adjusted via I _th) and refractive index contrast. Each point in Fig. 4(a) was calculated using a mesh of 32×32×32 and 773 k-points distributed throughout the first Brillouin Zone (FBZ). For silicon-air contrast, the largest gap was found at 17% dielectric fill, corresponding to a threshold intensity, I _th=0.8×I _max. The density of states for this gap-optimized structure is plotted in Fig. 4(b) using the same resolution and 7597 k-points. Using a 64×64×64 mesh and 1298 k-points distributed on the surface of the FBZ, a final gap width of 28% was verified.

 

Fig. 4. (a) Gap width between the 2nd and 3rd bands plotted as a function of dielectric fill and dielectric contrast for the structure shown in Fig. 1(b). Negative gap widths correspond to the 3rd band dropping below the 2nd. At a dielectric fill of 17%, the optimized holographic structure shown in Fig. 1(b) achieves a 28% gap for silicon-air dielectric contrast (ε=12.96). (b) Density of states, g(ν), for the corresponding holographic structure shown in Fig. 1(b) having a dielectric fill of 17%.

Download Full Size | PDF

In summary, we have demonstrated the use of genetic algorithms as a comprehensive approach to the design of holographic structures. In contrast to analytical methods that rely on symmetry considerations, this alternate approach allows for approximations to real-space structures, opening up a largely unexplored parameter-space for the design of holographic architectures. Using this technique, we have discovered beam parameters that produce a diamondlike structure having a PBG of 28%. The design’s unilateral beam geometry should allow the resultant structure to be fabricated on silicon substrates, making it a candidate platform for a host of all-optical devices integrated into existing silicon-based technologies.

We are grateful to Robert Shimmin and Vinayak Ramanan for fruitful discussions. This work is funded by the Army Research Office through a MURI grant (DAAD19-03-1-0227).