Angle-Based Parametrization with Evolutionary Optimization for OESCL-Band Y-Junction Splitters

Abstract

The design of passive photonic devices based on geometry optimization can lead to energy-efficient, small-footprint, and fabrication-ready geometries. In this work, we propose an angle-based parametrization method to optimize Y-junction splitters based on multimode interferometers. The selected figure of merit was the transmittance in the SCL and OESCL optical fiber communication bands. The performances of three optimization methods were compared: (i) particle swarm optimization (PSO), (ii) genetic algorithm (GA), and (iii) the covariance matrix adaptation-evolution strategy (CMA-ES). The results show that CMA-ES parametrization produces similar transmittance results (≤1.5% of difference) to PSO in the first 40 generations. The CMA-ES results are identical in the SCL (1460–1625 nm) and OESCL (1260–1625 nm) bands, whereas the GA and PSO executions are slightly different in terms of the rate and similarity of the figure of merit.

1. Introduction

Photonic integrated circuits (PICs) enable the interconnection of different types of photonic devices into a single platform. Designing energy-efficient on-chip optical interconnects remains a challenge when scaling the number of devices per chip [1]. Many optical passive devices work traditionally in the SCL optical communications band, while the OESCL band has gained attention due to ultra-wideband applications [2].

Among the passive components on a chip, the 3 dB power splitters are present in the order of dozens [3,4,5]. A 3 dB power splitter distributes the signal across multiple optical interconnects and is the essential building block for switches [6,7], modulators [8,9], and multiplexers [10,11]. Y-junctions based on multi-mode interferometers have become one of the essential building blocks in photonic integrated circuits [12]. Moreover, many inverse design methodologies for power splitters are based on multi-mode interferometer Y-junctions [10,12,13].

There are three main causes for excess loss in fabricated Y-junctions [12,14]: (i) a low tolerance to fabrication errors, (ii) minimum fabricated resolution available, and (iii) deviations in the sharp angles (see ${\alpha }_i$ in Section 2.1). Special mask techniques can to help overcome these problems, e.g., the electron beam technique [13,14]. In this regard, any device geometry, regardless of the type of optimization, has to be adjusted to the fabrication tolerances of the machinery available to the designer. This is especially true for non-intuitive geometries generated by optimization methods. Post-processes are usually added [13,15], in which the fabrication constraints are fed to a second algorithm to smooth and tweak the optimized geometry, making it fabrication-ready.

Figure 1 shows the timeline of some popular optimization algorithms used in the photonic design and relates them to their adoption (by the number of citations). We observe that the most popular algorithms are particle swarm optimization (PSO) [16] and genetic algorithms (GAs) [17]. These algorithms are usually used with default parameters due to the computational complexity of performing a hyperparameter optimization, which can limit their performance [18]. The covariance matrix adaptation evolution strategy (CMA-ES) [19] is a global optimization algorithm that has not been widely adopted. However, it does not need to optimize its parameters—only its population size. This characteristic is important because increasing the population size usually leads to a more exhaustive search. Furthermore, CMA-ES has been demonstrated to be efficient, robust, and easy for antenna optimization [20]. This motivates us to compare the performance and convergence of CMA-ES, GA, and PSO.

Figure 1 includes other optimization methods: (i) moving asymptotes (MMA) [21], (ii) limited memory Broyden–Fletcher–Goldfarb–Shanno bound constrained (L-BFGS-B) [22], and (iii) differential Evolution (DE) [23]. For a detailed review on optimization methods, please refer to [18,24,25].

Table 1. Comparison of photonic devices designed using iterative optimization algorithms.

