Hybrid Optimization of Photonics Crossing WaveguidesAbstract
This work investigates the optimization of silicon waveguide crossings using a hybrid GRASP-Simulated Annealing (GRASP-SA) algorithm coupled with two-dimensional finite element method (2D-FEM) analysis. An inverse design approach was adopted to systematically modify geometric parameters and enhance transmission efficiency at the design wavelength of 1.55 μm. Different values of the control parameter μ were evaluated, and the best-performing configuration was selected for further spectral analysis. The optimized structure achieved transmission efficiencies above 97%, corresponding to approximately 0.13 dB insertion loss. A wavelength sweep from 1.50 to 1.60 μm confirmed the stability of the optimized geometry. The results demonstrate that the proposed hybrid metaheuristic framework is an effective alternative for the design of high-performance photonic waveguide crossings.
Index Terms
Bio-Inspired Algorithms; Integrated optics devices; Hybrid GRASP-simulated annealing algorithm; Numerical approximation and analysis.
I. INTRODUCTION
Inverse design has become a powerful paradigm for the development of advanced photonic structures, particularly in scenarios where conventional intuition-based approaches are insufficient due to the complexity of electromagnetic interactions and multi-parameter dependencies. Early applications of metaheuristic optimization in photonics demonstrated the effectiveness of genetic algorithms in the design of two-dimensional photonic crystals and dispersion-engineered devices [1]-[3], [5]-[6], [12]. In parallel, alternative population-based and combinatorial strategies, such as scatter search [4], swarm intelligence methods including artificial bee colony algorithms [7]-[8] and particle swarm optimization [18], and artificial immune system approaches [9]-[10], further expanded the range of bio-inspired techniques available for inverse photonic design. These methodologies have also been successfully applied to plasmonic absorbers [11] and other complex electromagnetic structures.
More recently, inverse design has been recognized as a central tool in nanophotonic and intelligent material engineering [13]-[14], enabling the systematic discovery of innovative device geometries. In the specific context of waveguide crossings, comprehensive reviews and high-performance implementations have been reported in [15]-[17], highlighting the importance of minimizing insertion loss and crosstalk in optical integrated circuits. An inverse-designed binary waveguide crossing using particle swarm optimization combined with three-dimensional FDTD simulations was presented in [18], achieving compactness but requiring computationally intensive full-wave modeling.
From a methodological perspective, the Greedy Randomized Adaptive Search Procedure (GRASP) was originally introduced as a constructive metaheuristic incorporating adaptive randomized mechanisms [22]-[25]. Simulated Annealing (SA), on the other hand, is grounded in statistical mechanics principles and employs probabilistic acceptance criteria based on the Boltzmann distribution to escape local minima [26]-[29]. Hybrid metaheuristic frameworks and multi-level optimization strategies have been investigated in other engineering domains [19], and general theoretical foundations of metaheuristics are consolidated in [20]-[21]. Recent works continue to refine and extend GRASP- and SA-based approaches [34]-[37], demonstrating their adaptability to complex, high-dimensional optimization problems.
In parallel, state-of-the-art waveguide crossing designs have recently advanced through inverse and topology-based strategies. Zhou et al. [30] demonstrated ultra-broadband and low-loss waveguide crossings obtained through inverse design, with improved tolerance to fabrication imperfections. He et al. [31] introduced valley-Hall topological concepts to enhance robustness against structural variations. Zhong et al. [32] combined the adjoint method with direct binary search to balance footprint and optical performance, while Nan et al. [33] employed a three-dimensional Ta₂O₅-on-LNOI platform to achieve near-zero crosstalk and extremely low insertion losses. These works illustrate the diversity and sophistication of current optimization approaches in photonic integration.
Despite these significant advances, the direct integration of GRASP and Simulated Annealing within a unified framework specifically tailored for the inverse design of planar silicon waveguide crossings has not been systematically explored. Moreover, many high-performance approaches rely on gradient-based optimization or computationally demanding three-dimensional simulations. There remains a need for alternative hybrid strategies capable of efficiently navigating the design space while maintaining computational feasibility in planar geometries.
In this work, we propose a hybrid GRASP-SA algorithm combined with two-dimensional finite element method (2D-FEM) analysis in the frequency domain for the optimization of silicon-on-insulator waveguide crossings. The constructive randomized exploration characteristic of GRASP is integrated with the probabilistic refinement mechanism of SA, enabling a balance between diversification and intensification during the search process. The electromagnetic performance of candidate geometries is evaluated using 2D-FEM, providing a computationally efficient yet physically consistent framework for analyzing the targeted wavelength range. By optimizing both transmission efficiency and crosstalk at the design wavelength of 1.55 μm and subsequently performing spectral analysis in the 1.5-1.6 μm range, the proposed approach achieves transmission efficiencies above 97% (0.13 dB insertion loss) and crosstalk levels below -50 dB within a compact 4.9 × 4.9 μm2 footprint.
The proposed methodology establishes a coherent link between classical metaheuristic theory [22]-[29], recent hybrid developments [34]-[37], and contemporary photonic inverse design strategies [15]-[18], [30]-[33]. By providing a computationally efficient hybrid framework tailored to planar waveguide crossings, this work contributes an alternative optimization strategy for high-performance silicon photonic integration.
II. OPTIMIZATION STRATEGY
To solve combinatorial problems, in which there is a large set of variables, metaheuristics can be used. By an ideal adjustment of parameters, the search method can have its performance improved. Therefore, it is of fundamental importance to perform tests to validate parameters that can meet the requirements of the application of interest [19]-[20]. Topology optimization methods have proven to be a promising answer for solving problems in photonic system designs. The optimization procedure is a tool that can facilitate the design of optical components, contributing to the manufacturing process and refining the final performance, as parameter adjustments are made in the device structure [21].
To initiate the optimization process, Fig. 1 was used as the base structure (starting point of the hybrid algorithm), featuring variations in the geometric parameters of the waveguide's cross section. In addition, a guide width of w = 0.4 μm and the refractive indices of n1 = 3.476 and n2 = 1.444 were considered, for silicon and silica materials, respectively.
In Fig. 1, the input and output ports are identified by the optical powers P1 and P2. The excitation is applied at the left horizontal waveguide (P1), and the transmitted power is measured at the right horizontal waveguide (P2, corresponding to j = 2). The vertical arms are used for crosstalk evaluation.
Changes were implemented in the guide cross-section, in order to maximize the efficiency in the output power, thirty-two points of the structure highlighted in Fig. 1 were subjected to variations. The waveguide width was preserved without modifications at w = 0.4 μm.
The algorithm is adjusted to the device's geometry to achieve the desired performance, and transmission behavior is analyzed. Various photonic structure designs are executed using software, enabling the optimization of structures with millions of geometric parameters through small structural variations that can influence the objective function value and consequently update the design. Multiple geometries are tested in parallel, and the best ones are selected, speeding up the process. Among the advantages is the generation of non-intuitive designs that would be very difficult to create manually.
The purpose of the optimization strategy is to identify a set of points that maximize transmission efficiency. Due to the structure's 90° rotational symmetry, it is only necessary to find 4 points; the remaining 28 points are determined by 4-fold symmetry around the origin. To execute the iterative optimization process, the following coordinates associated with the points were considered and are represented in Table I.
The proposed structure exhibits fourth-order rotational symmetry (90°), a characteristic sufficient to guarantee the invariance of the device under rotations. Furthermore, the geometry has reflection symmetries with respect to the horizontal and vertical axes, such that the structure is mirrored: what occurs in the lower region is repeated in the upper region, as are the behaviors observed on the left and right. As an example, the regions identified by points 1 to 4 in Fig. 1 correspond to mirror images of regions 29 to 32. These symmetries arise naturally from the constraints imposed on the design and contribute to a balanced distribution of fields and equivalent optical responses along the intersecting waveguides. For reasons of clarity and conciseness, only rotational symmetry was explicitly mentioned in the original description.
The adapted GRASP-simulated annealing (SA) hybrid algorithm was used to optimize the power transfer given by Eff = Pj/P1 with j = 2 representing the output port. Only port 2 performs an efficiency assessment, as it measures the amount of energy reaching the destination. Therefore, at the other ports, it is called crosstalk. The higher the efficiency and the lower the crosstalk, the better the crossover performance. The selected coordinate interval is essential to guarantee high efficiency values.
The efficiency is defined as Eff = Pj/P1, with j = 2, where P1and P2denote the optical powers measured at ports 1 and 2, respectively. Throughout this work, the subscripts associated with Pj refer exclusively to port numbering and should not be confused with the indices used to label geometrical points in Fig. 1 and Table I.
III. GRASP ALGORITHM
GRASP is especially attractive when the project can be written as a combinatorial/discrete problem: it builds solutions by “greedy” steps with controlled randomness (avoiding getting stuck on the same choice) and then applies local search to refine. This logic combines very well with pixelated/binary inverse design (e.g., regions with/without material, holes/without holes, or discrete choices of geometries) typical of many integrated components, where each iteration tries to improve loss and crosstalk without exploding the computational cost. A current and well-established basis for the use of GRASP (including path relinking, which explores trajectories between “elite” solutions) is discussed in the review by Laguna et al. [36], which systematizes the variants and good implementation practices, useful when you want to adapt GRASP for multiple objectives and rigid constraints, exactly the type of scenario found in compact crossings. Additionally, Resende's entry [37] brings together the notedness and structure of GRASP, serving as a methodological reference when justifying the algorithm in the work's "method", and helping to explain why it is suitable for PIC architectures where the design of a crossing can be viewed as a discrete selection/arrangement of geometric elements under optical metrics and manufacturing limitations.
Grasp is a search algorithm designed to solve combinatorial optimization problems [22]-[23]. The algorithm helps in the more efficient construction of solutions, by integrating greedy, random and adaptive search strategies, with the aim of constructing a trajectory in search of optimal solutions. In contrast to evolutionary algorithms, in GRASP there is no parallel with nature regarding the adaptability of individuals, to enable survival in the environment; instead, at each iteration, viable solutions are constructed in order to obtain a greedy, random and adaptive process. On the other hand, at each iteration, feasible solutions are constructed so as to obtain a greedy, random and adaptive process [24].
Through a greedy function, an evaluation of the available element options is performed and the specificity of being "greedy" is associated with the search for the most profitable options, which is measured through the function [23].
The randomness process is employed by a probabilistic method associated to a parameter μ, which has the function of managing the Restricted Candidate List (RCL), which stores the best evaluated solutions. Candidates with higher quality have an increased probability of selection, similar to the roulette process used in Genetic Algorithm. Population heterogeneity is the result of the selection of the list of candidates, which may be linked to different types of solution quality, which may deviate from local optima [23].
The heuristic procedure has this peculiarity, in which each component is updated at each iteration, that is, the greedy function associated with the element is updated [23]-[24]. While the process of executing the steps of the algorithm is taking place, elements are selected to be inserted in the group of solutions and with that there is an adaptation for the emergence of new solutions, which will be able to maximize or minimize the greedy function, becoming intelligent and thus if an improved solution is found, there is an update of the solution [25].
The method follows a repetitive process, divided into two phases at each iteration. In the first phase, an initial viable solution is constructed progressively, where, at each iteration, a new element is inserted and the next ones are defined. This occurs through an ordering in a restricted list of candidates, according to an adaptive greedy function, where it is possible to evaluate the elements and select them. In a next phase, the local search is performed in relation to the neighborhood of the initially generated solution. A cost function is used to evaluate neighboring candidates and identify a local optimum. If there was an improvement in the solution found, it replaces the current solution and thus the best result is saved [23]-[25]. In each iteration of GRASP, the technique disregards the results of previous iterations. In other words, it starts from scratch each time, characterizing it as a repetitive sampling strategy [25].
The construction of the RCL employs minimum and maximum thresholds to establish the elements that will be selected for inclusion in the list [24]. Considering that the variables Emax and Emin, indicating the maximum and minimum efficiency, respectively, as defined in (1) and (2), identified by the set of solutions, where i symbolizes each element of the solution set of group G.
The RCL must be formed by the i elements associated with the group G and to manage the number of elements that will integrate the RCL, one of the mechanisms is based on the value [24]. To employ this procedure (3) was used to determine the variable k as the result of the operations with the minimum and maximum efficiencies achieved by the group of elements belonging to the G group.
The cardinality control of the elements of the RCL is done by the parameter µ, where µ ∈ [0, 1], the choice is made through a uniform discrete probability distribution [24], consequently the probability of the selection one value for µ belonging to [0,1] is the same. The optimal µ values obtained in the process are highlighted in the final result.
In order to quantify the elements that will be present in the restricted list, values of k will be obtained, and for a solution to belong to the list, it must have at least the value of k, that is, for the solution to meet this criterion, it must be present in the set E(i) ∈ [k, Emax].
The choice of GRASP in the proposed hybrid framework is justified not by the dimensionality of the search space alone, but by the intrinsic nonlinearity and multimodal nature of the optimization problem. Although the number of geometric parameters may appear moderate, the objective function depends on full-wave electromagnetic simulations performed through 2D-FEM analysis. In such problems, small variations in geometric dimensions can lead to highly nonlinear variations in field distribution, modal coupling, and interference patterns. Consequently, the resulting design landscape is typically non-convex and characterized by multiple local optima.
GRASP was selected because of its constructive randomized mechanism, which promotes diversified exploration of the solution space without requiring gradient information. This characteristic is particularly advantageous when the objective function is implicitly defined by numerical electromagnetic solvers, where analytical gradients are not readily available. Moreover, the adaptive nature of GRASP allows the progressive refinement of promising regions of the search space, making it suitable even when the parameter space is not extremely large but exhibits strong nonlinear behavior.
When combined with Simulated Annealing, which introduces probabilistic acceptance of uphill moves based on statistical mechanics principles, the hybrid GRASP-SA strategy enhances the ability to escape local minima and improve solution robustness. Therefore, the selection of GRASP is motivated by the nonlinear electromagnetic response of the structure and the need for a gradient-free, exploration-oriented optimization strategy, rather than solely by the apparent size of the search space.
IV. SIMULATED ANNEALING ALGORITHM
In photonic design problems, such as waveguide crossings, the search space is often highly nonlinear and full of local minima (discretized geometries, fabrication constraints, and multiple objectives such as insertion loss, crosstalk, and banding). In this scenario, Simulated Annealing (SA) is useful because it allows for the controlled acceptance of temporarily worse solutions to escape local optima and explore alternative configurations before "freezing" on a good final design. A recent example of validated SA use in optics appears in Chen et al. [34], who formulate a cost function and use SA to optimize parameters (phase in SLM) until they obtain a better match with the target pattern, showing in practice how well the method adapts to discrete variables and simulation-defined objective functions. Going more directly to the context of dense photonic circuits, design automation work for optical-on-a-chip networks also incorporate SA as a physical layout optimization step, precisely because of the impact of guide crossings on overall performance (insertion, accumulated loss and routing constraints), using SA to organize topology/routes with crossing awareness and reduce associated penalties [35].
The Simulated Annealing (SA) algorithm is inspired by the simulation of the behavior of physical systems due to its temperature [26]-[27]. If the system is at a high temperature, there is a disorder of a large number of atoms and consequently the energy of the system is also high [28] and vice versa, as the system has its temperature reduced, the movement of atoms tends to settle into lower energy states.
When there is no longer a probability of cost improvement, the system is said to be in a frozen state [27]. A system is said to be in thermal equilibrium at a temperature T, with probability of being with energy E, in state i is driven by a Boltzman distribution [28].
To simulate the behavior of atoms under thermal equilibrium, an algorithm was implemented assuming a fixed [29]. At the beginning of the process, we have an energy E0. When subjected to a random disturbance in the current configuration of the system, a change in the layout of the system is verified and with this a new resulting energy is calculated [28]. As a result of this change in the system, an energy variation (∆E) is generated which is given by (4), where Er is the energy resulting from the disturbance of the system and Ec is the current energy prior to the disturbance.
The symbol E0 denotes the initial energy of the system before the optimization process starts. During the iterative procedure, however, the energy prior to a given perturbation is represented by Ec, which corresponds to the current energy at that specific iteration. After applying a disturbance, the resulting energy is denoted by Er. Consequently, the energy variation is defined as ΔE = Er - Ec. This distinction allows a clear separation between the initial reference energy E0and the iteration-dependent energy Ec, avoiding possible ambiguity in the interpretation of the optimization steps.
If ∆E ≤ 0, it indicates that the current energy of the system is greater than or equal to the energy resulting from the perturbation. As a consequence, the new arrangement of the perturbed system is adopted as an entry point for the next step of the algorithm.
Otherwise ∆E > 0, the energy of the current arrangement of the system is greater than the energy of the new arrangement of the system, the variation must consider (1.5) as probability distribution to reach the Boltzman distribution [28] where kB is the Boltzman constant, given by 1.38064852 × 10-23 m2 kg s-2 K -1.
In the simulated annealing procedure, when a perturbation leads to an energy increase (ΔE > 0), the new configuration may still be accepted according to the acceptance probability given in Eq. (1.5), which is based on the Boltzmann factor exp(-ΔE/kBT). This probabilistic rule allows occasional uphill moves, preventing premature trapping in local minima and driving the stochastic evolution of the system toward a Boltzmann distribution at the prescribed temperature T.
In the simulated annealing procedure, when a perturbation leads to an energy increase (ΔE > 0), the new configuration may still be accepted according to the acceptance probability given in Eq. (1.5), which is based on the Boltzmann factor exp(-ΔE/kBT). This probabilistic rule allows occasional uphill moves, preventing premature trapping in local minima and driving the stochastic evolution of the system toward a Boltzmann distribution at the prescribed temperature T.
In the simulated annealing framework adopted in this work, the acceptance probability for energetically unfavorable moves is given by the Boltzmann factor, exp(-ΔE/kBT), where kB is the Boltzmann constant and T represents the control temperature parameter. This ensures that the introduction of kB is mathematically justified and consistent with the formulation of the simulated annealing algorithm.
In order to represent the cooling rate, a variable R is introduced to the system and is randomly generated within the range [0,1]. There is a comparison of the variable R with the value of P(∆E). If R ≤ P(∆E), the new configuration of the disturbed system is accepted and will be the input for the next step of the algorithm. Otherwise, the perturbed configuration is rejected and the current disposition is used as the entry point for the next step of the algorithm [28]. This process can be implemented in combinatorial optimization problems, where the objective is to define a set of variables (x1,x2,...,xN) that minimize or maximize an objective function. It is possible to create an analogy between this function and the energy for physical systems, since both can be decreased or increased [28].
In the simulated annealing implementation adopted in this work, the variable R is generated from a uniform probability distribution in the interval [0,1]. It is not drawn from a Gaussian or Boltzmann distribution. The Boltzmann distribution appears only in the acceptance criterion, through the Boltzmann factor exp(-ΔE/kBT), which defines the probability of accepting energetically unfavorable transitions when ΔE > 0. The random variable R is then compared with this acceptance probability: if R < exp(- ΔE/kBT, the new configuration is accepted; otherwise, it is rejected.
Initially, viable solutions are generated with the aim of optimizing the objective function. Through an iterative process, there is a movement towards disturbing the existing solutions so that they can generate new solutions, considering a previously established stopping criterion. In order to control the acceptance of solutions, there is a parameter called temperature that will determine which solutions should continue for the next iterations. If there is a reduction in the temperature value, there is a tendency to reduce the number of accepted solutions; otherwise, if there is an increase in the temperature value, the number of accepted solutions tends to increase. Therefore, it is very important that if there is a significant reduction in the temperature value, there is a proximity to the global optimal solution. In addition, there must be adequate management of the cooling so that it does not generate inappropriate computational effort time. At a given temperature, there is an approximation of the accepted solutions to the Boltzman distribution and it is possible to say that the system is in "thermal equilibrium" [28].
The purpose of the metaheuristic is to conduct a search within the neighbourhood of a single solution and movements in space will occur, analogous to concepts of thermodynamic systems and within this structure there are both components of diversification and intensification [28].
There is an improvement in local optimization when applying the SA algorithm, where there is a repetition of the initial solution and a mechanism of local changes until a better solution can no longer be generated. There is a disturbance on the generated solution, in the sense of making changes to it, in an attempt to avoid stagnation in a poor solution [27]. The idea is that it is possible to reach a better solution through a worse solution. As the iterative process progresses, there is a lower tendency for worse solutions to be accepted with the help of the probability conditioned on the Boltzman distribution, which will assist in the solution selection process and, at the end of the execution process, there is a tendency for the best solutions to be accepted.
By using this approach, it is expected to obtain significantly better results, due to it ability to escape from local optima. The advantage of the SA algorithm is not abandoning not so attractive solutions and returning to old movements and thus in the future advancing to a better cost of the objective function, this brings a differential in relation to local optimization [27].
V. CUSTOM GRASP AND SA HYBRID ALGORITHM
For the design of waveguide with curvature and a variation of 32 points in the crossing area, the GRASP and SA algorithms were integrated, in order to identify steps of both methods that were more aligned with the objective of the project, which is to find the structures that provide the best transmission efficiencies. The steps applied in the combination of the Grasp and SA algorithms are highlighted through the flowchart of Fig. 2.
1) Initialization of the temperature variable: A temperature variable was defined to control the acceptance or rejection of the available solutions. In the iteration process, a reduction in temperature occurs, it was established that the temperature starts with the value of 200 and ends with the value of 0, which in this case is the end of the program execution, coinciding with the execution of the total number of iterations and the number of generations. It is an analogy to the control of the thermodynamic system, which uses the variable to manage the optimization system, where it is possible to accept or reject objective function values, through probabilistic calculation.
2) Generation of solutions of the thermodynamic system: Creation of viable solutions, which represent the data of the waveguide crossing in analogy to the viable states of the thermodynamic system. The set of coordinates associated with the set of points of the crossing structure associated with Fig. 1 are created randomly, according to a previously defined interval.
3) Evaluation of solutions using the greedy function: Each solution generated in the previous step is evaluated using the objective function (Pout/Pin) - with the calculation performed using the 2D-FEM [4], [9], [10].
4) Construction of the candidate list: The parameters k and µ are used to control the presence of solutions in the candidate list, where the value of k is calculated according to (3) and the values of k are generated in the range from 0 to 1. For the solutions to be present in the list, it is necessary that the efficiencies associated with the solutions are in the range [k,Emax], where the value of k determines the minimum efficiency value to be maintained in the candidate list.
5) Ordering of the list of candidates: According to the values of the transmission efficiencies associated with the solutions analyzed previously, a classification of the solutions occurs and a queue is constructed, where the highest objective function values are inserted at the beginning of the queue and the lowest are inserted at the end of the queue.
6) Generation of neighborhoods for the candidate list: To expand the set of solutions, candidates neighboring the candidate list was generated. The neighboring solutions were generated through small movements in the coordinates of the candidate shortlist, which is an analogy to changes in the state of the thermodynamic system.
7) Solution acceptance control: There is a comparison between the existing solutions and the solutions generated in the neighborhood. If there is an improvement in the value of the greedy function associated with the new solutions generated, the neighboring solution is saved in the list. Otherwise, a probabilistic process is implemented so that there is an acceptance of the worst solutions, through an exponential function P(∆E) = exp(-[∆E /kBT]) associated with the Boltzman distribution, which manages the acceptance of the worst solutions, as there is a reduction in temperature, that is, the probability of accepting a bad movement decreases and when zero temperature is reached, there is no more acceptance of movements. In addition, there is a cooling rate R, generated, randomly, in the range [0,1]. This rate is compared with the value of P(∆E), if R ≤ P(∆E), the new configuration of the perturbed system is accepted and this new arrangement will be the input for the next step of the algorithm. ∆E represents the energy variation of the thermodynamic system, which analogously represents the variation of the transmission efficiency, resulting from the efficiency of the neighboring solution subtracted from the efficiency of the current solution.
8) Storage of the selected population: Upon reaching a temperature equal to zero, the state of the system is said to be frozen and there will be no further changes in temperature, thus the last set of solutions selected will be stored.
The hybrid GRASP-SA strategy does not decompose the electromagnetic problem itself into smaller physical subproblems; rather, it structures the optimization process into sequential constructive and refinement stages. The increase in the number of algorithmic stages per iteration does not imply a proportional increase in full-wave simulations. In the proposed implementation, each candidate geometry is evaluated using 2D-FEM, and the hybridization mainly affects how new candidate solutions are generated and refined. Since 2D-FEM is computationally less demanding than three-dimensional full-wave solvers, the overall computational cost remains manageable even with the additional refinement mechanism introduced by the SA stage.
Although the hybrid approach introduces extra decision steps within each iteration, these steps improve search efficiency by accelerating convergence toward high-quality regions of the design space. In practice, this reduces the number of ineffective explorations and mitigates premature convergence, which may occur in purely greedy or purely stochastic strategies. Therefore, the computational overhead associated with hybridization is compensated by a more structured and directed exploration process.
When compared to other optimization strategies reported in the literature, such as GA-based approaches [1]-[6], PSO combined with 3D-FDTD simulations [18], or adjoint-based and topology-driven methods [30]-[32], the proposed methodology offers a balance between computational efficiency and design flexibility. Gradient-based adjoint methods typically require sensitivity analysis and may be strongly dependent on initial conditions, while 3D full-wave approaches can be computationally intensive. In contrast, the present framework relies on a gradient-free metaheuristic combined with 2D-FEM analysis, which significantly reduces simulation time while maintaining physical consistency for planar SOI structures.
Therefore, the main advantages of the proposed technique are:
(i) a gradient-free hybrid optimization framework capable of handling nonlinear electromagnetic responses;
(ii) computational efficiency due to the adoption of 2D-FEM instead of full 3D simulations; and
(iii) the ability to simultaneously optimize transmission efficiency and crosstalk within a unified procedure. These aspects position the proposed approach as a competitive and computationally feasible alternative among existing inverse design strategies for waveguide crossings.
VI. NUMERICAL RESULTS
To optimize the waveguide crossing structure, the hybrid algorithm GRASP was used in conjunction with SA. In Fig. 3, it is possible to observe the transmission efficiencies for (a) µ = 0.34, (b) µ = 0.51 and (c) µ = 0.82, with values of 97.58% (0.11dB), 97.55%(0.11dB) and 97.56% (0.11dB), respectively. Table II presents the parameters associated with µ=0.34, responsible for generating the highest transmission efficiency value. To execute the algorithm, each iteration took an average of 5 minutes on a machine with an Intel (R) Core i7-7500 processor at 2.70 GHz, 16 GB of RAM, and a 2 TB HD.
Transmission efficiency as a function of the number of generations for (a) µ =0.34 (b) µ = 0.51 and (c) µ = 0.82.
Since the combination of algorithms was adopted to optimize the structure, it was necessary to execute more steps and therefore it is possible to verify in Fig. 3(a)-(c) that very high transmission efficiency values are obtained since the first generation, this is due to the number of filters implemented before finishing the execution of the generations, which made it possible to avoid very low efficiency values. Fig. 3 (a)-(c) show the evolutionary behavior of transmission efficiencies as a function of the number of generations, for each value of µ.
In Fig. 3, each color represents the set of individuals selected in a specific generation during the optimization process. A total of 200 generations were considered, and in each generation the 10 best individuals were selected according to the objective function. The color variation was adopted to facilitate visualization of the evolutionary behavior of the population throughout the optimization procedure, allowing the reader to observe the convergence trend of the algorithm over successive generations.
The parameter μ shown in the figure corresponds to a tested value of the algorithm control parameter and does not represent a physical quantity of the device. Therefore, the colored points are related to the evolutionary optimization process and not to spectral characteristics.
The spatial distributions, according to the value of µ obtained through the optimization process, are presented in Fig. 4 (a)-(c). It can be verified that 97.58% (0.11dB) of the power is transmitted. In Fig. 4 (a)-(c), the transmission efficiency in the wavelength range [1.50 µm,1.60 µm] was evaluated and it can be verified that a higher value power transfer occurs at 1.55 µm, as expected.
To execute the algorithm, 200 iterations were considered, where at each iteration a set of 10 viable solutions are generated and used for the next generation. In addition, at the end of the 200 iterations, the best 10 solutions of the entire process performed are generated. Each of these solutions underwent an evaluation process at wavelengths from [1.5 µm, 1.6 µm]. Each transmission efficiency value is associated with a set of solutions, since the solutions are associated with different geometries, this generates different efficiencies. In Fig. 5, the 10 best solutions obtained after the optimization process, the relationship between wavelengths in the range [1.5 µm, 1.6 µm] and the transmission efficiencies associated with each of these solutions were evaluated.
Transmission efficiency in the wavelength range from 1.50 μm to 1.60 μm for the 10 best solutions obtained for the value of μ that provided the highest transmission efficiency at 1.55 μm.
In the wavelength range from 1.5 to 1.6 µm, the dispersion of silicon and silica is weak, and the variation of their refractive indices is small. For this reason, constant refractive indices were used in the numerical simulations in order to simplify the modeling and focus on the influence of the device geometry. Moreover, both materials exhibit extremely low absorption in the near-infrared telecom band, allowing the imaginary parts of the refractive indices to be neglected. Under these conditions, material losses are negligible when compared to scattering effects introduced by the waveguide crossing itself.
Fig. 5 presents the spectral response of the optimized structure corresponding to the value of μ that yielded the highest transmission efficiency at the design wavelength of λ = 1.55 μm. Initially, different values of μ(0.34, 0.51, and 0.82) were tested during the optimization procedure, and the transmission efficiencies reported at the beginning of Section IV refer to this central wavelength. Among these cases, the value of μ that provided the best performance was selected, and the corresponding optimized geometry was then evaluated over the wavelength range from 1.50 to 1.60 μm. Therefore, Fig. 5 does not represent separate optimizations for each μ, but rather the broadband analysis of the single best-performing configuration. This clarification establishes a direct connection between the optimization stage and the subsequent spectral evaluation, ensuring coherence among the simulation steps described throughout the manuscript.
VI. CONCLUSIONS
This work presented the optimization of silicon waveguide crossings using a GRASP-SA hybrid algorithm combined with the two-dimensional finite element method (2D-FEM). The proposed strategy was employed to explore alternative geometrical configurations capable of improving transmission efficiency in integrated photonic devices.
Initially, a base geometry was defined and manually adjusted to ensure an adequate starting point for the optimization process. This preliminary stage contributed to accelerating convergence and allowed high transmission efficiencies to be observed from the early iterations. Subsequently, the inverse design procedure was carried out by introducing controlled modifications in the geometric parameters through the hybrid GRASP-SA framework.
The combination of GRASP and simulated annealing enabled an effective balance between diversification and intensification mechanisms during the search process. This hybridization allowed the algorithm to explore different regions of the solution space while avoiding premature convergence to local optima. When coupled with 2D-FEM analysis, the methodology efficiently evaluated candidate structures and guided the optimization toward improved geometrical configurations.
The optimized waveguide crossing achieved transmission efficiencies exceeding 97% (corresponding to approximately 0.13 dB insertion loss) at the design wavelength. Furthermore, the spectral analysis confirmed the robustness of the optimized structure within the investigated wavelength range. These results demonstrate that the proposed hybrid strategy is a viable and competitive alternative for the inverse design of high-performance photonic components.
Overall, the integration of metaheuristic optimization with numerical electromagnetic simulation provides a systematic and flexible framework for the development of advanced integrated optics devices, contributing to the diversification of design methodologies in silicon photonics.
ACKNOWLEDGMENT
The authors would like to thank CAPES Coordination for the Improvement of Higher Education Personnel (Grant 001), CNPq National Counsel of Technological and Scientific Development (Grants 303795/2022-0 and 407612/2025-4), FAPESB (Grant PIE0003/2022), UFBA and FUNDECT (Process: 83/049.252/2024).
DATA AVAILABILITY
The data supporting the findings of this study are available from the corresponding author upon reasonable request.
REFERENCES
- [1] G. N. Malheiros-Silveira, V. F. Rodriguez-Esquerre, and H. E. Hernandez-Figueroa, “Strategy of search and refinement by GA in 2-D photonic crystals with absolute PBG,” IEEE Journal of Quantum Electronics, vol. 47, no. 4, pp. 431-438, 2011.
- [2] G. N. Malheiros-Silveira and V. F. Rodriguez-Esquerre, “Photonic crystal band gap optimization by generic algorithms,” IEEE MTT-S International Microwave Symposium Digest, pp. 734-737, 2007.
- [3] S. Preble, M. Lipson, and H. Lipson, “Two-dimensional photonic crystals designed by evolutionary algorithms,” Applied Physics Letters, vol. 86, no. 6, pp. 061111-1-061111-3, 2005.
- [4] L. F. Vieira, V. F. Rodríguez-Esquerre, A. Dourado-Sisnando, and C. E. Rubio-Mercedes, “Inverse design of a taper by scatter search metaheuristic,” IEEE Photonics Journal, vol. 12, no. 4, pp. 1-9, 2020.
- [5] D. F. Rêgo, I. L. Gomes de Souza, and V. F. Rodríguez-Esquerre, “Ultra-broadband plasmonic groove absorbers for visible light optimized by genetic algorithms,” OSA Continuum, vol. 1, pp. 796-804, 2018.
- [6] J. P. da Silva, D. S. Bezerra, V. F. Rodriguez-Esquerre, I. E. da Fonseca, and H. E. Hernandez-Figueroa, “Ge-doped defect-core microstructured fiber design by genetic algorithm for residual dispersion compensation,” IEEE Photonics Technology Letters, vol. 22, no. 18, pp. 1337-1339, 2010.
- [7] D. Karaboga and B. Basturk, “A powerful and efficient algorithm for numerical function optimization: Artificial bee colony (ABC) algorithm,” Journal of Global Optimization, vol. 39, no. 3, pp. 459-471, 2007.
- [8] G. N. Malheiros-Silveira and F. G. Delalibera, “Inverse design of photonic structures using an artificial bee colony algorithm,” Applied Optics, vol. 59, pp. 4171-4175, 2020.
- [9] A. Dourado-Sisnando, V. F. Rodríguez-Esquerre, L. F. Viera, and C. E. Rubio-Mercedes, “Inverse design of tapers by bio-inspired algorithms,” Journal of Microwaves, Optoelectronics and Electromagnetic Applications, vol. 19, pp. 39-49, 2020.
- [10] A. Dourado-Sisnando, L. F. Vieira, V. F. Rodríguez-Esquerre, and F. G. S. Silva, “Artificial immune system optimization of complete bandgap of bidimensional anisotropic photonic crystals,” IET Optoelectronics, vol. 9, no. 6, pp. 333-340, 2015.
- [11] D. F. Rêgo, I. L. Gomes de Souza, V. F. Rodriguez-Esquerre, and G. N. Malheiros-Silveira, “Inverse design of broadband absorption in the visible with plasmonic multilayered planar structures,” Photonics, vol. 10, no. 9, pp. 922-1-922-9, 2023.
- [12] D. Correia, V. F. Rodriguez-Esquerre, and H. E. Hernandez-Figueroa, “Genetic-algorithm and finite-element approach to the synthesis of dispersion-flattened fiber,” Microwave and Optical Technology Letters, vol. 31, no. 4, pp. 245-248, 2001.
- [13] P. R. Wiecha, A. Y. Petrov, P. Genevet, and A. Bogdanov, “Inverse design of nanophotonic devices and materials,” Photonics and Nanostructures: Fundamentals and Applications, vol. 52, pp. 101084-1-101084-15, 2022.
- [14] Y. Zheng and Z. Wu, Intelligent Nanotechnology: Merging Nanoscience and Artificial Intelligence (Materials Today) Amsterdam, Netherlands: Elsevier, 2022.
- [15] S. Wu, X. Mu, L. Cheng, S. Mao, and H. Y. Fu, “State-of-the-art and perspectives on silicon waveguide crossings: A review,” Micromachines, vol. 11, no. 3, pp. 326-1-326-18, 2020.
- [16] Y. Ma, Y. Zhang, S. Yang, A. Novack, R. Ding, A. E.-J. Lim, G.-Q. Lo, T. Baehr-Jones, and M. Hochberg, “Ultralow loss single layer submicron silicon waveguide crossing for SOI optical interconnect,” Optics Express, vol. 21, no. 24, pp. 29374-29382, 2013.
- [17] L. Han, X. Ruan, W. Tang, and T. Chu, “Ultralow-loss waveguide crossing for photonic integrated circuits by using inverted tapers,” Optics Express, vol. 30, pp. 6738-6745, 2022.
- [18] K. Goudarzi and M. Lee, “Inverse design of a binary waveguide crossing by the particle swarm optimization algorithm,” Results in Physics, vol. 34, pp. 105268-1-105268-7, 2022.
- [19] C. Valenzuela, B. Crawford, R. Soto, E. Monfroy, and F. Paredes, “A 2-level metaheuristic for the set covering problem,” International Journal of Computers, Communications & Control, vol. 7, no. 2, pp. 377-387, 2012.
- [20] E. G. Talbi, Metaheuristics: From Design to Implementation Hoboken, NJ, USA: Wiley, 2009.
- [21] Y. Elesin, B. S. Lazarov, J. S. Jensen, and O. Sigmund, “Time-domain topology optimization of 3D nanophotonic devices,” Photonics and Nanostructures: Fundamentals and Applications, vol. 12, pp. 23-33, 2014.
- [22] J. P. M. Silva and K. A. Sakallah, “GRASP-A new search algorithm for satisfiability,” International Conference on Computer-Aided Design (ICCAD), pp. 220-227, 1996.
- [23] R. M. Aiex, M. G. C. Resende, P. M. Pardalos, and G. Toraldo, “GRASP with path relinking for three-index assignment,” INFORMS Journal on Computing, vol. 17, no. 2, pp. 224-247, 2005.
- [24] P. Festa and M. G. C. Resende, “GRASP,” Handbook of Heuristics Cham, Switzerland: Springer, pp. 465-488, 2018.
- [25] T. A. Feo and M. G. C. Resende, “Greedy randomized adaptive search procedures,” Journal of Global Optimization, vol. 6, no. 2, pp. 109-133, 1995.
- [26] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, “Optimization by simulated annealing,” Science, vol. 220, no. 4598, pp. 671-680, 1983.
- [27] D. S. Johnson, C. R. Aragon, L. A. McGeoch, and C. Schevon, “Optimization by simulated annealing: An experimental evaluation; part I, graph partitioning,” Operations Research, vol. 37, no. 6, pp. 865-892, 1989.
- [28] T. L. Friesz, H.-J. Cho, N. J. Mehta, R. L. Tobin, and G. Anandalingam, “A simulated annealing approach to the network design problem with variational inequality constraints,” Transportation Science, vol. 26, no. 1, pp. 18-26, 1992.
- [29] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, “Equation of state calculations by fast computing machines,” Journal of Chemical Physics, vol. 21, no. 6, pp. 1087-1092, 1953.
- [30] Y. Zhou, Y. Du, P. Zeng, F. Jiao, J. Wang, J. Song, Y. Tian, M. Yu, and X. Guan,“Compact inverse-designed ultra-broadband and low-loss waveguide crossing for photonic integration,” Optics Letters, vol. 50, pp. 7296-7299, 2025.
- [31] L. He, H. Ji, Y. Dong, and X. Zhang,“An inverse-designed topological waveguide crossing on valley-Hall photonic crystals,” Optics Communications, vol. 587, pp. 131948:1-131948:27, 2025.
- [32] Z. Zhong, Q. Liu, J. Yang, X. Wu, M. Geng, K. Wei, and Z. Zhang, “Design of low-loss and low-crosstalk compact waveguide crossing based on the adjoint method and direct binary search algorithm,” Physics Letters A, vol. 563, pp. 131060:1-131060:15, 2025.
- [33] B. Nan, Y. Ren, R. Wu, L. Song, R. Liu, Y. Zheng, M. Wang, and Y. Cheng, “Near-zero crosstalk and ultra-low loss waveguide crossings enabled by three-dimensional Ta2O5-on-LNOI integrated photonic platform,” arXiv preprint, 2025.
- [34] Y.-J. Chen, K.-C. Chang, T.-L. Yang, and S.-C. Chu, “A hybrid simulated annealing approach for loaded phase optimization in digital lasers for structured light generation,” Photonics, vol. 12, no. 10, pp. 1005:1-1005:12, 2025.
- [35] Physically Aware Wavelength-Routed Optical NoC Design (ASPDAC’25), “Simulated Annealing (SA) optimization,” ASPDAC 2025, Tokyo, Japan, 2025.
- [36] M. Laguna, R. Martí, A. Martínez-Gavara, S. Pérez-Peló, and M. G. C. Resende, “Greedy randomized adaptive search procedures with path relinking: An analytical review of designs and implementations”, European Journal of Operational Research, vol. 327, no. 3, pp. 717-734, 2025.
- [37] T. A. Feo, M. G. C. Resende, “Greedy Randomized Adaptive Search Procedures”, Journal of Global Optimization, v. 6, pp. 109-133, 1995.
-
Editor:
Carlos E. Capovilla
