Stochastic optimization of broadband reflecting photonic structures.

Photonic crystals (PCs) are built to control the propagation of light within their structure. These can be used for an assortment of applications where custom designed devices are of interest. Among them, one-dimensional PCs can be produced to achieve the reflection of specific and broad wavelength ranges. However, their design and fabrication are challenging due to the diversity of periodic arrangement and layer configuration that each different PC needs. In this study, we present a framework to design high reflecting PCs for any desired wavelength range. Our method combines three stochastic optimization algorithms (Random Search, Particle Swarm Optimization and Simulated Annealing) along with a reduced space-search methodology to obtain a custom and optimized PC configuration. The optimization procedure is evaluated through theoretical reflectance spectra calculated by using the Equispaced Thickness Method, which improves the simulations due to the consideration of incoherent light transmission. We prove the viability of our procedure by fabricating different reflecting PCs made of porous silicon and obtain good agreement between experiment and theory using a merit function. With this methodology, diverse reflecting PCs can be designed for any applications and fabricated with different materials.

Introduction

The unique optical features of confining light or controlling the propagation of electromagnetic waves in photonic crystals (PCs) enable many important device applications. Waveguides[1], Bragg-mirrors[2], surface plasmonic structures[3], nanocavities and nano-lasers[4,5] are some examples of optical devices that can be fabricated using PCs. These structures have in common that they are formed of periodically arranged dielectric materials, and their functionality depends highly on the distribution of those materials. Photons entering a PC interact with its periodic dielectric constant and are consequently organized into photonic bands. Analogously to electrons in a crystal, their propagation will be limited by photonic band gaps (PBGs) where transmission states are forbidden. Hence, different designed PCs constrain the propagation of the photons (electromagnetic waves) within the structure depending on their periodic arrangements.

The simplest PCs are one-dimensional and are composed of alternating high and low refractive index layers. These structures can be designed to interact with the electromagnetic waves in such a manner that photonic mirrors (also called Bragg-reflectors) are formed, where “perfect” reflection for specific wavelengths can be achieved. Broadband mirrors can also be engineered if the reflection of a broader wavelength range is needed. These broader mirrors can be achieved by overlaying different Bragg reflectors where each one reflects around a specific central wavelength, resulting in a wider PBG. However, choosing the proper central wavelengths and hence the layer configuration for these Bragg reflectors to superpose adequately is no simple task; furthermore designing such reflecting structures may involve complicated and unpractical methods.

In the context of this work, previous reports have been focused on the PBG enlargement based on staggered structures[6,7] or chirped gratings[8-10], amongst others. Commonly, predetermined structure parameters such as layer number, the amount of Bragg reflectors or layer thickness functions constraint to experimental parameters, limit the versatility of these methods. In this paper, we present a practical method to custom design and then evaluate the fabricated broadband PC mirrors. This procedure can be applied to different dielectric materials and easily modified to design other PC structures such as filters for example.

Generally selecting the optimal number of layers and determining their thickness, related to the reflectance wavelength range, represents the main difficulty when designing broadband reflecting PCs. For practical purposes it is pursued that only the wavelength range for the desired application needs to be defined; also, the size of the multilayer correlates directly to the fabrication time which is sought to be the lowest possible. Hence, in this work, we present a stochastic optimization method that finds the best PC layer configuration to build any reflecting structure within the selected wavelength range. This procedure is based on evaluating all possible PCs using their theoretical reflectance spectra and choosing the most reflective PC as the optimized structure. Here we used the equispaced thickness method[11] to simulate the PCs reflectance spectra considering the effect of incoherent light transmission through the layers based on the transfer matrix formalism. With this approach, we are able to model more realistic multilayered structures which may present imperfections or interface roughness of the fabricated structures leading to scattering effects which need to be considered.

