Gradient Probabilistic Algorithm for Compact Lithium Niobate Integrated Photonic Devices

Abstract

Compact photonic devices are highly desired in photonic integrated circuits. In this work, we use an efficient inverse design method to design a 50/50 beam splitter in lithium niobate integrated platforms. We employ the Gradient Probability Algorithm (GPA), which is built upon traditional gradient algorithms. The GPA utilizes the adjoint method for the comprehensive calculation of the electric field across the entire design area in a single iteration, thereby deriving the gradient of the design area. This enhancement significantly accelerates the algorithm’s execution speed. The simulation results show that an ultracompact beam splitter with a footprint of $13\mathsf{\mu }$m × $4.5\mathsf{\mu }$m can be achieved in lithium niobate integrated platforms, where the insertion loss falls below 0.5 dB within the 1500 nm to 1700 nm range, thus reaching its lowest point of 0.15 dB at 1550 nm.

1. Introduction

With the rapid development of advanced communication technologies such as 5G, the Internet of Things (IOT), quantum communication, satellite internet, and smart city infrastructure, the demand for efficient data transmission, low latency, high security, and broad coverage has become increasingly urgent. These technologies require communication systems to have exceptional transmission efficiency, stability, and reliability to support the rapid transfer and processing of vast amounts of data [1]. Lithium niobate materials provide crucial support in this regard, with their superior electro-optic properties, high optical transparency, and strong piezoelectric effect, thus making them an ideal choice for high-performance optical devices. To address the need for efficient LN photonic devices, there is a pressing demand for high-quality lithium niobate thin film materials. The breakthrough in this direction came in 1998 when Levy et al. pioneered the fabrication of LN thin films using the “ion slicing” method [2]. In recent years, advancements in technology have led to the development of single-crystal LN on insulator (LNOI) through ion implantation and bonding techniques, thereby providing a promising platform for integrated optics. Next, high-performance devices such as 50/50 LN beam splitters will further enhance the capabilities of communication systems, thus providing stronger support for the future development of communication technologies. This splitter can equally divide incoming optical signals into two parts, with each carrying half of the original power, thus making it indispensable in building optical networks, signal processing, and various sensing systems. The 50/50 beam splitter not only improves the performance of existing communication systems but also supports the construction of more complex and efficient optical communication networks in the future. As technology continues to advance, the 50/50 beam splitter will play an increasingly important role in high-speed data transmission, optical signal processing, and integrated optical circuits, thus driving modern optical communication technology towards greater efficiency and intelligence.

Presently, the advancement of 50/50 beam splitters on the LNOI platform is in progress. Qiyuan Yi et al. designed an ultrawide band 50/50 splitter using lithium niobate thin films, thus achieving an 800 nm bandwidth. The device exhibits insertion losses of 0.2 dB, 0.16 dB, and 0.53 dB at wavelengths of 1310 nm, 1550 nm, and 2000 nm, respectively [3]. Meanwhile, D. Li et al. developed a high-performance asymmetric multimode interferometer (MMI) splitter on an X-cut LNOI platform. This MMI splitter can achieve splitting ratios from 50:50 to 95:5. The device size measures at $5.8\mathsf{\mu }\mathrm{m}\times (26.5\sim 35.6\mathsf{\mu }\mathrm{m})$, with insertion losses ranging from 0.1 dB to 0.9 dB [4].

Furthermore, compact photonic devices play a pivotal role in achieving high integration and low power consumption for photonic applications [5,6,7,8,9,10]. Specifically, their significance lies in their ability to accommodate a greater number of functional components within a small footprint, thereby enhancing the integration and functionality of photonic integrated circuits. The small size results in shorter signal transmission paths, thus leading to reduced energy dissipation and faster response speeds. These characteristics are particularly advantageous for high-speed optical communication and data processing.

However, LN photonic devices face challenges in dense integration, since their relatively low refractive index cannot make strong light confinement and light control within small footprints [11].