Work	Device	Method	Bandwidth	Loss	Size
Mak 2016 [28]	2  × 2 splitter	PSO	1530–1570 nm	S: 4.11 dB E: 6 dB	4.8 × 4.8 $\mathsf{\mu }$m${}^2$
Lu 2018 [29]	Beam splitter	PSO	1500–1600 nm	S: no losses	2 × 2 $\mathsf{\mu }$m${}^2$
Chang 2018 [10]	3 dB Y-splitter	Subwavelength inverse design	1520–1580 nm	E: 1.5 dB	2.88 × 2.88 $\mathsf{\mu }$m${}^2$
Xu 2022 [30]	3 dB Y-splitter	ADTO	1588–2033 nm	S: 0.83 dB	5.4 × 2.88 $\mathsf{\mu }$m${}^2$
Xu 2017 [31]	3 dB Y-splitter	Nonlinear fast search	1530–1560 nm	E: 0.96 dB	3.6 × 3.6 $\mathsf{\mu }$m${}^2$
Piggott 2017 [13]	1 × 3 splitter	Inverse design	1400–1700 nm	E: 0.64 dB	3.8 × 2.5 $\mathsf{\mu }$m${}^2$
Wang 2016 [32]	3 dB Y-splitter	Taper design	1530–1600 nm	E: <0.19 dB	5 $\mathsf{\mu }$m length taper
Tahersima 2019 [33]	3 dB Y-splitter	DNN inverse design	1450–1650 nm	S: 0.45 dB	2.6 × 2.6 $\mathsf{\mu }$m${}^2$
Our work	3 dB Y-splitter	CMA-ES	1260–1625 nm	S: 0.15 dB	2 × 1.2 $\mathsf{\mu }$m${}^2$

Losses (S) means the simulation results whereas response (E) signifies the experimental results. PSO: particle swarm optimization. ADTO: analog and digital topology optimization. DNN: deep neural network.

shows examples of photonic splitters designed using different algorithms to those presented in this work. Our CMA-ES-based optimized splitter presents a wider bandwidth along with a shorter footprint, and we also expect lower losses.

Our work extends an angle-based parametrization [26]. It defines the device’s internal angles and requires fewer parameters than conventional topological optimization because its number of parameters depends on one dimension of the design space. In contrast, conventional topology optimization relies on two design space dimensions. Then, the optimization of the device is reduced to a medium-dimensional problem, allowing us to use global optimization [18]. The usage of global optimization is a key difference compared to other approaches, such as level-set and density topology, which usually rely only on local optimization [9,24,27].

To the best of our knowledge, we propose the first strict comparison of (i) the waveband analysis, (ii) convergence rate, and (iii) similarity of the figure of merit between three different executions of GA, PSO, and CMA-ES to provide PIC designers with adequate tools for fabrication-ready ultra-broadband 3 dB MMI (multi-mode interferometer) power splitters.

2. Angle-Based Parametrization

Our goal is to propose a parametrization that achieves functional devices with a high tolerance to fabrication errors (i) using a lower quantity (∼15) of parameters compared to the conventional methods, and (ii) by avoiding sharp angles.

2.1. Materials and Methods

We used a silicon-on-insulator (SOI) platform, which is extremely attractive because of its compatibility with microelectronics, introducing now-popular scalable photonic integrated circuits.

In the inverse design, we need to specify the permittivity distribution via parametrization. The most used methods are (i) the density topology and (ii) level-set. In both methods, the design region is divided into pixels, and a permittivity is associated to each pixel [24]. For more methods, please refer to [34].

A drawback of these methods is that they usually require a high number of parameters to have the flexibility of exploring desired fabricable devices. The higher the degrees of freedom, the better we can explore the designs. However, this process can require considerable computing resources (i.e., CPU- or memory) [18], leading to significant execution times.

We propose the angle-based parametrization depicted in Figure 2. The pink region depicts the silicon core and the white background color represents silicon oxide. We divided a rectangular MMI splitter design region into z equally spaced vertical-aligned segments. The segments have different widths. Therefore, the division gives us a set of angles between consecutive segments, where ${\alpha }_i{|}_1^p$ are the angles and optimization parameters, and $p=z-1$.

The waveguide width $W_g$, MMI body length, $L_{MMI}$, and gap $\delta$ can be adjusted to fit multiple outputs or operation frequencies. For simplicity, we demonstrate a 1 × 2 splitter working at 1550 nm.

A 1 × 2 MMI splitter parametrization can be summarized in three main steps:

(i)