Given the above requirements, to design high quality PCs, topology optimization must be conducted. When addressing topology optimization, previous studies have focused on the use of gradient-based algorithms through sensitivity analysis[12-15]. Moreover, it should be noted that in previous reports, the optimization was done over linearized PCs for short bandwidths. The optimization of PCs for short bandwidths does not present the extra complexity of being multimodal (i.e. having many local minima), the current study explores a much wider range of bandwidths which presents severe nonlinearities resulting in multimodality of the solution space (i.e. many local minima). While gradient-based algorithms are efficient in many cases, they present two main drawbacks for the current research. First, gradient-based algorithms are of a local search nature, which is ideal when confronted with convex optimization tasks; however, the problem presented in this study is highly nonlinear and nonconvex. Second, the combination of derivative calculations and derivative evaluations is computationally intensive for the current application. Furthermore given that the simulations of the reflectance spectra consist of functions with sharp gradients (e.g. exponential, cosines, sines), first order approximations might be ineffective and can lead to convergence problems from the line search[16] framework. Therefore, the topological optimization in this research is performed through a stochastic optimization algorithm particularly tailored for the current application.

In order to test the performance of this optimization procedure, one dimensional PCs structures need to be fabricated. Production of one dimensional PCs can be achieved using several experimental techniques limited to the fabrication of multilayers with large refractive index contrast. Porous silicon (P-Si) is a nanostructured material widely used to fabricate one-dimensional PCs due to the tunability of its porosity and simple fabrication procedure. In this study, we prove the use and viability of our optimization method by fabricating three different reflecting PCs and using a merit function we measure the differences between theory and experiment finding excellent agreement.

The paper is composed as follows. In section 1 we present a mode to calculate the central wavelength distribution of the Bragg mirrors which form a broadband mirror. In section 2 the stochastic optimization method developed for the search of an optimized reflecting PC structure is shown. The results of the simulations and the fabrication of the optimized broadband mirrors are shown in section 3; where the agreement between theory and simulations is also discussed and evaluated using a merit function. In section 4 the concluding remarks are presented followed by the Methods section where the calculations of the theoretical reflectance spectra, the optimization algorithm and the fabrication of P-Si are described.

Central wavelength distribution

Designing a broadband photonic mirror requires an adequate superposition of individual Bragg reflectors (named submirrors) that conform it. Each submirror (n = 1, …, n, …, f ) creates a PBG that forbids EM transmission around the central wavelength (Λ(n)) it is designed for. In this manner, to reflect over a wide wavelength range it is necessary to find an appropriate distribution of Λ(n) for the submirrors to superpose in an optimized way. In recent work, we proposed a method based on the Padé approximant to find each central wavelength[17], where we defined its distribution as:here a, a1 and b1 are unknown coefficients; to find them we propose the following system of equations:

The central wavelength of the first Λ(1) and the last Λ(f) submirrors determine the reflected wavelength range that is required and consequently define the PBG of the broadband mirror. To solve the system of equations 2 we need to know which submirror is the intermediate n and its corresponding central wavelength Λ(n). These two parameters represent two of the optimization parameters to find the best PC structure to construct a broadband mirror.

The alternating layers that form each submirror need to satisfy the quarter wavelength condition which relates the layer thickness d and its refractive index η to the central wavelength as ηd = Λ(n)/4. The number of layers that form the PC structure is also important to optimize because it defines the width of the entire multilayer and correlates directly to the fabrication time. In general, if the number of submirrors is increased the reflectance spectrum of the broadband mirror presents fewer oscillations, and the quality of the PC reflector is enhanced. However, if the structure is formed of too many layers, the fabrication time increases substantially achieving an ineffective manufacturing process. For example, when fabricating PC mirrors with P-Si and the etching time is large (>2 hours) a porosity gradient in the layers is formed due to the change of the electrolyte concentration. This effect generates a widening of the PBG towards lower frequencies, and although we have boarded this issue previously in a qualitatively way[17], we will not address this matter in this work. Nevertheless, we need an adequate number of submirrors to design and fabricate optimized PC mirrors. Thus we finally consider three optimization parameters in our stochastic optimization method: Λ(n), n and n. Each set of these parameters represent the coordinates of particles in a three-dimensional space (the solution space) for all possible PCs at a defined wavelength range. Here we search for the particle with the most reflective spectrum by using the stochastic optimization method presented below. At the end, we use the coordinates of this particle to fabricate the optimized PC mirror with P-Si and prove the usability of our procedure.

Stochastic optimization Method