The emergence of optical inverse design has improved this situation; optical inverse design, as a methodology for crafting photonic devices and leveraging calculation techniques and optimization algorithms, enables precise tailoring over the structure and parameters of photonic devices, which facilitates the design of ultracompact devices [12]. Furthermore, inverse design offers greater flexibility in adjusting the structure with various creative degree of freedoms. In the realm of inverse design, there are many prevalent optimization algorithms, including the direct binary search (DBS), genetic algorithm (GA), particle swarm optimization (PSO), and gradient probability algorithm. The direct binary search algorithm is a robust yet straightforward approach with broad applicability, thus often meeting various design requirements effectively [13]. Its thorough exploration of the entire parameter space, however, renders it sensitive to the number of parameters in the design space. When more parameters are involved, this algorithm will become much more time-consuming. The genetic algorithm and particle swarm optimization are heuristic methods simulating populations where prominent individuals are retained. Despite their efficacy, these algorithms entail substantial computational demands [14], especially in scenarios with intricate design parameters, thus potentially diminishing efficiency.

The gradient probability algorithm, leveraging the adjoint method, conducts forward and inverse simulations of the design area to derive the forward and inverse electric fields, thus subsequently calculating gradient information. Remarkably, this method obtains gradient information for the entire design area through only two simulations, thus making it independent of the parameter complexity in the design area. In comparison to both violent and heuristic algorithms, the gradient probability algorithm not only ensures optimal solutions but also boasts faster computational speeds [15].

To achieve a more compact device design, J. Xu et al. utilized the DBS algorithm to design a dual-mode 50/50 beam splitter, thus achieving an insertion loss of 0.83 dB over a wavelength range from 1588 nm to 2033 nm [16]. Although the DBS algorithm provided high bandwidth and low loss, it required significant time, with the entire optimization process taking 72 h. To save on time costs, K. Wang et al. employed the digital adjoint method to design a single-mode 3 dB power splitter with a footprint of $2.6\mathsf{\mu }\mathrm{m}\times 2.6\mathsf{\mu }\mathrm{m}$ on a 220 nm thick top silicon layer of the silicon-on-insulator (SOI) platform. Within a 40 nm bandwidth range (1530–1570 nm), the insertion loss varied between 0.33 dB and 0.3 dB [17]. By using the adjoint method, a significant amount of simulation time was saved, with the optimization process taking 1.2 h. To further reduce loss while optimizing time efficiency, we used the GPA algorithm to design a 50/50 beam splitter on the LN integrated platform [18]. Our algorithm, based on the adjoint method, differs in that it further processes the gradient information. Finally, we achieved an X-cut LN ultracompact beam splitter with dimensions of $13\mathsf{\mu }\mathrm{m}\times 4.5\mathsf{\mu }\mathrm{m}$ on the LN integrated platform. The optimization process took 10.5 h. Across a wavelength range from 1500 nm to 1700 nm, the insertion loss remained below 0.5 dB and reached 0.15 dB at 1550 nm.

2. Gradient Probability Algorithm

Figure 1 shows the schematic of the beam splitter structure. The design in based on a thin-film LN integrated platform. Its etching depth is $h=200$ nm. We took a rectangular region for the inverse design with a length of $M=13\mathsf{\mu }\mathrm{m}$ and a width of $N=4.5\mathsf{\mu }\mathrm{m}$. The input and output waveguides have a width of $w=0.8\mathsf{\mu }\mathrm{m}$, and the spacing of the two output waveguides is $c=1.6\mathsf{\mu }\mathrm{m}$. The design area is partitioned into $i\times j=65\times 22$ square pixels, with each possessing identical geometric dimensions. The individual pixel length is b = 200 nm, with the central circular region representing the etched hole. Considering the fabrication practice, the sidewalls of the waveguide, as well as the holes, have an angle of $\theta ={80}^{\circ }$. In geometric structures with sharp edges or corners, electric or magnetic fields tend to concentrate at these points, thus creating singular fields [19]. To effectively mitigate the impact of these singular fields, dielectric caps [20] can be used to smooth the field distribution. Alternatively, employing smooth transition curves and rounded edges can help avoid sharp edges and corners. Additionally, gradient material techniques can be utilized to gradually change the dielectric constant or conductivity at material interfaces, thus smoothing the field distribution. In this work, due to manufacturing and design considerations, numerical techniques were employed to reduce the impact of singular fields. High-precision mesh refinement and appropriate boundary conditions were used to minimize numerical instability. The boundary conditions were set to a Perfectly Matched Layer (PML), and fixed-size locally refined meshes were applied within the design area to enhance the accuracy of the numerical simulations. At the sharp edges of the geometric structure, the mesh size was set to 20 nm to capture the rapid variations in the electromagnetic field. Compared to larger meshes in other regions, this local mesh refinement significantly smooths the field distribution, reduces field concentration, and effectively mitigates the numerical instability caused by singular fields, thereby improving the stability and accuracy of the simulation results.