First, we divide the initial topology into $W_i[i=1\cdots z]$ longitudinal equally spaced segments 1 based on [12,35,36], and we set symmetry on the transverse direction of the propagation. We also set the values of $\delta$ = 0.2 $\mathsf{\mu }$m and $L_{MMI}$ = 2 $\mathsf{\mu }$m. This allows us to reduce the number of angles ${\alpha }_i$ to optimize because the angles on the upper edge of the device match those at the bottom edge.

(ii)

Second, we modify the dimensions of the current width $W_i$ 2 based on the previous width and angle, as shown in Equation (1). In particular, we can set $W_1=W_g=0.5\mathsf{\mu }$m for the 1550 nm working window.

$$W_i=W_{i-1}+2d·tan\left({\alpha }_{i-1}\right)\mathrm{for}i>1$$

(1)

(iii)

Third, we place a power monitor 3 at one waveguide output to measure the transmittance of the fundamental TE (transverse electric) mode in order to establish the figure of merit.

2.2. Figure of Merit (FOM)

The choice of an FOM is the most critical parameter on the optimization setup. The FOM is a value used to characterize the performance of the design. We evaluated the Y-junctions efficiency by the transmittance over a set bandwidth. To improve transmittance without affecting the area on-chip, we limited the size of the MMI by restricting the parameter range of the angles and setting an FOM directly proportional to the transmittance of the fundamental TE mode over a frequency window. We ran the optimizations over SCL (conventional) and OESCL (wide) bands to obtain MMI geometries and evaluated the performance.

Our goal is to maximize the FOM. Let $T\left(f\right)$ be the transmittance evaluated at the frequency f. We propose the FOM as the transmittance accumulated from $f_1$ to $f_2$ shown in Equation (2).

$$FOM={\int }_{f_1}^{f_2}T\left(f\right)df$$

(2)

where $f_1=\frac{c}{{\lambda }_1}$, $f_2=\frac{c}{{\lambda }_2}$, the frequency is in Terahertz, and c is the speed of light in the vacuum. The values of ${\lambda }_1$ and ${\lambda }_2$ depend on the operating band. We define two FOMs:

(i)

For the $SCL$ band, ${\lambda }_1=1.625$ $\mathsf{\mu }$m and ${\lambda }_2=1.46$ $\mathsf{\mu }$m. The ideal value for the FOM is 10.425 assuming rectangular transmittance windows with a maximum amplitude of 50% on one output waveguide along the SCL band.

(ii)

For the $OESCL$ band, ${\lambda }_1=1.625$ $\mathsf{\mu }$m and ${\lambda }_2=1.26$ $\mathsf{\mu }$m. The ideal value for the FOM is 26.721 assuming a transmittance of 50% on one output waveguide along the OESCL band.

Although it is more geometry-restrictive than density topology and level-set parametrization, our approach shows promising results and potential for improvements (see Section 4). We produced fabrication-ready devices without any additional post-processing steps for smoothing or tweaking the optimized geometry. In the same fashion, other fabrication constraints can be morphed into optimization parameters, as we propose here, to produce faster and direct results.

3. Optimization Algorithms

This section details the three evolutionary algorithms used in this work: (i) GA, (ii) PSO, (iii) CMA-ES. Our implementations of these algorithms are used to independently optimize the two FOMs defined in Section 2.

3.1. Genetic Algorithm

GA is a metaheuristic inspired by the process of natural selection [17], and GA is used for the design of novel optic devices such as (i) polarization beam splitters [37,38] and (ii) Y-branch waveguides [27,33,39].

As depicted in Algorithm 1, GA has four main steps. First, it starts with n individuals represented by p parameters. Second, it will simulate the evolution of these individuals during k generations. In each generation, the simulation starts selecting a set of individuals to survive; those are called parents. Third, parents are combined following a crossover policy to produce n new individuals that share some parameters with the selected parents. Finally, the individuals randomly mutate some of their parameters.

Algorithm 1 depicts our implementation, with the following method specifications:

Algorithm 1: Framework of GA