The optimization task performed in this work is to find the best possible PC layer configuration with a maximum reflectivity in the selected wavelength range. To this aim a stochastic optimization algorithm is used to explore the solution space of all possible PCs, with respect to the performance criteria presented in equation 5 (obtained in the Methods section): . The highest value of RP relates to the maximum reflectivity of the PC mirror that we want to optimize through an algorithm. The algorithm developed for this purpose should possess good exploration and exploitation capabilities to seek an optimal solution efficiently. In simple words, the exploration process refers to the ability of an optimization algorithm to seek for the global optimum in the solution space of an unknown optimization problem, while the exploitation process refers to its ability of applying the knowledge of previous and current solutions to look for better ones. However, the processes of the exploration and the exploitation have conflicting goals, the two must be well balanced to achieve a good optimization[16]. These processes can be better understood when a gradient-based approach and a random search algorithm are compared. On one hand, a gradient-based approach has great exploitation capabilities, as it very well exploits the information of previous and current solutions, however, it easily gets trapped by local minima. On the other hand, a random search algorithm can find a global solution if space is searched widely enough, but its complete lack of exploitation makes it an inefficient algorithm in general.

To address the aforementioned challenges and obtain good optimization performance, three stochastic optimization methods are combined to obtain a hybrid algorithm. This hybrid algorithm uses a number N of particles with coordinates (Λ(n), n, n) as its agent and searches the solution space through a combined form of: a) Random Search (RS), b) Particle Swarm Optimization (PSO) and c) Simulated Annealing (SA). Also, when the algorithm determines that there are no good solutions outside a specific domain, the search space is reduced to focus efforts in areas which are more likely to have a high-quality solution. The algorithm, as well as the reasons behind this combination, are described in the Methods section.

Results

In this study, we present three different optimized reflecting PCs each one designed to inhibit electromagnetic wave transmission over diverse and wide wavelength ranges. The first mirror M1 is designed to reflect from Λ(1) = 400 to Λ(f) = 2000 nm, the second M2 from Λ(1) = 600 to Λ(f) = 1200 nm and the third M3 from Λ(1) = 800 to Λ(f) = 1800 nm. We restrict our study only for p = 5, which seems to be an adequate number of periods and hence used this value for all the calculations. We ran the optimization algorithm twice for each wavelength range and obtained slight different optimized particles (each one with different coordinates) for all mirrors but with very similar RP values between every pair (A and B). Thus the optimization algorithm finds optimal PC configurations in every search, even though they have small differences in their central wavelength distribution. The optimized parameters of the three pairs of mirrors are summarized in Table 1 where the similarities can be compared. Given these values we solved equation 1 and found the central wavelength distribution Λ(n) of each PC configuration. Insets in Figs 1–3 are the Λ(n) of mirrors A and B respectively where no significant differences can be observed. From the quarter wavelength condition, we were able to calculate the layer thicknesses of the submirrors which form the PC structures and afterward we used this information to finally fabricate the mirrors.

Table 1

Optimization parameters (Λ(n), n, n) of the reflecting PCs that we designed in this study which were obtained from the stochastic optimization method. Their reflectance performance criteria RP are also enlisted.

PC	Λ(ni) (nm)	n i	n s	RP
M 1A	966	5	12	136’506
M 1B	1114	5	12	136’507
M 2A	838	3	8	55’022
M 2B	929	4	8	55’021
M 3A	1076	3	8	94’521
M 3B	1508	6	8	94’519

Figure 1

Theoretical and Experimental Reflectance spectra of mirrors M1.Theoretical (T, T) and Experimental (E and E) reflectance spectra of broadband mirrors A and B designed to reflect from 400 to 2000 nm. Their central wavelength distribution Λ(n) (inset figure) are obtained from the optimized parameters.

Figure 3

Theoretical and Experimental Reflectance spectra of mirrors M3.Theoretical (T, T) and Experimental (E and E) reflectance spectra of broadband mirrors A and B designed to reflect from 800 to 1800 nm. Their central wavelength distribution Λ(n) (inset figure) are obtained from the optimized parameters.

Theoretical and Experimental Reflectance spectra of mirrors M2.Theoretical (T, T) and Experimental (E and E) reflectance spectra of broadband mirrors A and B designed to reflect from 600 to 1200 nm. Their central wavelength distribution Λ(n) (inset figure) are obtained from the optimized parameters.