Within inverse design methodologies, devices are typically segmented into pixels, thus aiding computer processing. In prevalent algorithms, pixels generally exhibit two states: 0 and 1. Here, 0 signifies no etching, while 1 signifies the etched state. In this study, the gradient probability search algorithm employs more than two states for pixels. Specifically, we utilized three pixel states—type-I, type-II, and type-III—denoting a different etching radius, as illustrated in Figure 1. Type I, II, and III correspond to hole diameter of $D_1=200$ nm, $D_2=150$ nm, and $D_3=0$ nm. Given that the purpose of the 50/50 beam splitter is to evenly divide the input beam into two beams of equal intensity, we employed a symmetrical etching method to accomplish this objective. This method entails etching the upper design area and simultaneously replicating the same etching process at the corresponding location in the lower design area, thus ensuring symmetry between the upper and lower design areas. After the above processing, both output terminals will have the same output, and only half of the parameters need to be simulated, which is $i\times j=65\times 11=715$ pixels. $Q(i,j)$ is utilized to denote the pixel parameters of row i and column j, where $i=$ 1, 2, 3, …, 65 and $j=$ 1, 2, 3, …, 11. Figure 2a demonstrates the pixel values after initialization.

Due to the implementation of symmetrical design, the transmission efficiency of one of the outputs was utilized as the figure of merit (FOM), which is given by

$$FOM=T_{out}.$$

(1)

The adjoint method proves instrumental in calculating gradients. The gradient of the FOM concerning the spatial dielectric constant $e\left(x,y,z\right)$ within the design area can be expressed as

$$\frac{dFOM}{d{\epsilon }_r}=R=real\left(E_{input-fore}\ast E_{ouput-back}\right),$$

(2)

where ${\epsilon }_r$ represents the spatial dielectric constant of the pixel, $E_{input-fore}$ denotes the forward electric field of the light source, and $E_{output-back}$ corresponds to the accompanying electric field. Each pixel is uniformly subdivided into a $10\times 10$ grid, and the electric field information of each pixel is acquired through monitoring. The electric field information is sampled at the grid intersections. Consequently, the electric field point information matrix for each pixel is $11\times 11$, and the collective electric field point information matrix amounts to $651\times 111$.

After employing the adjoint method for gradient calculation, the overall gradient value R for the design area is derived. Subsequently, R is uniformly partitioned into a $65\times 11$ matrix, followed by aggregation to compute the average value. Figure 2b illustrates the gradient values corresponding to each pixel after initialization. This process yields the gradient $r\left(i,j\right)$ corresponding to each pixel, thus elucidating the overall refractive index change trend of the pixel. Following the adjoint method, an initial step involves conducting a forward simulation to acquire the forward electric field. Subsequently, only the inverse light source is activated to initiate the inverse simulation, thereby leading to the acquisition of the adjoint electric field [21]. In accordance with Equation (2), the real part obtained by multiplying the forward electric field and the inverse electric field yields the gradient information for the entire design area. Consequently, the gradient of the FOM concerning each pixel can be expressed as

$$dp\left(i,j\right)=\frac{dFom}{dQ\left(i,j\right)}=Average\left(E_p\left(i,j\right)\right).$$

(3)