• $generate\_population(n,p)$: returns n vectors defined by p angles (${\alpha }_i,i\in [1,p]$ defined in Section 2 with random values from a uniform distribution in the range of the constrains of each parameter.
• $select\left(population\right):$ chooses q individuals (i.e., we use roulette wheel selection) depending on the probability value $pro{b}_i$ given by Equation (3):

  $$pro{b}_i=\frac{FO{M}_i-min\left(FOM\right)}{\sum (FO{M}_j-min\left(FOM\right))}$$

  (3)

  where $FO{M}_i$ is associated with the i-th individual.
• $crossover(population,parents)$: returns n vectors of dimension p. The i-th vector is the combination of two random parents $p{a}_i$ and $p{b}_i$ selected from $parents$. Its d-th parameter is defined at random by $p{a}_{i,d}$ or $p{b}_{i,d}$ with equal probability (i.e., we are using uniform crossover). This generates new combinations of angles to be evaluated.
• $mutation\left(children\right)$: returns the new combinations of the previous step, but each individual may have some of its parameters added a random value $u\in U(-r,r)$. This updates the angles ${\alpha }_i$ of each individual.

3.2. Particle Swarm Optimization

In [16], the authors proposed PSO to simplify the social behavior of a flock of birds. It has shown success in discontinuous multidimensional problems in electromagnetism applications [18,20,27]. The PSO algorithm is similar to GA, but the i-th individual is related to a number ${\nu }^{\left(i\right)}$ and a p-dimensional vector $x_{best}^{\left(i\right)}$.

We can think of the i-th individual as a particle defined by: (i) its position $x^{\left(i\right)}$ (the p-dimensional vector associated with the individual), (ii) its velocity ${\nu }^{\left(i\right)}$, and (iii) the best position the particle has found so far, named $x_{best}^{\left(i\right)}$. Each particle accumulates speed in a favorable direction given by the best position found by the particle and the best position globally found. As a result, particles can move independently of local perturbations. By adding random walks, the particles can move in unpredictable ways that can potentially find better directions.

Using the framework of Algorithm 1, we can describe PSO methods as follows:

• $generate\_population(n,p)$: returns n particles defined by p angles (${\alpha }_i,i\in [1,p]$) defined in Section 2, where its attributes have random values from a uniform distribution in the range of the constraints of each parameter.
• $select\left(population\right)$: returns the particle with the best FOM.
• $crossover(population,parents)$: returns the population updated with the use of Equation (4) (this updates the angles ${\alpha }_i,i\in [1,p]$ of each of the n particles) and Equation (5):

  $$x^{\left(i\right)}\leftarrow x^{\left(i\right)}+{\nu }^{\left(i\right)}$$

  (4)

  $${\nu }^{\left(i\right)}\leftarrow \omega {\nu }^{\left(i\right)}+c_1{r}_1\left(x_b^{\left(i\right)}-x^{\left(i\right)}\right)+c_2{r}_2\left(x_b-x^{\left(i\right)}\right)$$

  (5)

  where $x_b$ is the best position globally found, $\omega$, known as inertia weight, represents the will of the particle to conserve its current velocity, $c_1$ and $c_2$ quantify the relative attraction of $x_b^{\left(i\right)}$ and $x_b$, respectively, and $r1,r2\in U(0,1)$ represent the unpredictable behavior.

3.3. Covariance Matrix Adaptation Evolution Strategy

CMA-ES [19] is an evolutionary strategy that has shown promising results for the design of antenna devices [20], splitting/combining [40], and grating filters [41]. The general idea of CMA-ES, shown in Algorithm 2, is to maintain: (i) an p-dimensional vector $\mu$ where p is the number of parameters (i.e., number of angles), (ii) a matrix C, and (iii) a number $\sigma$ to generate n sample points from a distribution $\mathcal{N}(\mu ,{\sigma }^{2}C)$.

Algorithm 2: CMA-ES algorithm

Taking points from this distribution limits the search space to a hyperellipse. Then, the algorithm evaluates points from the hyperellipse. Using the obtained values, it can decide to: (i) move the hyperellipse to another search space region and (ii) expand or reduce the covered region. CMA-ES works iteratively on this idea until the hyperellipse almost degenerates to a point, which is potentially the global optimum. For a detailed description of the algorithm, please refer to [42].

In Algorithm 2, we consider global variables for the sake of simplicity. We summarize CMA-ES in five steps. In line 1, we repeat the same procedures k times. In iteration t, we start with line 2 sampling n random $p-$dimensional (each dimension corresponds to an angle ${\alpha }_i$ defined in Section 2) points $x_i$, where x is sorted in descending order by evaluating the FOM. The points are obtained from the distribution $\mathcal{N}(\mu ,{\sigma }^{2}C)$. The next steps update the parameters $\mu ,\sigma$, and C such that we can evaluate different angle parametrizations (potentially with better FOM values) in the next iterations. In line 3, we update the mean $\mu$ using a weighted average of the sample points; see Equation (6).

$${\mu }^{(t+1)}\leftarrow \sum _{i=1}^n{w}_i{x}_i$$

(6)

The coefficients $w_i$ are fixed and chosen to proportionate larger contributions to the best-performing points; then, the mean $\mu$ is moved to a favorable region. Next, we need to properly change $\sigma$ to reduce or expand the hyperellipse of the next iteration. Therefore, in line 4, we control $\sigma$ using Equations (7) and (8).

$${\sigma }^{(t+1)}\leftarrow {\sigma }^{\left(t\right)}exp\left(\frac{c_{\sigma }}{d_{\sigma }}\underset{\mathrm{evolution}\mathrm{path}\mathrm{comparison}}{\underbrace{\left(\frac{\left|\right|p_{\sigma }\left|\right|}{\mathbb{E}\left|\right|\mathcal{N}(0,I)\left|\right|}-1\right)}}\right)$$

(7)

$$\mathbb{E}\left|\right|\mathcal{N}(0,I)\left|\right|=\sqrt{2}\left(\frac{\Gamma \left(\frac{p+1}{2}\right)}{\Gamma \left(\frac{p}{2}\right)}\right)$$

(8)

where $p_{\sigma }$ is a cumulative variable that tracks steps over time, $c_{\sigma }\in [0,1]$ is a variable that determines the cumulation time for $p_{\sigma }$, and $d_{\sigma }\approx 1$ is a damping parameter that determines the possible rate of change in ${\sigma }^{(t+1)}$. The key part of Equation (7) is the evolution path comparison; here, the length of $p_{\sigma }$ is compared to its expected length under random selection. From this comparison, we can control $\sigma$ as increasing, decreasing, or keeping its value. Finally, in line 5, we change $\Sigma$ in a favorable direction using Equation (9).

$$\begin{array}{c}{C}^{(t+1)}=\stackrel{\mathrm{cumulative}\mathrm{update}}{\overbrace{\left(1-c_1{c}_c(1-h_{\sigma })(2-c_c)-c_1-c_{\mu }\right)C^{\left(t\right)}}}\hfill \\ \hfill +\underset{\mathrm{rank}-\mathrm{one}\mathrm{update}}{\underbrace{c_1{p}_C{p}_C^T}}+\underset{\mathrm{rank}-\mu \mathrm{update}}{\underbrace{c_{\mu }\sum _{i=1}^n{w}_i^{\prime }{\delta }^{\left(i\right)}{\left({\delta }^{\left(i\right)}\right)}^T}}\end{array}$$

(9)

where $c_{\mu }\le 1$ is the learning rate for the rank-μ update, $c_1\le 1-c_{\mu }$ is the learning rate for the rank-one update, $c_c\in [0,1]$ is the learning rate for the cumulative update, $h_{\sigma }$ is the evaluation of a unit step function used to properly update the evolution path, $p_C$ is a cumulative vector used to update the covariance matrix, $w_i^{\prime }$ is a modification of the coefficients $w_i$, and ${\delta }^{\left(i\right)}$ is the sampled deviations.

In Equation (9), the first term (cumulative update) retains the information of the previous covariance matrix. The second term (rank-one update) allows for elongating the distribution along a favorable axis. The third term (rank-μ update) increases the search in areas where good solutions are likely to exist. The combination of the three terms changes C to move the hyperellipse in a favorable direction.

3.4. Algorithm’s Implementation

Our angle-based method (Section 2) was optimized with each of the algorithms using the parameters in

Table 2. Algorithm’s parameters.

Parameter	Value	Description
n	40	population size
k	100	number of iterations
p	11	number of characteristics
q	10	number of selected parents
r	3	range of random change in parameters
$\omega$	1	inertia weight
$c_i$	1	relative attraction coefficients

following the variables described in the three previous sections. We used $p=11$ angles to parametrize the splitter. For the optimization, we executed the algorithms $k=100$ times and each step works with $n=40$ exploratory designs. The parameter r is described in Section 3.1 and the parameters $\omega ,c_i$ in Section 3.2.

Let us impose the non-optimized geometry limits of our 3 dB Y-junction: a 2000 nm long and 1500 nm wide rectangular multi-mode interferometer (when ${\alpha }_i=0$ and all $W_i$ are the same). Thus, using $p=11$ the feature size of each subsegment is well above 120 nm, which is the limit of fabrication of many UV-lithography processes even by multi-project wafer standards [43].

4. Results

We executed each algorithm with a fixed number of iterations using the parameters in Table 2. The executions were run three times in the bands SCL for commonly optical interconnected applications and OESCL for ultra-wideband applications. We used this differentiated frequency window approach to test the robustness of the transmittance integral defined in Equation (2). Lumerical’s 2.5D variational FDTD was run on the solutions. All optimized designs depend on the initial population and thus on the initial guess. Before our work evaluation, we performed an exploratory analysis for validation, executing an optimization for all algorithms four times each. This was performed individually with a different initial random guess because we do not want to make additional assumptions about the final design and wanted to simply let the optimization algorithms find some patterns.

Our results are presented in three main terms:

(i)

Convergence rate is the algorithm’s speed in achieving FOM stability, with a variation <2%. To calculate it, we used the mean (Geomean) of the FOM trend.

(ii)

To evaluate the similarity of the FOM, we defined the deviation factor $\Delta FOM$ as the absolute difference between the geometric mean, Geomean, and the actual FOM trend value, $\Delta FO{M}_{i,j}=|{\mathrm{Geomean}}_i-FO{M}_{i,j}|$, where j is each algorithm execution, i.e., PSO${}_1$, PSO${}_2$, and PSO${}_3$, and i is the type of algorithm, i.e., Geomean${}_{PSO}$, Geomean${}_{GA}$, or Geomean${}_{CMA-ES}$. Additionally, the maximum FOM was also considered to characterize the performance of each method.

(iii)

The waveband analysis was performed over a defined band, either SCL or OESCL. We analyzed the maximum and minimum transmittance values to characterize a planar response of the device over the given frequency. We considered variations ≤1% of the transmittance to be under the fabrication error limit.

The average optimization total execution time was 12 h for a 100-iteration analysis for each algorithm. Our experiments were performed on an individual workstation with an Intel i7 processor 4th gen and 32 GB of RAM.

4.1. Convergence Rate

Figure 3 depicts the FOM trend for each optimization run over the SCL (top) and OESCL (bottom) bands. From Figure 3 (top), we make two key observations regarding the Geomean. First, ${\mathrm{Geomean}}_{CMA-ES}$ shows that CMA-ES is the fastest algorithm, converging at the 40th generation. Notice that PSO has a high initial growth rate; however, the <2% criteria is only met in the 71st iteration. Second, GA exhibits the slowest convergence rate; it converges at the ∼90th iteration, obtaining the lowest FOM.

From Figure 3 (bottom), we make three observations. First, in the OESCL band, the fastest method, regardless of its maximum FOM (25.6), is the GA, which converges at the 20th generation. However, as seen here and in the SCL analysis, the FOM values constantly get stuck. This is a drawback for the GA, even after 100 generations. Second, ${\mathrm{Geomean}}_{CMA-ES}$ achieves a maximum FOM at the 42nd generation; however, its maximum is just behind that of ${\mathrm{Geomean}}_{PSO}$. Lastly, the PSO algorithm achieved a maximum FOM over 26, which is greater than FOM${}_{\mathrm{GA}}$ at the ∼50th generation.

We conclude that the CMA-ES algorithm is the fastest and has a similar FOM to PSO in the SCL-band analysis. In the OESCL-band, its maximum FOM decays ∼2% in relative amplitude, and its convergence rate stays almost the same at the ∼40th generation. This consistency is a desirable feature for a global optimization algorithm.

4.2. Similarity of the FOM

We expect that $\Delta FOM$ decreases as the number of iteration increases; however, this is not the case for the GA in the SCL band depicted in Figure 4 (top). From this figure, we make two key observations. First, GA has the most dissimilar FOM trends among its three executions. Each GA execution finds a different best parametrization multiple times through its 100 generations. Second, CMA-ES and PSO have a clear stabilization of $\Delta FOM$. This is especially true for CMA-ES after the 40th generation. The similarity of the FOM is down to 0.05, and this result constitutes the best one for the SCL band analysis.

From Figure 4 (bottom), we make the following observations. First, the PSO has the best similarity of FOM in the OESCL analysis, with $\Delta FOM\le 0.1$ after convergence at the 18th generation. Second, GA and CMA-ES algorithms have around the same similarities as FOM in the convergence region. However, CMA-ES keeps the same pattern as the SCL analysis, while the GA cannot be directly compared. We conclude that PSO and CMA-ES results are more congruent in their pattern of convergence than the GA, even after evaluating them with different FOM integrals.

4.3. Waveband Analysis

Zhang et al. [12] proposed a longitudinal segmentation approach for the inverse design of Y-junction splitters using PSO over the 1500–1580 nm band. This segment-based approach resembles our proposal. However, it requires additional post-processing, such as geometry smoothing. Figure 5 depicts the waveband comparison between [12] and our CMA-ES best result. Our approach shows an improved transmittance of <1550 nm wavelengths. The footprint of our device is only 5% larger than Zhang’s device, and it does not need additional processing. These features motivate us to evaluate our method’s performance on a broader band.

Based on our definition of the FOM, which considers the integrals over the SCL and OESCL frequency domain, we expect planarity from the SCL band, given that it represents the narrower band between the two.

At first glance, the results shown in Figure 6 confirm our expectation. Note that the figure depicts the best and worst executions for each algorithm. Additionally, we calculated the difference between the minimum and maximum values of transmittance in the waveband for each of the three different runs. We make two key observations. First, in the SCL band, PSO is ≤1%, GA is ≤1.5%, and CMA-ES is also ≤1%. All executions show similar transmittances, and their differences are around the fabrication-error limit. Second, the OESCL results show that the GA is up to 4%, whereas PSO is ≤1.4% and CMA-ES is ≤2%, respectively.

Figure 7 depicts the electric field distribution at 1330 nm and 1550 nm. It shows that the gap between the arms ($\delta$ in Figure 2) is more constraining than the angles of the junction. This is because of our angle-based parametrization method. We observe a smooth transition from the MMI region through the two spline-shaped arms in both bands. Fabrication is feasible using photolithography.

We conclude that even when the FOM integral of the OESCL band defined in Equation (2) goes up to more than 2.5 times the value of the integral of the SCL band, only the transmittance of the GA executions is affected. The effect on the other two methods remains under the fabrication error limit.

5. Conclusions

We proposed the angle-based parametrization of a 3 dB Y-junction splitter that produces foundry fabrication-ready devices and optimized new geometries evaluating three evolutionary algorithms: PSO, GA, and CMA-ES. We evaluated the algorithms considering their convergence rate, the similarity of FOM, and the performance of their waveband analysis.

We found that although CMA-ES is a global optimization algorithm, it can be applied to photonic devices’ topological optimization just by defining the population size. Our analysis shows that CMA-ES parametrizations produce similar transmittance results to PSO but in fewer generations. This result persists in optimizations executed in the SCL and OESCL-bands, whereas GA and PSO executions fluctuate.