Fast numerical methods for the design of layered photonic structures with rough interfaces

: Modiﬁed boundary conditions (MBC) and a multilayer approach (MA) are proposed as fast and efﬁcient numerical methods for the design of 1D photonic structures with rough interfaces. These methods are applicable for the structures, composed of materials with arbitrary permittivity tensor. MBC and MA 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 ﬁeld 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.


Introduction
Photonic crystals are of high interest for many practical applications [1,2,3]. One-dimensional photonic crystals are known for a long time and have been studied thoroughly [4,5,1,6]. A lot of effort was also put into the optimization of finite multilayer films [7,8,9,10]. Different optimization procedures, such as the quasi-Newton algorithm [8] 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 antireflection in more than one frequency band. Reentry [11] 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. In order to avoid time-consuming 3D simulations, we propose a multilayer approach (MA) or modified boundary conditions (MBC) to take roughness into account. These approaches are 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 [12,13,14].
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 Section 4 we describe MA and MBC for rough interfaces in detail. In Section 5 we provide numerical results of MBC and MA for rough interfaces. A sinusoidal and a random hilly profiles of roughness are taken as examples. Results for isotropic and anisotropic uniaxial crystals are shown.

Optimal design of a dual band omnidirectional mirror
In most cases it is not possible to find an optimal design of a structure analytically. Various optimizers [15,16] can be helpful to find good solutions numerically.
If a layered structure consists of N layers, one has (at least) N real-valued parameters (thickness of each layer) to be optimized. The search space should be defined for each parameter within a reasonable interval [d min , d max ]. 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 high, but also when the fitness function has a complicated behavior, e.g. when it has many local optima. According to our experience [17] 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.
To illustrate the optimization procedure, we optimize the omnidirectional reflection of the structure in two bands. Hereafter, we denote the reflection coefficients: R ss , R sp , R ps , and R pp . The first and second indices correspond to polarization of the incident and reflected light respectively. The s-and p-polarized light correspond to electric or magnetic field parallel to layers respectively. The multilayered structure 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, we consider flat interfaces between the layers. In case of rough interfaces MBC or MA can be used (see Section 3 and Section 4). MBC and MA 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, we must define the fitness function, i.e. the goal of the optimization. We optimize the omnidirectional reflection of the unpolarized light in two bands B 1 := 800 − 1300 meV and B 2 := 2000 − 2500 meV. Therefore, we define the fitness function as the averaged reflection: where the reflection coefficients are integrated over the angles of incidence θ ∈ [0, π 2 ] and energy bands B 1 and B 2 . The thickness of the i-th layer (i = 1, ..., N) is d i . We now need to maximize the fitness function Eq. (1) of N variables.  The results of ES the optimization are shown in Fig. 1. One may see that a high reflectivity is achieved in the B 1 and B 2 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 multilayered 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 roughness into account at least approximately. We propose two ways for doing this: 1) multilayer approach (MA) 2) modified boundary conditions (MBC). These methods work very fast, this is essential for the optimization problems.

Model
Geometry. The rough boundary is located between two subspaces: z < −δ , ε = ε 1 , µ = 1, and z > +δ , ε = ε 2 , µ = 1. The region |z| < δ is occupied with the rough boundary, see Fig. 2. Both sizes of the roughness 2δ and l are supposed to be small in comparison with the wavelength. The volume fractions of two components are Rough boundary between interfaces with permittivities ε 1 and ε 2 .
We can divide the region −δ < z < δ into M layers and assign an effective permittivity ε i eff (i=1,...,M ) to each of these layers. However, standard formulas for effective permittivity, such as Maxwell-Garnett or Bruggeman, are not valid in the general case of arbitrary anisotropic materials and geometric profile of roughness. Therefore, we can only use the bounds for effective permittivity: These bounds for effective permittivity were first derived by Wiener [18] and they 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:

Model
Geometry. The rough boundary is located between two subspaces: z < −δ , ε = ε 1 , µ = 1, and z > +δ , ε = ε 2 , µ = 1. The region |z| < δ is occupied with the rough boundary, see Fig. 2. We suppose that the dielectric permittivity is a function of z in this region. Thus, The wave is incident from z = −∞ and propagates in the XZ plane, thus (∂ /∂ y = 0). The angle of light incidence is θ with respect to Z axis.
Let us choose the rectangular l × 2δ contour at the rough boundary. Then the Maxwell's equations for the circulation around this contours are: We can write the left part of the first equation as: because B n is continuous at the boundary and can be considered as constant there. Thus we obtain: The right part of the first equation (5) can be written as: Assuming δ λ and l λ , we obtain Here all fields are written at z = 0. The first and second equations originate from the first equation (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 Eq. (6) can be simplified, if we consider a wave of a certain polarization. For E-waves we have E = (0, E y , 0) H = (H x , 0, H z ).
The second and third equations are satisfied identically, while the first and the last equations reduce to Similarly, we have for H-waves and we obtain Here The generalization of the boundary conditions Eq. (7) and Eq. (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 boundary conditions Eq. (7) and Eq. .
Here the superscripts (1) and (2) correspond to the first and second bordered material, respectively. 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. [19]. The wave equation for the E-wave in this region is: The integration of this equation over the δ -region at the boundary (see Ref. [19] for the details of calculation) leads to the following boundary conditions for the electric fields E 1,2 and their derivatives E 1,2 = ∂ E 1,2 /∂ z at each side of the boundary (we omit the subscript y for reasons of simplicity): Here α and β are the parameters, which are independent of the wavelength and incidence angle. The wave equation for the H-wave in δ region is If ε 33 is independent of x, then Eq. (11) can be written as: Its integration leads to the following BC for H y field: Expansion similar to Eq. (6), but to the second order on δ /λ , leads to α = t 1 , β = (ε 2 − ε 1 )δ 2 /2 and γ = t 2 . Here the parameters α and β are the same as in Eq. (10).

Numerical results of modified boundary conditions and multilayer approach
We demonstrate the application of MBC and MA for rough interfaces by the following examples.

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 boundary between air and 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: sin( πy d y ) Period and height of the grating are much smaller than the considered wavelength. As a consequence, the dependence of the reflectivity on the azimuthal angle ϕ is almost negligible, despite the absence of the symmetry with respect to ϕ.
We calculate the reflection spectra of the structure in four different ways. The results are shown in Fig. 4.a and Fig. 4.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. One can see 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], β [cm 2 ] and γ [cm].
In Fig. 4 one can see, that MA is in very good agreement with the reference solution for both polarizations and for all angles of incidence θ . MBC also is 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 ). While 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 show a very good agreement. Thus the application of MA and MBC to the considered sinusoidal roughness profile is justified.

Random hilly roughness
In order to demonstrate that MA and MBC methods are not restricted to sinusoidal approximation of 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 h max = 10 nm. The roughness profile is periodic, with periods d x = d y = 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 in Sec. 5.1 for the sinusoidal roughness, we consider separately isotropic and anisotropic materials. The results of MA and MBC for an isotropic material with ε = 5 are shown in ig. 5. Reflectivity of the structure shown in Fig. 3.b. (a) and (b) correspond to an isotropic material with ε = 5, (c) and (d) correspond to an anisotropic material with ε 11 = ε 22 = 7, ε 33 = 9. Red curve -reference solution (CST Studio), black dash-dotted curve -flat slab with thickness 105 nm, blue curve -multilayer approach (MA), green curve -MBC between air and slab.   Fig. 3.c for normally incident light. Permittivities of materials are ε 1 = 2 and ε 2 = 7 (on top). Red curve -reference solution (CST Studio), black dash-dotted curve -flat interfaces, blue curve -multilayer approach (MA) at interfaces, green curve -MBC at interfaces.

Conclusions
We presented modified boundary conditions (MBC) and a multilayer approach (MA) 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 scattering matrix technique as field solver and an evolutionary strategy as numerical optimizer these methods reveal a fast way of optimizing 1D photonic structures for various applications.