where $E_p\left(i,j\right)$ is the product of the forward and backward electric fields of each pixel. Subsequently, the gradient value of each pixel is mapped into the probability space. The absolute value of the gradient for each pixel is then input into the hyperbolic tangent function $G=tanh\left(x\right)$ as the abscissa, thus resulting in the normalized gradient probability G, which is expressed as

$$G=\alpha \times \frac{\left[exp\left(x\right)-exp\left(-x\right)\right]}{\left[exp\left(x\right)+exp\left(-x\right)\right]},$$

(4)

$$x=\left|\nabla \frac{dFom}{dQ\left(i,j\right)}\right|.$$

(5)

a represents the learning rate. We control the number of pixels to be selected later by changing the value of a. The value of a is determined by repeated adjustments The normalized gradient probability for the entire design area is partitioned for each pixel, thus resulting in the corresponding gradient probability for each pixel denoted as $g\left(i,j\right)$. Figure 2c illustrates the normalized gradient probability corresponding to each pixel.

Following the algorithmic procedures described above, we achieve the computation of the normalized gradient probability. Subsequently, the gradient probability associated with each pixel is calculated. By incorporating this probability alongside the original gradient value of the pixel, we proceed to adjust the parameter value $Q\left(i,j\right)$ based on the modification rule stipulated, which is given by