Once the layer thicknesses of each PC configuration were defined, we fabricated all the aforementioned mirrors with P-Si following the fabrication procedure described in the Methods section. The experimental and theoretical spectra of each pair of mirrors are compared in Figs 1–3 where no differences between the theoretical spectra are observed, but slight disagreements in the experimental counterparts are found. In general, there is a very good agreement between simulations and experiments, but we need a quantitative measure to calculate the differences which we present in the next section.

Discussion

The theoretical reflectance spectra of the PC mirrors represent the performance criteria of each optimized structure. In practice these calculated spectra are attempted to be reproduced with the experiments; hence we validated the accuracy of the fabricated mirrors using a merit function as a quantitative quality measure between theory and experiment. Here we define the merit function aswhere T(λ) and E(λ) are the simulated and experimental reflectance spectra, respectively. The best agreement between experiments and theory requires N to be small. The optimization method that we present in this report is based on evaluating the most reflective theoretical spectra of each PC structure, furthermore, we use this method and find an adequate experimental agreement of theoretical and experimental spectra. We emphasize that to have quantitative comparison we calculated the merit function for all the mirrors that we fabricated within this study and find good agreement for all the samples. In Table 2 we resume the obtained values of N which show differences of order 10−5 between the mirrors A and B of M1 and M3 and differences of order 10−4 between the mirrors M2. Thus, these results show the viability of our optimized methodology to design and fabricate broadband high reflective PC structures.

Table 2

Values of the merit function N and fabrication time (t) of the reflecting PCs presented in this report.

PC	N (×10−3)	tf (h)
M 1A	4.41	1.78
M 1B	4.36	1.78
M 2A	3.83	0.93
M 2B	3.23	0.93
M 3A	1.11	1.30
M 3B	1.20	1.30

n class="Chemical">Values of the merit function N and fabrication time (t) of the reflecting n class="Chemical">PCs presented in this report.

Besides fabrication and reproduction quality, the time diminution of the process is important to be considered. In previous work, we have reported P-Si mirror structures which are as reflective as the optimized mirrors we present here, but their multilayer sizes are much larger because of their large submirror number (n = 20)[17]. These mirrors take 3.5 h at least to be fabricated which double the time used to produce the optimized structures presented in this paper (see Table 2). Therefore by using the stochastic optimization method to design and then fabricate the PC mirrors we have improved the efficiency of the custom mirror making process.

We stress that this optimization method which balances exploration and exploitation and combined with the solution space reduction can also be employed to optimize 2D or 3D PCs efficiently. The only necessary adaptation to optimize 2D and 3D structures would be to define an adequate form of constructing the possible PCs along with customized performance criteria for the desired application, as we did with the Padé central wavelength distribution for the 1D case.

Concluding Remarks

In this research, we have developed a practical design procedure to find the most reflective broadband PC within a desired wavelength range. First, we consider the alternating layers of each submirror formed with two materials of different refractive indexes, which can be wavelength dependent. Second, to form the broadband PC, we vary three design parameters: the central wavelength of each submirror, the distribution of these central wavelengths, and the number of Bragg mirrors to cover the complete wavelength range of reflection. To find an optimized high reflective PC, we combined three stochastic optimization methods and obtained a hybrid algorithm which furthermore includes a space reduction strategy to search for the optimized PC. Also, we validate this procedure by fabricating several PCs applying the optimization method, and we measure the agreement between theory and experiment by means of a merit function. The results show good agreement and also prove efficient fabrication times which add to the improvement of the optimization method to design reflecting PCs. We want to stress tree main contributions of this paper: first, the novel optimization algorithm constructed in this work can be used to optimize other highly nonlinear and nonconvex systems; second the validation against experimental performance of the optimized theoretical PCs; and third the results are not only valid for PC made with porous silicon multilayers. The last implies that with this procedure any reflecting PC structure can be designed for different applications and fabricated with either material adequate for PC manufacturing.

Methods

Simulation of theoretical reflectance spectra

