1. Introduction

Integrated photonics plays an indispensable role in optical communication, computing, and sensing. Historically, integrated optics uses opto-geometric parameters in order to optimize the performances of the devices, based-on of either analytical [1] or numerical models [2,3]. With the constant growth of data traffic between the consumers and remote data centers, the use optical transceivers is now well established, which will lead in one way or another to an increase of the photonics device density on chip. Thus, reducing the footprint of integrated optics will becomes a further challenge in the next decade, which brings to the vast concept of nanophotonics, taking place when the dimension of the photonic component, or a part of it, becomes smaller compared to the wavelength (<$\lambda /2n$) [4]. While nanophotonics devices exhibiting a single longitudinal or transversal period subwavelength pattern can be optimized with a parametric scanning, most of the optimization work in this domain rely on shape or topology optimization, also known as inverse design techniques [5–11]. Such methods, lead to efficient designs with non-intuitive shapes for which the manufacturing processes remains challenging. Furthermore, these methods are based on the computation of the gradient of the objective function. Consequently, they might be stuck in local minima/maxima rather than achieving the best design in the search parameter space. In this work, we concentrate on optimizing a silicon mode converter which has been identified as a critical component being realized in all the integrated photonic systems, while playing a fundamental role in wave division multiplexing and demultiplexer systems [12–17].

Here, we apply a global optimization method based on a statistical learning strategy in order to design a compact metamaterial mode converter structure at the telecommunication wavelength. Our optimization method is related to the class of Bayesian optimization methods and is known as Efficient Global Optimization (EGO) [18,19]. Contrary to the traditional common global optimization strategies like Genetic Algorithms (GAs) [20], EGO is not based on adaptive sampling, but on a surrogate model constructed on the basis of available objective function evaluations. This surrogate model utilizes statistical learning criteria related to the optimization target (usually called merit function) in order to identify which design (set of parameters) should be tested in the next iteration that would provide better results close to the predefined goal.

In general, the EGO is based on two phases. The first one is the Design Of Experiment (DOE), which aims at generating an initial database of designs. In essence, a uniform sampling strategy (e.g. the Latin Hypercube Sampling) [21] is deployed in order to generate different designs. Each of these designs is evaluated using an electromagnetic solver, and the value of a cost function is deduced from the numerical EM field. In the second phase, using the data obtained from the DOE, a Gaussian Process (GP) model is constructed to fit these data. This GP model allows us to predict the values of the cost function in the parameter space without the need to perform additional electromagnetic simulations. Once this GP model is determined, one can estimate at any point of the design space, the objective function (mean of the GP model) and an uncertainty value (variance of the GP model). The mean and the variance are used together to determine a statistical merit function. In our case, we rely on the expected improvement, which is a function whose maximum defines the next design parameters set to be evaluated. That is to say, in the search parameter space where this function is maximized, we extract the corresponding parameter values, and the corresponding design will be simulated using our electromagnetic solver. Then the database is updated accounting for this new observation (construction of a new GP model based on the updated database). We repeat this process until a predefined convergence criterion is reached, or when the expected improvement is sufficiently small. The reader can refer to Ref. [19] for more details about our EGO method where we provide an analytical example to show the performance of the method and compare it with a state of the art stochastic global optimization method. We emphasize that EGO is a global optimization method, which means that it will not be stuck in a local minima/maxima like the classical gradient methods [5,6,22,23]. Besides, in many cases, this method is computationally cheaper than classical global evolutionary strategies like GAs [20,24]. We also note that it has been recently shown that Bayesian based optimization techniques outperform most of the classical optimization methods [19,25].

Besides, it is worth discussing the main differences between EGO and the Artificial Neural Networks (ANNs) given the increasing popularity of ANN-based modeling approaches in the photonics research community. In fact, ANNs are not really designed for optimization, but for the construction of surrogate models, which can be then employed for different tasks. In general, the main difficulty with ANNs is related to the large amount of input data needed to train the network in order to be accurate for all possible configurations. In general, this limits strongly the number of parameters.

In contrast, the EGO algorithm is a metamodel-based optimization method, which focuses the search in promising regions for the minimizing a cost function. As a consequence, it is less sensitive to the dimension than ANNs. However, for EGO, there is anyway a limitation due to the fact that it relies on the construction of an internal Gaussian Process model, to decide which new point should be evaluated. This limitation is difficult to define precisely in general, because it depends on the studied problem. In the literature, it is reported that EGO should in general not be employed when dealing with more than one hundred parameters. Nevertheless, high dimensional Bayesian optimization is a very active research topic. In particular, preconditioning methods are studied, in order to reduce the number of active parameters effectively used during optimization.

