Preface #

In recent years, inverse design has become an important methodology in the development of photonic devices. Unlike traditional forward design workflows—which rely heavily on experience, intuition, or extensive parameter sweeps—inverse design reformulates device development as an optimization problem with clearly defined performance metrics. By specifying the desired electromagnetic response and allowing the algorithm to freely explore the design space, inverse design can automatically generate high-performance photonic structures that are often difficult to derive manually.

In this case, we demonstrate a Y-branch splitter and how automatic geometric optimization can be achieved using parameterized structural descriptions. The algorithm automatically adjusts the control points of the parameterized structure, enhancing both design efficiency and device performance.

Simulation settings #

The complete simulation workflow in this case is automatically driven by the attached Python script simworks_ybranch_opt3D.py, without the need for manual configuration of the structure or light source in the software interface. The script is built on the SimWorks optimization framework and includes the following key modules:

• Base Environment Setup

  The script first configures the software interface paths and imports necessary Python scientific libraries (e.g., NumPy, SciPy). The paths for os.add_dll_directory and sys.path.append must point to the \bin and \api\python folders under the software installation directory, respectively.

  import os
  import sys
  
  # This is the default installation path. If you install the software in a different directory, please update the path accordingly.
  os.add_dll_directory(r"C:\Program Files\SimWorks\SimWorks FD Solutions\bin")
  sys.path.append(r"C:\Program Files\SimWorks\SimWorks FD Solutions\api\python")
  
  import numpy as np
  import scipy as sp
  import simworks
  from debug_simworksybranch import simwork_ybranch_init
  

  Required libraries can be installed using the following command, If you encounter any issues, please refer to the official documentation for each library (NumPy、SciPy、Matplotlib、PyQt6).

  pip install numpy scipy matplotlib PyQt6
  

  Then, the script calls the simwork_ybranch_init function to load the predefined base simulation model file debug_simworksybranch.py. This model defines the complete simulation environment for the SOI platform and includes all objects required for the optimization process, specifically: a broadband mode source spanning $1.3-1.8\ \mu m$, an FDFP monitor (opt_fields) used for shape optimization, an output power monitor (fom), and a custom mesh for enhanced simulation accuracy.

• Parameterized Geometry Modeling

  The core of inverse design is converting a continuous physical structure into a set of discrete, optimizable parameters. This case uses a “control points + spline interpolation” parameterization: a polygonal structure defined by $20$ control points is placed in the device's splitting region, where the $X$ coordinates are fixed and evenly distributed, while the $Y$ coordinates serve as the variables to be optimized.

  # Define the span and number of points
  initial_points_x = np.linspace(-1.0e-6, 1.0e-6, 10)
  initial_points_y = np.linspace(0.25e-6, 0.6e-6, initial_points_x.size)
  
  
  def splitter(params):
      ''' Defines a taper where the paramaters are the y coordinates of the nodes of a cubic spline. '''
  
      ## Include two set points based on the initial guess. The should attach the optimizeable geometry to the input and output 
      points_x = np.concatenate(
          ([initial_points_x.min() - 0.01e-6], initial_points_x, [initial_points_x.max() + 0.01e-6]))
      points_y = np.concatenate(([initial_points_y.min()], params, [initial_points_y.max()]))
  
      ## Up sample the polygon points for a smoother curve. Some care should be taken with interp1d object. Higher degree fit
      # "cubic", and "quadratic" can vary outside of the footprint of the optimization. The parameters are bounded, but the
      # interpolation points are not. This can be particularly problematic around the set points.
      n_interpolation_points = 20
      polygon_points_x = np.linspace(min(points_x), max(points_x), n_interpolation_points)
      interpolator = sp.interpolate.interp1d(points_x, points_y, kind='cubic')
      polygon_points_y = interpolator(polygon_points_x)
  
      ### Zip coordinates into a list of tuples, reflect and reorder. Need to be passed ordered in a CCW sense 
      polygon_points_up = [(x, y) for x, y in zip(polygon_points_x, polygon_points_y)]
      polygon_points_down = [(x, -y) for x, y in zip(polygon_points_x, polygon_points_y)]
      polygon_points = np.array(polygon_points_up[::-1] + polygon_points_down)
  return polygon_points
  

  During each optimization iteration, the current Y coordinates are used to generate a smooth, continuous boundary curve via cubic spline interpolation. This curve, mirrored about the axis, forms the complete polygon geometry of the splitter.

  

  To ensure manufacturability, each Y-coordinate is constrained within the physical range of $0.2$ to $0.8\ \mu m$. This parameterization balances design freedom (adjustable via the number of control points) and optimization efficiency.

• Optimization Problem Formulation

  A clear optimization goal is defined: maximize the total transmission from the input to the output ports. In the script, this is implemented via the ModeMatch module, with the ideal target set to unity transmission across the wavelength range.

  ######## FIGURE OF MERIT ########
  fom = ModeMatch(monitor_name='fom',
                  mode_number='fundamental mode',
                  direction='Forward',
                  target_T_fwd=lambda wl: np.ones(wl.size),
                  norm_p=1)
  
  ######## OPTIMIZATION ALGORITHM ########
  scaling_factor = 1.0e6
  optimizer = ScipyOptimizers(max_iter=30,
                              method='L-BFGS-B',
                              scaling_factor=scaling_factor,
                              pgtol=1.0e-5,
                              ftol=1.0e-5,
                              scale_initial_gradient_to=0.0)
  

  The optimization algorithm uses the efficient L-BFGS-B method (quasi-Newton), with a maximum of 30 iterations. To accelerate gradient computation, the script employs the adjoint method, which computes the gradient of the objective with respect to all control points using only two electromagnetic simulations (one forward, one adjoint), reducing complexity from O(N) in traditional perturbation methods to O(1) and significantly speeding up optimization.

• Automated Optimization Loop

  Once the modules are integrated, the script runs a fully automated optimization loop:

  1. Geometry Update: Generate a new polygon boundary based on current control points.

  2. Electromagnetic Simulation: Automatically run FDTD simulations in SimWorks.

  3. Performance Evaluation: Calculate transmission at output ports and compare to the target.

  4. Gradient Computation: Compute gradients using the adjoint method.

  5. Parameter Update: Update control point positions based on the optimizer.

  This loop continues until convergence or the maximum iteration number is reached, requiring no manual intervention.

Simulation results #

Before presenting the final optimization results, we first provide several intermediate iterations to illustrate how the structure and performance converge throughout the process. The figure below shows the device geometry and field propagation distribution in the splitter region during the first iteration. At this early stage, the total transmission is relatively low, and significant radiation loss occurs within the branching transition region.

As the optimization iterations progress, the geometry of the splitting region is gradually refined, resulting in a significantly smoother optical field transition. Upon completion of the optimization, the program automatically generates and saves the final 3D project file, along with two optimization process records. The file optimization_report.txt logs the FOM value and optimized parameters for each generation, while convergence_report.txt stores the complete iteration history data for convergence analysis.

The final generated 3D project file is shown in the figure below and can be directly opened in the SimWorks software for subsequent analysis:

The final optimized structure demonstrates outstanding performance across the operational band: the total transmission efficiency at the output increases from an initial value of approximately $0.65$ to about $0.92$, with a significant reduction in losses. This result fully showcases the powerful capability of inverse design in generating high-performance photonic device geometries.