One of the most common tools to calculate multilayer structures scattering spectra (transmission, reflection, and absorption) is the transfer-matrix method (TM)[18]. However, the form of the conventional TM assumes coherent light propagation through a plane and homogeneous layers which may result in unrealistic spectra calculations. In practice, layers formed of P-Si, for example, have nonuniform interfaces presenting rough surfaces as shown in the SEM image of Fig. 4. The interaction of light with these imperfections lead to destructive interference effects that can be introduced in the simulations by considering incoherent layers in the system. The equispaced thickness method (ETM) proposes a form to incorporate these effects by including equidistant variations in the thickness of the layers[11]. Using the TM formalism the transmission or reflectance spectra for each variation is calculated and then averaged over the equispaced thicknesses.

Figure 4

SEM image of mirror. SEM image of a porous silicon photonic mirror where the high and low refractive index layers are shown.

SEM image of mirror. SEM image of a porous n class="Chemical">silicon photonic mirror where the high and low refractive index layen class="Chemical">rs are shown.

In this study we use the ETM to simulate the reflectance spectra of each desired PC structure; moreover, these spectra represent the evaluation function of our stochastic optimization method. Within the TM theory, the electromagnetic field propagation in the structure formed of l layers is represented by the matrixwhere E, H, E(, H( are the electrical and magnetic fields in the first and the last interface, respectively. The transfer matrix is defined by the product of the characteristic matrices of each layer from j = 1, 2, … l as

For a transverse electric field at normal incidence, the transfer matrix is a function of the complex refractive index η, the thickness d of the layer j and the wavelength λ of the electric fieldHere Y is determined by , where ε0 and μ0 are the vacuum permittivity and the permeability.

Because the PC mirrors are formed of n number of submirrors (n = 1, ..., f ), with p number of periods composed of layers of high (H) and low(L) refractive indexes, the total transfer matrix for this structure is constructed by multiplying iteratively the submirror matrices

The reflection coefficient can be written in terms of the elements m of M, here Y stands for the incident media (air):from here the reflectance is obtained as

The reflectance spectrum of a PC multilayer is simulated by calculating R for all wavelengths within the defined range. To incorporate the equispaced thicknesses in the previous description we redefine the thickness of each submirror layer as d(n) = Λ(n)/4η(n) + δ where k = H, L and δ represents the variance of the thickness due to the surface roughness or inhomogeneities which cause incoherent light propagation. In this study, we define δ as a function of the measured layer rugosity S (see the Porous silicon fabrication Section for detailed explanation of the determination of S of P-Si) in the formwhere i takes values from 1 to X, the number of equidistant values. For each given value of i, we calculate the corresponding reflectance spectrum of the multilayer which forms the PC. Then we have a collection of X reflectance spectra and we average them over the equispaced thickness values to finally obtain the simulated spectrum of the PC mirror. In this study, we used X = 5 after finding that a larger number did not significantly change the simulations.

Following on we perform a numerical integration of the averaged reflectance spectrum to determine its reflectance performance criteria (RP) that is needed for the stochastic optimization method as the evaluation function. For the numerical integration we first apply a linear interpolation to the spectrum and then integrate over the wavelength range of the PC mirror:

In this manner we have defined the mathematical problem which is nonlinear and presents nonconvex features. Therefore we require an optimization method able to deal with these characteristics. The main objective is to evaluate RP for each different PC structure in the parametric three dimensional space comparing and classifying the optimization algorithm described below to find the most reflective PC which corresponds to the PC with the highest RP value.

Optimization algorithm

The algorithm proposed in this research uses three main phases: Initialization, Evaluation and Classification and Space Exploration.

Initialization

In this phase, N particles are randomly placed in the solution space of the optimization problem.

En class="Chemical">valuation and Clasn class="Chemical">sification

The value of RP for each different PC structure in the solution space is compared and classified accordingly to the optimization algorithm described below to find the most reflective PC which corresponds to the PC with the highest RP value. In this phase, the performance criteria of each particle is evaluated, after which particles are classified into three groups, and the group they are assigned to will determine the search strategy they will follow.

N particles with the highest en class="Chemical">valuated function will follow an n class="Chemical">SA search strategy for the following n iterations.

N particles with the lowest en class="Chemical">valuated function will follow an n class="Chemical">RS search strategy for the following n iterations.

The remaining N particles will follow a n class="Gene">PSO search strategy for the following n iterations.

Space Exploration

The above classification works in such a way that the position of each particle is better exploited depending on the information it can supply to the swarm as a whole.

The N particles should be a small number of particles which are assumed to be near high-quality solutions and therefore should intensively explore their neighborhood. To this end, SA, an efficient stochastic optimization algorithm of a local nature is used. The particular implementation in this work follows from[19].

A similar reasoning is applied to the N particles. In the solution space, these particles are positioned such that they might not be close to high-quality solutions, however by using the knowledge of other particles in the swarm they can search for better solution areas which have not yet been found by other particles. PSO is mainly designed for this purpose, as it is inspired in a flock of birds collectively exchanging information and foraging for food, likely to move closer to an optimum of a fitness function. The implementation in this work follows from[20].

Finally, the N particles are those assumed to be relatively far from any high-quality solution, and hence exploring the neighborhood around them would not bring benefit to the swarm. The No Free Lunch Theorem[21] establishes a theoretical proof that any optimization algorithm is in average as good as the Random Search strategy, and that the advantage of specialized algorithms is in using and obtaining knowledge about the problem at every iteration. Given that particles that are too far away from optimal solutions can hardly provide useful information, then particles classified as N are reinitialized randomly.

Repeat

After n iterations the algorithm either returns to the En class="Chemical">valuation and Clasn class="Chemical">sification phase or terminates.

Space reduction

During the algorithm, the best pon class="Chemical">sition of each particle in the current cycle is recorded. After N iterations, the search space is reduced to the smallest hypercube that contains the best pon class="Chemical">sition of all particles in the current cycle, this terminates the cycle. Subsequently, the algorithm returns to the Initialization phase.

In this way, given that RS, PSO, and SA are in ascending order of exploration and in descending order of exploitation a balance in the exploration-exploitation paradigm can be achieved. Furthermore, the space reduction allows the algorithm to focus efforts in the areas that are most likely to have a high quality (or possibly the global) optimum. Figure 5 shows a schematic representation of the optimization algorithm.

Figure 5

Diagramm of optimization algorithm. Diagrammatic representation of stochastic optimization algorithm.

It is demonstrated in the Results Section that this algorithm is highly efficient in finding the best solution for the den class="Chemical">sign of n class="Chemical">PCs.

Porous silicon PC fabrication

Porous silicon (P-Si) PCs were prepared by electrochemical anodization of highly boron-doped p+-type (100) crystalline silicon wafers with resistivity <0.005 Ω ⋅ cm in a hydrofluoric acid solution composed of ethanol, HF and glycerin in a volume ratio of 7:3:1 (70 ml of ethanol, 30 ml of HF and 10 ml of glycerol). We fabricated high and low refractive index layers by alternating the current density from 3.0 mA/cm2 to 40.0 mA/cm2 respectively. The corresponding refractive index values, needed for the design of the PC structure, were obtained by using the Bruggeman approximation as reported in[22]. Using etching rates of v = 14.49 nm/s and v = 1.72 nm/s we calculated the time at which each current density needed to be applied to form the desired thickness of the layers. Also, to avoid HF concentration decrease which causes a porosity gradient in the formed layers, we implemented 1 second long pauses to the etching time. After electrochemical etching, the samples were rinsed in ethanol for 10 minutes and dried under a nitrogen stream. We subsequently oxidized the samples for stabilization of the P-Si at 300 °C during 15 minutes.

The morphology of P-Si varies with porosity, and as shown in the SEM image 4 the high and low porosity layers (low and high refractive index respectively) fabricated in this study resemble a coral-like structure. The surfaces between layers are not smooth and present rugosity S. One way to determine S is to measure the maximal variations in the layer thicknesses which can be obtained from SEM imaging. Here we measured the layer thicknesses of our samples using a Hitachi S5500 electron microscope and averaged over the measurements finding S = 7.6 ± 0.1 nm.

Experimental reflectance spectra were measured using an UV-VIS-IR Spectrophotometer (Shimadzu UV1601) in three different places of each sample to asess homogeneity. The average of these measurements (R) are normalized with regard of a standard aluminium reflectance spectrum (R) using the relation: E = (R ⋅ R)/100. These normalized reflectance spectra E are presented in Figs 1 to 3 as the experimental reflectance spectra (E and E) of each fabricated PC mirror.