Here we exploit a in-house fullwave 3D electromagnetic solver based on the Discontinuous Galerkin Time-Domain (DGTD) method from the DIOGENeS [26] software suite. This solver has been specifically designed for modelling and simulating nanoscale light-matter interaction problems [27] with complex geometrical shapes. Shortly, the DGTD method [28] can be seen as a mixture of a classical (continuous) Finite Element Time-Domain (FETD) and Finite Volume Time-Domain (FVTD) method. Consequently, the DGTD method combines attractive features of FETD and FVTD methods. In particular, a DGTD method relies on a high order approximation of the unknown field with a polynomial interpolation method ($\mathbb {P}_p$ where $p$ is the interpolation degree). It is formulated on a fully unstructured tetrahedral mesh allowing for a local refinement of the mesh and the possibility of using curvilinear cells for a high order approximation of sophisticated geometry features with strong discontinuity. These modelling features of the DGTD makes it an ideal electromagnetic solver for nanophotonic devices. For more details about our solver, the reader can refer to Refs. [26,27].

2. Design and results

In Fig. 1, we present a 3D schematic view of the structure under consideration. We consider a graded index configuration acting as a mode converter from a multimode waveguide of width $W$ excited on its fundamental mode towards a single-mode waveguide of width $w$ at wavelength $\lambda =1310~$nm. In the end, one can consider such a device as a metalens, mimicking the behavior of a lens at the nanophotonics scale, and applied to an optical mode. This metalens consists of several Si strips of varying widths with the same length $L$ along the propagation direction $z$. The height of the strips is fixed as $h=310$ nm (in $y$ direction). The Si regions (green color) is placed on top of a semi-infinite substrate made of $SiO_{2}$ (red color) and is coated from above by a semi-infinite air region (not shown here). The structure is terminated from all sides by Perfectly Matched Layers (PMLs) (not shown). The permittivities of the Si and SiO$_{2}$ regions at the operating wavelength are $3.506^{2}$ and $1.4468^{2}$, respectively. The principal aim of this design is to convert the fundamental mode injected from the surface $S_{in}$ (solution to the 2D waveguide problem with Si width $W=3000~$nm) to a single-mode waveguide at surface $S_{out}$ (solution to the 2D waveguide problem with Si width $w=300~$nm) with maximum efficiency. Certainly, in order to accomplish this target, one needs to optimize the widths, and the length $L$ of the Si strips simultaneously. In other words, we search for the optimal configuration to maximize the following Figure Of Merit (FOM)

 

Fig. 1. Schematic view of the 3D photonic metalens. The structure consists of Si region (green part) on top of subtrate made of SiO$_{2}$ (red part). The Si region is divided into three parts; the input port with width $W=3000~$nm, the output port with width $w=300~$nm, and in-between, several Si strips with lenght L and height $h=310~$nm. The input mode is injected from the surface $S_{in}$ and the FOM (see Eq. (1)) is computed at the output surface $S_{out}$. The widths of the strips are denoted by $e_{i}$, $i \in \{0,6\}$.

Download Full Size | PDF

In Eq. (1), $E_{in}$ and $H_{in}$ represent the input mode, which is the solution of a 2D waveguide of width $W=3000~$nm at the input surface S$_{in}$. Similarly, $E_{out}$ and $H_{out}$ represent the solution of a 2D waveguide with $w=300~$nm at the output plane S$_{out}$, while $E$ and $H$ are the solution we obtain after the propagation through the strips. The FOM in Eq. (1) combines the maximization of the total transmission and the maximization of the overlap between the solution $E$, $H$ and the given output solution $E_{out}$, $H_{out}$ at the output surface $S_{out}$.

In order to maximize the FOM provided in Eq. (1), we utilize the EGO mutually with our DGTD solver. In the optimization problem we consider $8$ parameters that represent the length of the strips $L$ and their widths $e_{i}$, $i \in \{0,6\}$ as shown in Fig. 1. We respect the fabrication constraints in the search parameter space, such that the widths of the strips are allowed to vary from $50$ nm to $200$ nm, and the period is fixed to be $250$ nm. Besides, the length $L$ is ranging from $1\mu$m to $4\mu$m in order to realize a compact mode converter device. The optimization results are shown in Fig. 2. In Fig. 2(a), the blue points represent the first phase of the optimization, i.e., the design of experiments (DOEs); here we examined only 50 designs. As we have outlined above, based on this initial database, a surrogate model is then constructed, and the second optimization phase (represented by the black points in Fig. 2(a)) is performed. During the second phase, the current point is compared to the previous point, and the best point is kept (see green curve). As it can be noticed, the optimizer obtains several interesting points during the first ten iterations of the second phase (first ten black points). After that the green curve remains constant meaning that no further improvement can be anticipated in this region. This illustration can also be inferred in Fig. 2(b) where the merit function is plotted as a function of the iterations of the second phase (starting from 50). One notices that the merit function remains close to zero after iteration 85. This is an indication that no major improvement is expected and that the optimization process can be resumed.

 

