A multilayer approach (MA) and modified boundary conditions (MBC) are proposed as fast and efficient numerical methods for the design of 1D photonic structures with rough interfaces. These methods are applicable for the structures, composed of materials with an arbitrary permittivity tensor. MA and MBC are numerically validated on different types of interface roughness and permittivities of the constituent materials. The proposed methods can be combined with the 4x4 scattering matrix method as a field solver and an evolutionary strategy as an optimizer. The resulted optimization procedure is fast, accurate, numerically stable and can be used to design structures for various applications.
© 2011 Optical Society of America
Photonic crystals are of high interest for many practical applications [1–3]. One-dimensional photonic crystals are known for a long time and have been studied thoroughly [1,4–6]. A lot of effort was also put into the optimization of finite multilayer films [7–10]. Different optimization procedures, such as the quasi-Newton algorithm  or genetic algorithms [9, 10] have been applied to improve the performance of finite multilayer mirrors.
Some applications require the optimization of the structure for high reflection or anti-reflection in more than one frequency band. Reentry  is an example of such an application, where reflection in two ranges is desired.
In practice, roughness of the interfaces is inevitable and causes severe numerical problems, because it destroys the symmetry of the structure. Regarding the problem of interface roughness in 2D photonic crystals, we address the reader to . Optical response from rough multilayers has been studied in a number of works. For example, perturbation theory and scalar theory of diffraction have been applied in [13, 14]. We propose two alternative methods to take into account the roughness. Namely, a multilayer approach (MA) and modified boundary conditions (MBC). The main advantages of these methods are 1) MA and MBC are not restricted to smooth roughness profiles 2) both methods are also applicable to anisotropic materials 3) MA and MBC can be easily combined with the 4x4 scattering matrix method as field solver. The resulting procedure is fast, numerically stable, and easy to implement. Details of the scattering matrix formalism can be found in [15–17]. The drawback of MA and MBC is that the scattering losses are neglected. However, this is essentially justified by the assumption that the size of roughness is small compared to wavelength.
This paper is organized as follows. In Section 2, we illustrate the optimization procedure on a dual band mirror with flat interfaces. In Section 3 and 4, we describe MA and MBC for rough interfaces in detail. In Section 5, we provide numerical results of MA and MBC for rough interfaces. A sinusoidal and a random hilly profile of roughness are taken as examples. Results for isotropic and anisotropic uniaxial crystals are shown.
2. Optimal design of a dual band omnidirectional mirror
If a layered structure consists of N layers, one has (at least) N real-valued parameters (thickness of each layer) that need to be optimized. The search space should be defined for each parameter within a reasonable interval [dmin, dmax]. For such optimizations, various algorithms, such as genetic algorithms (GA), micro-genetic algorithms (MGA), or evolutionary strategies (ES) may be applied. The optimization task becomes difficult not only when N is large, but also when the fitness function has a complicated behavior, e.g. when it has many local optima. According to the experience , ES is very powerful in real parameter optimization problems and outperforms genetic algorithm (GA), particle swarm optimization (PSO) etc. in most cases. Therefore, we applied it to the multilayer problem without extensive testing. We used an (m+n) ES with adaptive mutation for the optimization in the following example. Here m is the initial number of parents, n is the number of children created in each generation. In the following example we used m = 5 and n = 40.
In order to illustrate the optimization procedure, the omnidirectional reflection of the structure is optimized in two bands. Hereafter, the reflection coefficients are denoted: Rss, Rsp, Rps, and Rpp. The first and second indices correspond to the polarization of the incident and reflected light, respectively. The s- and p-polarized light corresponds to an electric or magnetic field parallel to the layers, respectively. The multilayered structure that is to be optimized consists of N=14 layers, embedded between air interfaces. It is composed of two alternating materials with permittivity tensors ɛ̂1 and ɛ̂2:
We take these permittivities just as examples. Any other materials with full permittivity tensors could be used instead. In this example, flat interfaces between the layers are considered. In the case of rough interfaces, MA or MBC can be used (see Section 3 and 4). MA and MBC are applicable when the characteristic size of the roughness is much smaller than the wavelength in the material (δ << λ).
As a necessary step for the optimization, the fitness function must be defined, i.e. the goal of the optimization. We optimize the omnidirectional reflection of the unpolarized light in two bands B1 : = 800 – 1300 meV and B2 : = 2000 – 2500 meV. Therefore, we define the fitness function as the averaged reflection:Eq. (1) of N variables.
The results of ES the optimization are shown in Fig. 1. It can be seen that a high reflectivity is achieved in the B1 and B2 bands. Figure 1(b) shows the geometric profile of the obtained optimal structure.
In the optimization procedure, we limited the number of fitness calls to 14200. As a consequence it is not guaranteed that the global optimum is found. As long as a sufficiently good solution is found, this is acceptable in practice. The required optimization time was several hours in Matlab on Intel(R) Core(TM) i7 CPU 2.67GHz with 8 GB of RAM.
The interfaces of fabricated multilayer films are not ideally flat. Roughness between interfaces can be a result of fabrication inaccuracies or interdiffusion, if the structure operates at high temperatures. Therefore, it is desired to take into account the roughness - at least approximately. We propose two ways for doing this: 1) multilayer approach (MA) 2) modified boundary conditions (MBC). These methods work very fast, which is essential for the optimization problems.
3. Multilayer approach for rough interfaces
The rough interface is located between two subspaces: z < −δ, ɛ= ɛ1, μ = 1, and z > +δ, ɛ= ɛ2, μ = 1. The region |z| < δ is occupied with the rough interface, see Fig. 2. Both sizes of the roughness 2δ and l are supposed to be small compared to the wavelength. The volume fractions of two components are f1 = V1/V and f2 = V2/V, where V is the total volume of the rough region.
We can divide the region −δ < z < δ into M layers and assign an effective permittivity to each of these layers. However, standard formulas for the effective permittivity, such as Maxwell-Garnett or Bruggeman, are not valid in the general case of arbitrary anisotropic materials and geometric profiles of roughness. Therefore, we can only use the bounds for the effective permittivity:21] and are valid for arbitrary anisotropic materials and shapes of the inclusions.
As an effective permittivity for each layer, the mean value of the two bounds can be taken:
4. Modified boundary conditions for rough interfaces
The rough interface is located between two subspaces: z < −δ, ɛ= ɛ1, μ = 1, and z > +δ, ɛ= ɛ2, μ = 1. The region |z| < δ is occupied with the rough interface, see Fig. 2. We suppose that the dielectric permittivity is a function of z in this region. Thus,
Let us choose the rectangular l × 2δ contour at the rough interface; then the Maxwell’s equations for the circulation around this contour are:Eq. (7), Eq. (8), Eq. (10), and Eq. (12), are valid in the both, CGS and SI, unit forms.
We can write the left part of the first equation (5) as:5) can be written as: 5) being written for the different components of the magnetic field. The third and forth equations are derived from the second equation of (5).
The boundary conditions in Eq. (6) can be simplified, if we consider a wave of a certain polarization. For E-waves, we have
The generalization of the boundary conditions in Eq. (7) and (8) for the anisotropic structures leads to very lengthy expressions. However, they are essentially simplified for uniaxial crystals, where the dielectric permittivity tensor is diagonal and ɛ11 = ɛ22 ≠ ɛ33. Then the boundary conditions Eq. (7) and (8) hold with , and
The second order correction on δ/λ needs an additional expansion of the fields in the integrands of Eq. (5). This can be done also by integration of the wave equations over the interface region |z| < δ, see Ref. . The wave equation for the E-wave in this region is:22] for the details of calculation) leads to the following boundary conditions for the electric fields E1,2 and their derivatives E′1,2 = ∂E1,2/∂z at each side of the boundary (we omit the subscript y for reasons of simplicity):
The wave equation for the H-wave in the δ region isEq. (11) can be written as: Eq. (6), but to the second order on δ/λ, leads to α = t1, β = (ɛ2 − ɛ1)δ2/2 and γ = t2. Here the parameters α and β are the same as in Eq. (10).
5. Numerical results of multilayer approach and modified boundary conditions
We demonstrate the application of MA and MBC for rough interfaces by the following examples.
5.1. Sinusoidal roughness profile
Consider the reflection from a slab with a sinusoidal grating on the top, see Fig. 3(a). The sinusoidal grating emulates a rough interface between air and the slab. The permittivities of the slab and of the sinusoidal grating are the same ɛ = 5. The structure is embedded in air. The dielectric permittivity in the region −h/2 < z < h/2 is:
We calculate the reflection spectra of the structure in four different ways. The results are shown in Fig. 4(a) and (b). The red curve represents the solution obtained with CST Studio (Frequency Domain Solver). We refer to it as the “reference solution”. The black dash-dotted curve represents the reflection from a flat slab with thickness L + h/2 = 105 nm. It can be seen that the reflection from a flat interface deviates significantly from the “reference solution”. This is especially pronounced for higher energies. The blue curve corresponds to the multilayer approach, with M = 6 layers in the sinusoidal region −h/2 < z < h/2. The green curve corresponds to the MBC solution (Eq. (10) and Eq. (12)), with optimal parameters α [cm], β [cm2] and γ [cm].
In Fig. 4 it can be seen, that MA is in very good agreement with the reference solution for both polarizations and for all angles of incidence θ. MBC is also in good agreement. The parameters α, β and γ of the MBC are independent of the frequency and incidence angle. These parameters can be found either by fitting to the reference solution, as done in this paper, or by fitting to experimental reflection/transmission spectra of multilayers.
Now we consider the reflection from the same structure, but composed of anisotropic material with ɛ11 = ɛ22 = 7, ɛ33 = 9. It should be noted that the application of MBC is very much simplified in the case of isotropic, or uniaxial crystals (ɛ11 = ɛ22 ≠ ɛ33), whereas MA is equally easy to apply for arbitrary permittivity tensors.
The results of MA and MBC are shown in Fig. 4(c) and Fig. 4(d). It can be seen that MA (blue curve) is in very good agreement with the reference solution (red curve). MBC also shows a very good agreement. Thus the application of MA and MBC to the considered sinusoidal roughness profile is justified.
5.2. Random hilly roughness
In order to demonstrate that MA and MBC methods are not restricted to a sinusoidal approximation of a rough interface, we consider the reflection from the slab with random hills (polygons) on top, see Fig. 3(b). Maximum height of the hills does not exceed hmax = 10 nm. The roughness profile is periodic, with periods dx = dy = 50 nm. The characteristic size of the roughness is much smaller than the wavelength, therefore the dependence of the reflectivity on the azimuthal angle φ is practically negligible.
As described in Sec. 5.1 for the sinusoidal roughness, isotropic and anisotropic materials are considered separately. The results of MA and MBC for an isotropic material with ɛ = 5 are shown in Fig. 5(a) and (b). The results for an uniaxial crystal with ɛ11 = ɛ22 = 7, ɛ33 = 9 are shown in Fig. 5(c) and (d). Evidently, a good agreement of MA and MBC with the reference solution is observed.
If a structure consists of N layers with rough interfaces, then either MA or MBC can be applied to each of the N-1 interfaces. This is demonstrated in the next example.
5.3. Multilayers with rough interfaces
In this final example, we apply MA and MBC to a multilayered structure with rough interfaces. The structure is shown in Fig. 3(c). Structure is composed of two materials with permittivities ɛ1 = 2 and ɛ2 = 7, whereby the high index material is on top. Interfaces 1–3 are rough, with random hilly roughness profiles, as depicted in Fig. 3(b). The roughness of each individual interface was generated randomly. The maximum height of the hills does not exceed hmax = 10 nm. The structure is periodic with periods dx = dy = 50 nm.
The reflection spectra of a normally incident light are shown in Fig. 6. The red curve corresponds to the “reference solution”, obtained with the CST Studio. Some resonances are observed for energies above 9 eV. The reason is that the size of the roughness becomes comparable with the wavelength at these energies. The black dash-dotted curve corresponds to flat interfaces. The blue curve corresponds to MA with M = 6 layers in the interface regions. The green curve on the plot corresponds to MBC at each rough interface. It can be seen that MA and MBC in general demonstrate a better agreement with “reference solution”. This is especially pronounced in two regions, which are marked with arrows in Fig. 6.
It should be noted that MA and MBC combined with scattering matrix method work very fast. Calculation of one frequency point of the structure shown in Fig. 3(c) requires only about 1 ms with the scattering matrix code and about 1 min with the CST Studio. Calculations were performed on Intel(R) Core(TM) i7 CPU 2.67GHz with 8 GB of RAM.
We presented a multilayer approach (MA) and modified boundary conditions (MBC) as fast and efficient numerical methods for the design of layered structures with rough interfaces. These methods are applicable for arbitrary anisotropic materials. MA and MBC provide a good agreement with the reference solutions for sinusoidal and random hilly roughness profiles. Combined with the S-matrix method as a field solver and an evolutionary strategy as numerical optimizer these methods reveal a fast way of optimizing 1D photonic crystals for various applications.
References and links
1. J. Joannopoulos, S. Johnson, J. Winn, and R. Meade, Photonic crystals: molding the flow of light (Princeton Univ Pr, 2008).
2. A. Mekis, J. Chen, I. Kurland, S. Fan, P. Villeneuve, and J. Joannopoulos, “High transmission through sharp bends in photonic crystal waveguides,” Phys. Rev. Lett. 77, 3787–3790 (1996). [CrossRef] [PubMed]
3. V. Shklover, L. Braginsky, G. Witz, M. Mishrikey, and C. Hafner, “High-Temperature Photonic Structures. Thermal Barrier Coatings, Infrared Sources and Other Applications,” J. Comput. Theoret. Nanosci. 5, 862 (2008).
4. D. Berreman, “Optics in stratified and anisotropic media: 4 × 4-matrix formulation,” J. Opt. Soc. Am. 62, 502–510 (1972). [CrossRef]
5. P. Yeh, Optical waves in layered media (WileyNew York, 1988).
6. I. Abdulhalim, “Analytic propagation matrix method for linear optics of arbitrary biaxial layered media,” J. Opt. A: Pure Appl. Opt. 1, 646 (1999). [CrossRef]
7. F. Kärtner, N. Matuschek, T. Schibli, U. Keller, H. Haus, C. Heine, R. Morf, V. Scheuer, M. Tilsch, and T. Tschudi, “Design and fabrication of double-chirped mirrors,” Opt. Lett. 22, 831–833 (1997). [CrossRef] [PubMed]
8. T. Yonte, J. Monzón, A. Felipe, and L. Sánchez-Soto, “Optimizing omnidirectional reflection by multilayer mirrors,” J. Opt. A: Pure Appl. Opt. 6, 127–131 (2004). [CrossRef]
9. M. del Rio and G. Pareschi, “Global optimization and reflectivity data fitting for x-ray multilayer mirrors by means of genetic algorithms,” in “Proceedings of SPIE ,”, vol. 4145 (2001), vol. 4145, p. 88. [CrossRef]
11. D. Bose, E. McCorkle, C. Thompson, D. Bogdanoff, D. Prabhu, G. Allen, and J. Grinstead, “Analysis and model validation of shock layer radiation in air,” VKI LS Course on hypersonic entry and cruise vehicles, Palo Alto, California, USA (2008).
13. P. Bousquet, F. Flory, and P. Roche, “Scattering from multilayer thin films: theory and experiment,” J. Opt. Soc. A 71, 1115–1123 (1981). [CrossRef]
14. I. Ohlídal, “Approximate formulas for the reflectance, transmittance, and scattering losses of nonabsorbing multilayer systems with randomly rough boundaries,” J. Opt. Soc. Am. A 10, 158–171 (1993). [CrossRef]
15. L. Li, “Formulation and comparison of two recursive matrix algorithms for modeling layered diffraction gratings,” J. Opt. Soc. Am. 13, 1024–1035 (1996). [CrossRef]
16. D. Whittaker and I. Culshaw, “Scattering-matrix treatment of patterned multilayer photonic structures,” Phys. Rev. B 60, 2610–2618 (1999). [CrossRef]
17. S. Tikhodeev, A. Yablonskii, E. Muljarov, N. Gippius, and T. Ishihara, “Quasiguided modes and optical properties of photonic crystal slabs,” Phys. Rev. B 66, 45102 (2002). [CrossRef]
18. A. Fallahi, M. Mishrikey, C. Hafner, and R. Vahldieck, “Efficient procedures for the optimization of frequency selective surfaces,” IEEE Trans. Ant. Prop. 56 (2008). [CrossRef]
19. C. Hafner, C. Xudong, J. Smajic, and R. Vahldieck, “Efficient procedures for the optimization of defects in photonic crystal structures,” J. Opt. Soc. Am. A 24, 1177–1188 (2007). [CrossRef]
20. J. Fröhlich, “Evolutionary optimization for computational electromagnetics,” Ph.D. thesis, ETH Zurich, IFH Laboratory (1997).
21. O. Wiener, “Die Théorie des Mischkörpers für das Feld des stationären Strömung. Abh. Math,” Physichen Klasse Königl. Säcsh. Gesel. Wissen 32, 509–604 (1912). [PubMed]
22. M. Mishrikey, L. Braginsky, and C. Hafner, “Light propagation in multilayered photonic structures,” J. Comput. Theor. Nanosci. 7, 1623–1630 (2010). [CrossRef]