$$Q\left(i,j\right)=\left\{\begin{array}{c}if\frac{df}{dQ\left(i,j\right)}>0,Q\left(i,j\right)=D_{K+1}\\ if\frac{df}{dQ\left(i,j\right)}<0,Q\left(i,j\right)=D_{K-1}.\end{array}\right.$$

(6)

The specific steps are outlined as follows: (1) Generate a random number L, where L ranges between 0 and 1. (2) Compare $g\left(i,j\right)$ with L, and select pixels $Q\left(i,j\right)$, where $g\left(i,j\right)>L$. (3) Modify these selected pixels based on the gradient value of $Q\left(i,j\right)$. If the gradient value of $Q\left(i,j\right)$ is positive, increase the parameter: $D_k=D_{k+1}$, where $K\le 3$. Conversely, if the gradient value of $Q\left(i,j\right)$ is negative, decrease the parameter: $D_k=D_{k-1}$, where $K\ge 1$. Here, a represents the learning rate, which is set to $a=10$. The number of selected pixels is controlled by a to be 10% of the total, thereby mitigating the risk of falling into local optima. Figure 2d displays the selected pixels earmarked for modification, which are accompanied by the respective gradient values. Following the prescribed modification rules outlined above, Figure 2e presents the pixel values after modification.

The initial iteration process of the algorithm is outlined as follows:

Step 1: Initialize the divided pixels with random integers ranging from one to three, thus forming a 65 × 11 initialization matrix, which is expressed as $Q\left(i,j\right)=D_k,i\in \left(1,65\right),j\in \left(1,11\right),D_k\sim random\left(1,3\right),k=1,2,3$.

Step 2: Calculate the FOM.

Step 3: Utilizing the adjoint method, calculate the electric field information for the entire design area, and compute the gradient information for the entire design area using Equation (3). Subsequently, calculate the gradient value $r\left(i,j\right)$ for each pixel. Figure 3a,b illustrate the gradient information R for the entire design area and the gradient value $r\left(i,j\right)$ for each pixel, respectively.

Step 4: Derive the normalized gradient probability $g\left(i,j\right)$ for each pixel based on Equations (4) and (5). Compare the gradient probability of each pixel with a random number L, and select pixels where $g\left(i,j\right)>L$ for modification. Figure 3c depicts the normalized gradient probability $g\left(i,j\right)$.

Step 5: Identify the pixels to be modified, and adjust their values according to the positive and negative initial gradient values $r\left(i,j\right)$ of these pixels. The modified pixel value distribution state can be obtained as per the modification rule outlined in Euqation (6). With this, the initial iteration concludes, and the updated pixel distribution values are utilized to proceed to the subsequent iteration.

The sampling methodology embedded in our algorithm manifests itself with discernible advantages, thus underscoring its efficacy. Capitalizing on gradient information throughout the entire spatial domain, our algorithm outperforms traditional methodologies, thereby yielding more potent enhancements in device performance. In stark contrast to the deterministic outcomes generated by conventional gradient algorithms under specified initial states and parameters, our algorithm adopts a probabilistic sampling strategy that injects an element of unpredictability. While this momentary uncertainty might entail a transitory dip in the FOM value, it plays a pivotal role in facilitating the algorithm’s evasion of local optima. This distinctive feature significantly augments stability, especially when confronted with diverse initial conditions. Next, we will verify this view through simulation.

This section has employed an iterative approach to elucidate the operational principles and optimization intricacies of the gradient probability algorithm. The next section will unveil the ultimate optimization outcomes, thus conducting a comprehensive analysis of the algorithm’s performance through a comparison with traditional reverse design algorithms.

The above algorithm has been implemented together with the Lumerical 3D-FDTD solution [22]. Throughout each iteration, the optimization interacts with the FDTD solution to retrieve the electromagnetic field in space and subsequently obtain gradient information. Following each iteration, the FDTD updates the device structure [23].

3. Simulation and Results

The algorithm stability plays a pivotal role in ensuring the generation of dependable and uniform outcomes across varying iteration counts or conditions. A stable algorithm exhibits resilience to alterations in input data, initial conditions, or parameters, thereby ensuring consistent outputs in diverse scenarios. This stability not only bolsters the reproducibility of research but also enables the algorithm to maintain steadfast performance in different experiments or application contexts. This enhances the reliability and applicability of the algorithm in scientific research and engineering applications. To assess the algorithm stability, we conducted three optimizations, each utilizing distinct initial structures, with the iteration count set to 100 for each optimization.

This approach ensures consistent algorithmic performance across a specific range of datasets [24]. In accordance with the aforementioned algorithmic execution steps, Figure 4 shows the ultimate device structure, light field distribution, and performance subsequent to three rounds of optimization [25]. It illustrates the variation in FOM values throughout the optimization process [26].

In the FOM evolution plot presented in Figure 4c, the variations in FOM values across three distinct optimization processes are depicted by the black, red, and blue lines. The plot unveils a noteworthy pattern where, during the initial optimization stages, the FOM experienced substantial growth, thus culminating in its peak after surpassing 50 iterations. Following this, the rate of the FOM increase diminished, and at around 70 iterations, the FOM achieved its highest values for the three optimization iterations: 0.4798, 0.4831, and 0.4783, with corresponding iteration counts of 68, 70, and 69, respectively. Subsequent to reaching these peak values, the FOM levels stabilized. The results from three optimization iterations underscore the algorithm’s consistent performance under diverse initial conditions, thus affirming its robust stability [27].

Figure 5 shows the changes in the insertion loss for three optimizations at different wavelengths. It is observed that within the wavelength range from 1260 nm to 1380 nm, the insertion loss exceeded 1 dB, while in the range from 1390 nm to 1490 nm, it surpassed 0.5 dB. Between 1500 nm and 1700 nm, the insertion loss remained below 0.5 dB, thus reaching its minimum value at 1550 nm. In the context of the three distinct initial structural optimizations, the corresponding minimum values of the insertion loss were 0.179 dB, 0.15 dB, and 0.192 dB.

Algorithms need to be considered not only for their accuracy but also for their operational efficiency. An efficient algorithm can save a lot of time for the designer. To verify the algorithm’s speed, we compared it with the traditional DBS algorithm, which uses the same design area and initial conditions as the probabilistic gradient method and takes the output power as the optimization objective. This algorithm has only two parameter states: type-I for an etching diameter of 200 nm and type-II for no etching. The optimized structure of the DBS algorithm is shown in Figure 6:

In Figure 6c, it is evident that the FOM gradually increased with the growing number of iterations, thus reaching its peak at the 1414th iteration with a recorded value of FOM = 0.4842. Subsequently, the FOM stabilized.

Figure 7 illustrates the insertion loss at various wavelengths. Notably, within the wavelength range of 1260 nm to 1380 nm, the insertion loss exceeded 0.5 dB. In the intervals from 1390 nm to 1520 nm and from 1620 nm to 1700 nm, the insertion loss surpassed 0.2 dB. Between 1530 nm and 1610 nm, the insertion loss remained below 0.2 dB, thus reaching its minimum value of 0.13 dB at 1550 nm.

Upon comparing

Table 1. Comparison of optimization algorithms for 50/50 beam splitter.

Method	Number of Iterations	Simulation Time	Insertion Loss
DBS	1430	48 h	0.13 dB
GPA
(First optimization)	100	10.7 h	0.18 dB
GPA
(Second optimization)	100	10.5 h	0.15 dB
GPA
(Third optimization)	100	10.6 h	0.19 dB

, it is evident that the final FOM values of the grobability gradient algorithm and the DBS algorithm are very close. However, the probability gradient algorithm exhibited significantly faster speed, which is primarily attributed to its sampling method within the design area. By leveraging gradient information to guide the sampling process, this algorithm effectively reduces computational complexity. In contrast, the DBS algorithm gradually approached the optimal solution by continuously adjusting parameters, thus necessitating more simulation time and resulting in a slower algorithm speed [28]. Consequently, the probability gradient algorithm demonstrates higher efficiency and performance in practical applications than DBS. Through such optimization algorithms, quicker optimization results can be obtained, thereby allowing for more efficient decision making in the design process [29].

4. Discussion

Comparing our photonic device with those previously reported is of paramount importance. This comparative analysis allows us to effectively evaluate the performance and advantages of our novel intelligent algorithm design. By juxtaposing our device against existing ones, we can objectively assess performance metrics such as efficiency and stability. Specifically, to demonstrate the efficacy of our device engineered through the gradient probability algorithm, we compared it with the reported lithium niobate 50/50 beam splitter, as shown in

Table 2. Comparison table with reported lithium niobate 50/50 beam splitters.

Design Methods	Insertion Loss	Device Length	Reference
MMI	0.02 dB	$223\mathsf{\mu }\mathrm{m}$	[30]
Y-branch	0.3 dB	$118\mathsf{\mu }\mathrm{m}$	[31]
Adiabatic tapers	0.2, 0.16, 0.53 dB	$60\mathsf{\mu }\mathrm{m}$	[3]
GPA	0.15 dB	$13\mathsf{\mu }\mathrm{m}$	this work

.

As shown in Table 2, the device utilizing the MMI structure exhibited an insertion loss of 0.02 dB, but it also has a relatively large size, with a length of $223\mathsf{\mu }\mathrm{m}$. In contrast, the device using the Y-branch structure has a designed length of $118\mathsf{\mu }\mathrm{m}$ and yielded an insertion loss of 0.3 dB. The beam splitter designed with adiabatic tapers demonstrated losses of 0.2 dB, 0.16 dB, and 0.53 dB at wavelengths of 1310 nm, 1550 nm, and 2000 nm, respectively, with a device length of $60\mathsf{\mu }\mathrm{m}$. In comparison, our design shows significant improvements, thus featuring a much smaller size with a length of just $13\mathsf{\mu }\mathrm{m}$ and a lower insertion loss, as detailed in Table 2. These results highlight the efficiency and compactness of our approach, thereby making it highly suitable for advanced photonic integrated circuits.

5. Conclusions

This study employed a gradient-based probability gradient algorithm to design an ultracompact thin-film lithium niobate 50/50 beam splitter, with a small footprint of $13\mathsf{\mu }\mathrm{m}\times 4.5\mathsf{\mu }\mathrm{m}$. We demonstrated the compactness of the GPA inverse-designed 50/50 beam splitter by comparing it with previously reported LN 50/50 beam splitters. The stability of the design was validated through varying iteration times. Across three experiments, the FOM values exhibited a consistent trend of progressively increasing and stabilizing after reaching the maximum value, thereby indicating algorithm convergence. In the range from 1500 nm to 1700 nm, the insertion loss was below 0.5 dB, thus reaching its minimum value of 0.15 dB at 1550 nm. Furthermore, we conducted a comparison with traditional DBS algorithms using identical design parameters. The results indicate similar final effects between the two methods, with the gradient probability algorithm demonstrating faster performance.