Fig. 2. Optimization results: (a): represents the 1 - FOM (to maximize the FOM given by Eq. (1)) as a function of the solver calls. The blue points represent the DOE (50 points), the black ones represent the optimization iterations, and the green curve indicates the best results obtained along the optimization process. (b): the merit function versus the number of solver calls (for the optimization iterations, i.e, after 50 iterations as indicated by the green curve in (a)).

Download Full Size | PDF

As it can be seen in Fig. 2(a), the best design obtained provides nearly $80.2\%$ of efficiency. Yet, the two preceding points yield almost the same value of the objective function. In Table 1, we extract the parameters of the two interesting designs (provideing approximately $80\%$ of efficiency). We observe from the two first rows in Table 1 that the widths of the central strips are approximately the same for the two designs. Nevertheless, the two designs exhibit different values of $L$. This suggests that the EGO obtained at least two global optimal designs in the search parameter space, which give nearly the same value of the objective function. This feature has already been illustrated in our previous work [19], using an analytical example, where we have shown that the EGO is able to capture effectively almost all the global minima/maxima in the search parameter space. It is worth stating that during the optimization process, a high order polynomial order $\mathbb {P}_{4}$ is used in our DGTD solver to approximate the unknown field, together with a coarse mesh (around 45 000 cells), which is enough to obtain accurate results and obtain the optimization results in a reasonable time (typically with 48 cores, one fullwave solver requires almost 40 minutes). However, more precise results are presented in Table 2 for the best optimized design (first row in Table 1) as a function of the mesh size and polynomial order. In this table, we note that the efficiency of the device is approximately $80\%$. The corresponding field map in $x-z$ plane is given in Fig. 3(a), where the focusing is clearly demonstrated. For the second design (the second row in Table 1), and using a refined mesh with $177~000$ cells, we obtain $78.5\%$ of efficiency, which is slightly smaller than the best design that was obtained previously using the coarse mesh (see Table 3). This illustrates the accuracy of our optimization results.

 

Fig. 3. (a): Field map of $\Re e(H_{y})$ for the optimized design in the $x-z$ plane obtained from our DGTD solver using $\mathbb {P}_{4}$ interpolation and a mesh with $334~000$ cells. The parameters are given in the first row in Table 1. (b): numerical convergence assessment for the TE case using the values provided in Table 3 as a function of the length of the strips L using different polynomial orders, and different mesh sizes. As it can be seen from the shaded region, there are three global points similar to what we have found in Table 1 and the best FOM is obtained at $L=3700$ nm.

Download Full Size | PDF

Table 1. Optimized deigns obtained from the EGO (L and $e_i$ in nm). The results show that there are two global points in our design with nearly the same FOM $\approx 80\%$., however, the parameters are slightly different. Here, we considered PML thickness of $500$ nm, which is proved to be enough.

View Table | View all tables in this article

Table 2. The convergence of the optimized design using EGO (first row in Table 1) as a function of the mesh size and the polynomial order in the DGTD solver for the total transmission T (third column) and for the FOM (last column).

View Table | View all tables in this article

Table 3. Values of the strip widths to obtain a gradient index lens profile for TE polarization ($e_i$ in nm). These widths are determined by computing the effective mode index for each strip [14], see the text for more details.

View Table | View all tables in this article

In order to render further physical insights on the results achieved from the optimization method, we provide an estimation of the widths of the strips using the graded index lens profile [14]. Following the same strategy presented in our former work for the TM case [14], the widths of the strips are roughly estimated by calculating the effective mode index for each strip, by considering a single strip placed in a periodic configuration, with period $250$ nm. Changing the width of this single strip and computing the corresponding effective mode index creates a library for effective mode indices as a function of the corresponding widths. Afterwards, only strips that provide a graded index lens profile [14] will be chosen from the library. In Table 3, we present the values based on this approximated technique [14], which are similar to the best design obtained from our optimization approach (see the first row in Table 2). Nevertheless, the graded index profile approach does not give any information about the optimal length along the propagation direction. Therefore, one has to sweep over the length L in order to obtain the optimal length and maximize the FOM depicted in Eq. (1). Using the values given in Table 3, we sweep over the length L to obtain the maximum FOM. Interestingly, we found that there are at least two values of the length L that yield the same value of the FOM around $L=3700$ nm and $L=3400$ nm (see the highlighted region in Fig. 3(b)). This result validates the outcome from the optimization method (see Table 1) and infers that we rigorously obtained the best parameters that provide the maximum FOM as described in Eq. (1) for the mode converter metamaterial depicted in Fig. 1 with fixed length L. We would like to emphasize that due to the feature of the TE mode (electric field polarized perpendicular to the strips), the near field coupling must be taken into account in the calculations. Furthermore, most of the electromagnetic energy is located inside the small gaps between the strips (the smallest gap between strips is in the centre, and it is approximately $50.5$ nm) as it is indicated in Fi. 3(a). Consequently, a high order polynomial order is required to capture the strong discontinuity between the strips, as it can be seen in Fig. 3(a). Accordingly, numerical methods that rely on low order approximations like Finite Difference Time-Domain (FDTD) might encounter serious problems in achieving the convergence, as it is indicated in Fig. 4(a), where the FDTD method does not converge even with a uniform mesh size as small as 10 nm (approximately 30 million cells) for the TE case (we refer to Fig. 3(b) for the demonstration of the numerical convergence of the DGTD method using different mesh sizes, different mesh types, and different polynomial orders). It is worth mentioning that for the TM case (parameters are given in Ref. [14]), where the electric field is polarized along the same direction of the strips, the results obtained with the FDTD method globally agree with the ones acquired from the DGTD solver (the maximum FOM is almost around $L=2000$ nm as it is shown in Fig. 4(b)). The DGTD method outperforms the FDTD method in particular because we use a fourth order interpolation of the electromagnetic field in each mesh cell to achieve a high accuracy using a relatively coarse mesh. In Fig. 3(b) we can see the global validation of DGTD solver on one hand, and the importance of the high accuracy even in the simple TM case on the other hand.

 

Fig. 4. FDTD convergence study (a): TE case where the parameters are identical to the ones used in Fig. 3(b) for the DGTD solver. (b): TM case considered in Ref. [14], where all the opto-geometrical parameters are identical to the ones used in Ref. [14].

Download Full Size | PDF

Having demonstrated the capacity of our optimization method together with the importance of our DGTD solver in obtaining a converged solution for the TE case, we go one step further by investigating the influence of varying the lengths of the strips on the FOM. We start from the optimized widths given in the first row in Table 1, while changing the length of each individual strip and keeping the lengths of the three central ones as in Table 1 (see Fig. 5). In this case, we have only five optimization parameters, which correspond to the lengths of the five outer strips (symmetry properties are applied). Interestingly, using our EGO method, we found that changing the length of the strips yields an enhanced efficiency of the device as high as $87\%$ as indicated in Table 4 and Fig. 5. As it may be noticed, we have a gradient change in the length of the strips from the centre to the edges, which ensures the change of the phase and enhancing the focusing efficiency as it is demonstrated in Fig. 5.

 

Fig. 5. (a): 3D mesh for the optimized design with different lengths for the strips. (b): Field map of $\Re e(H_{y})$ for the optimized design with changing length in the $x-z$ plane obtained from our DGTD solver using $\mathbb {P}_{4}$ interpolation and mesh with $334~000$ cells. The strips width parameters are given in the first row in Table 1.

Download Full Size | PDF

Table 4. Optimized values of the strip lengths obtained from the EGO method (L$_i$ in nm) with1 FOM$=87\%$. The values of the strip widths are similar to the ones in the first row in Table 1. Here we fix the the three length of the center strips (value shown in first row in Table 1) and change the lengths of the five outer strips (symmetry is applied).

View Table | View all tables in this article

3. Conclusion

By combining a global optimization method and a high order fullwave solver, we achieved optimal designs for a 3D photonic graded index mode converter. We tackled the most demanding situation where the strips are perpendicular to the incident electric field polarization. In this case, the strong near field coupling between the strips must be considered, and the small gaps between the strips need to be optimized in order to provide high efficient designs. Our 3D fullwave solver based on the DGTD method is an ideal solver to tackle this problem, due to the possibility to use a high order polynomial interpolation for the numerical field approximation, and to deal with structures with strong variatiobn sof space scales. We considered two cases. The first one when all the Si strips are sharing the same length, but with different widths. In this classical case we have shown that our statistical learning-based optimization technique EGO is able to capture nearly all the global designs in our parameter space (all the designs that provide the maximum efficiency, even if the parameters are different) and the efficiency can reach 80%. Secondly, based on the widths of the strips obtained, we investigated the impact of varying the length of the strips on the focusing efficiency. We have shown that maximum efficiency of 87% can be obtained using fully CMOS compatible and easily fabricated design.