Efficient topology optimization of optical waveguide devices utilizing semi-vectorial finite-difference beam propagation method

: An eﬃcient gradient-based topology optimal design approach using beam propagation method (BPM) for a 3-D semi-vectorial design problem of optical waveguide devices is proposed. A semi-vectorial ﬁnite-diﬀerence BPM (SVFD-BPM) based on an alternating direction implicit method (ADIM) is employed for the wave propagation analysis in our design approach. In comparison to conventional 3-D topology optimal design approaches, it is expected that our approach can reduce its computational cost and design an optical waveguide which is too long to analyze using a ﬁnite element method or ﬁnite-diﬀerence time/frequency-domain method. In this paper, we design several waveguide devices (S-bend, power splitter, and mode converter), and it is shown that our topology optimal design approach surely works.


Introduction
A topology optimal design method has great potentials to achieve further miniaturization of optical waveguide devices and to improve device performance.The topology optimization has much design freedom than a sizing optimization or a shape optimization, and it can optimize the refractive index distribution itself.Thus, it is expected that optical devices with higher performance can be obtained by using the topology optimization.The applications of the topology optimal design to optical devices have been reported and the effectiveness has been demonstrated in previous studies [1][2][3][4][5][6][7][8][9][10][11].
Topology optimal design methods reported so far employ a finite element method (FEM) [1][2][3][4][5][6][7][8], a finite-difference time-domain (FDTD) method [9,10], or a finite-difference frequency-domain (FDFD) method [11] for wave propagation analysis.A beam propagation method (BPM) [12][13][14][15] is widely used as an efficient numerical analysis method of optical waveguide devices.An envelope of the unknown electromagnetic field varies slowly for the propagation direction in waveguides without backward reflection.Thus the BPM realizes the cost effective numerical analysis because larger step size is allowed in the propagation direction.It is expected that the devices whose device length is too long to analyze characteristics with FEM, FDTD method, or FDFD method, can be designed by employing the BPM.We have already reported the topology optimal design approach employing the BPM in our previous study [16].However, this approach is applied only for 2-D design problem.3-D design is necessary to consider vertical confinement and radiations to a substrate and a cover.
In this paper, our design approach using the BPM is extended for the case of a semi-vectorial 3-D design problem.We employ a semi-vectorial finite-difference BPM (SVFD-BPM) based on an alternating direction implicit method (ADIM) [14,15] for efficient calculation.Sensitivity analysis is necessary in gradient-based topology optimization, then we employ an adjoint variable method (AVM) as the sensitivity analysis method.Also, we employ a density method which is a way to represent refractive index distribution in the design region.In the density method, the design parameter called a normalized density parameter is allocated to each of discretized domains, and the refractive index in the domain is determined depending on the value of the normalized density parameter.The AVM in the case of using the 2-D scalar or 3-D full-vectorial BPM for a sizing optimization has already been reported in [17,18].We extend the AVM reported in [17,18] for the case of using the density method and the SVFD-BPM based on the ADIM.
In section 2, we begin with a 3-D scalar Helmholtz equation, then we review the basic equations for SVFD-BPM based on the ADIM.After that, we describe the sensitivity analysis method ∂ ∂y in the case of using the SVFD-BPM and the density method.In section 3, we design several weakly guiding waveguide devices to verify the validity of the sensitivity analysis and our design approach.Section 4 is the conclusion.

Formulation of the SVFD-BPM based on ADIM
We assume core and cladding materials are linear, lossless, isotropic, and frequency-independent.From Maxwell's equations in frequency domain, the Helmholtz equation of light in the waveguide shown in Fig. 1 can be obtained, where k 0 is a free space wavenumber, n(x, y, z) is a distribution of refractive index, Φ(x, y, z) is a scalar electric or magnetic field.E x or H y is used for x-polarization and E y or H x is for y-polarization.The definitions of the operators, D xx and D yy , are shown in Table 1.
In the BPM, Φ is expressed in the form [12], where n 0 is a reference refractive index.By substituting Eq. (2) into Eq.( 1) and employing Fresnel approximation (∂ 2 φ/∂z 2 = 0), the basic equation of the BPM is obtained, where By discretizing z-direction using Crank-Nicolson scheme and employing ADIM, the update equations of the BPM can be expressed in following form [15]: where m is a step number in z-direction, ∆z is the step size, and φ m denotes φ(x, y, m∆z).The matrix-vector form of Eqs. ( 5) and ( 6) is obtained by descritizing the transverse direction as follows: where [Γ i ](i = 1, 2, 3, 4) can be expressed by square matrices.In this paper, we employ Stern's method [19] to discretize the transverse directions.

Sensitivity analysis
In the density method, the refractive index distribution in the design region is determined by the distribution of a normalized density parameter ρ(x, y, z).We assume uniform core height in the design region, thus in this case, n and ρ do not depend on y in the design region.The refractive index distribution in the design region is represented by using the density method as follows: where k is a step number in x-direction, n 1 and n 2 are core and cladding refractive indices respectively, n k,m and ρ k,m denote n(k∆x, m∆z) and ρ(k∆x, m∆z) respectively (∆x is a step size in x-direction), and H(ρ) is a modified Heaviside function defined in [16].The gray region, which has an intermediate refractive index between n 1 and n 2 , can be controlled by a penalty parameter.
The refractive index distribution is binarized by making the penalty parameter infinite.In the design examples shown in this paper, the penalty parameter is taken to be a lower value in the initial phase of the optimization, then the value of the parameter is increased as the optimization process is iterated.The penalty parameter is eventually large enough, thus the effect of the binarization is of no significant.The characteristics of optical waveguide devices can be usually expressed by S-parameters.A normalized output power into port-n is represented by following overlap integral in semivectorial problems, where N z is a step number at an output port in z-direction, * denotes complex conjugate, and ψ n is an eigenmode field in port-n.
This overlap integral can be expressed in following matrix-vector form by using rectangle rule.
where † denotes Hermitian transpose, ∆y is a step size in y-direction.From Eqs. ( 7) and ( 8), Eq. ( 11) can be rewritten as where T denotes transpose, and λ n,m+1 T is defined by By differentiating Eq. ( 12) with respect to ρ k,m , we get the following relation: When m N z − 1, the vectors λ n,m and λ n,m+1/2 can be calculated as follows: If m = N z − 1, Eq. ( 15) is replaced with The derivative of the matrices [Γ i ] (i = 1, 2, 3, 4) can be calculated analytically because the refractive index is represented by Eq. ( 9).
In weakly guiding waveguides, the S-parameters can be expressed approximately by only one component. where and φ n is an eigenmode field in port-n.In this case, ψ n is replaced with pφ n in Eq. (11).
Since the devices treated in this paper are weakly guiding waveguides, we employ Eq. ( 18) to evaluate the S-parameters.

Design examples
We show the validity of our design approach and the sensitivity analysis by designing several weakly guiding waveguide devices.
In the design examples, the density parameters are updated every iteration by a steepest descent method.

S-bend waveguide
The design model of an S-bend waveguide is shown in Fig. 2. The structural parameters shown in Fig. 2 are the following: n 1 = 1.45, n 2 = n 3 = 1.445, w = 8 µm, h = 4 µm, l = 100 µm, d 1 = d 2 = 30 µm, s = 15 µm, L x = 80 µm, and L z = 900 µm.In order to maximize the normalized output power in port 2 when the fundamental TE-like or TM-like wave with the wavelength of 1.55 µm is launched into port 1, the objective function is set as follows: The transverse and longitudinal step sizes are taken to be ∆x = ∆y = 0.2 µm and ∆z = 1 µm, respectively.n 0 is taken to be the effective index of the fundamental TE-like or TM-like mode in port 1.
First, we show that our design approach surely works regardless of which field component is used as φ, and the BPM is available in a very low-contrast-index waveguide even though a waveguide structure is complicated to some extent.In this optimal design, a uniform medium (ρ = 0.3 in the design region) is given as an initial structure.In the topology optimization using the density method, fine structures typified by a checker-board pattern often emerge.Such fine structures can be removed by a structural smoothing filter (the detail is described in our previous paper [16]).In this design example, however, the structural smoothing filter is not used in order to show the results without an influence of the smoothing filter.
The normalized output power as a function of the iteration number is shown in Fig. 3.According to Fig. 3, we can see that the normalized output power is improved every iteration, and the power converges in several iterations in the case that φ = E x , H y , E y , or H x .Moreover, an inset in Fig. 3 shows the output power converges quickly.The convergence speed may be strongly affected by the property of an objective function.The composition of the objective function shown in this paper is relatively simple, thus the normalized output power seems to converge relatively quickly.For instance, in the design of an optical triplexer, the convergence speed may be lower because composition of the objective function is more complicated [8].The objective function is expressed as sum of objectives at three different wavelengths when the triplexer is optimized.Since the structural variation to improve the device performance at a certain wavelength may degrade the performance at another wavelength, the objective function may have complicated dependency on the design variables.The binarized optimized structures and the propagation waves are shown in Fig. 4. The term "binarized" in this paper means that ρ is rounded off to the closest whole number, then the refractive index in the design region is determined to be n 1  (core) or n 2 (cladding).The normalized output power in the binarized optimized structure in the case that φ = E x , H y (TE-like wave), E y , or H x (TM-like wave) are 0.929, 0.930, 0.924, 0.929, respectively.It is noted that the binarized optimized structures shown in Fig. 4 are almost similar, and the difference of the output properties is not significant even though mosaic-like structures emerge slightly.The normalized power |S n1 | 2 is calculated using Eq. ( 18) in this paper.If the approximation in the formulation of the BPM is not valid for a specified waveguide, the output power is largely different depending on which field component is used as an unknown variable (φ) even though the polarization is the same [20].Therefore, the numerical results indicate that the BPM is available when refractive index difference is very low (about 0.34%) even though waveguide structure is complicated to some extent.
After this design example, the smoothing filter is employed every iteration in order to avoid a very complex profile and obtain a more practical structure.The smoothing filter itself has been widely used in the topology optimization.Topology optimization using the smoothing filter may lead to a local optimum structure that is simple and has good property to some extent.Although the smoothing filter certainly hinders optimization, considering practical manufacturing, it is necessary to employ this filter.
Next, we will show the dependence of optimization results on an initial structure.A uniform medium (ρ = 0.3 in the design region) or a tilted waveguide (ρ = 0.8 in core region, 0.3 in cladding region) is given as an initial structure.
Figure 5 shows that the normalized output power as a function of the iteration number in the case that a uniform medium or a tilted waveguide is given as an initial structure.In Fig. 5, although there are ripples due to the smoothing filter, the output power is improved almost every iteration.The binarized optimized structures are shown in Fig. 6.Although the profiles of the optimized structures are not identical, the operation principles are probably the same.According to Figs. 6(a) and (c), we can see that offsets emerge in the vicinity of the input and output ports in order to reduce mode mismatching.Moreover, in the design region, the waveguide widths are wider than those of the input and output ports to enhance the confinement of light.The normalized output power in the binarized optimized structure shown in Figs.6(a) and (c) are 0.947 and 0.941, respectively.Certain optimized structures based on the same principle of the operation are automatically obtained whether a special structure is given as an initial structure, or not.In the binarized initial tilted waveguide, the output power is 0.769.Thus the performance of the optimized S-bend waveguides is superior to the binarized initial structures.

Y-branch type power splitter
In this design, a = b = 0.5 and a symmetric condition is imposed to prevent imbalance of power.
The dashed line shown in Fig. 7 represents the symmetric axis.The transverse and longitudinal step sizes are taken to be ∆x = ∆y = 0.2 µm and ∆z = 1 µm, respectively.n 0 is taken to be the effective index of the fundamental TE-like mode in port 1.A uniform medium (ρ = 0.3 in the design region) or a straight Y-branch waveguide (ρ = 0.8 in core region, 0.3 in cladding region) is given as an initial structure, then we will show that the dependence of optimization results on an initial structure.
Figure 8 shows that the normalized output power as a function of the iteration number in the case that a uniform medium or a straight Y-branch waveguide is given as an initial structure.Only the normalized output power in port 2 is shown in Fig. 8 because there is no imbalance of the output power.According to Fig. 8, we can see that the normalized power is improved almost every iteration to realize the 1:1 power splitter in both cases.The binarized optimized structures and the propagation wave are shown in Fig. 9.Although the optimized structures shown in Figs.9(a) and (c) differ, the offset waveguide junctions appear in the vicinity of the output ports and Y-branch structures which can reduce radiation loss into forward direction emerge.A shape of Y-branch like the structures shown in Figs.9(a) and (c) has been presented, e.g. in [21].Therefore, the operation principle of both structures seems to be the same.The normalized output power in the binarized structure shown in Figs.9(a) and (c) are 0.471 and 0.478, respectively.The output power in the binarized initial straight Y-branch waveguide is 0.413, thus our design approach provides the superior structure compared to the initial structures.

Mode converter
Figure 10 shows the design model of a mode converter.This optimal design aims to obtain a converter which converts the 1st higher order TE-like (TE 1 ) wave into the fundamental TE-like (TE 0 ) wave.The structural parameters shown in Fig. 10 where S 21(TE1→TE0) 2 is a normalized output power of the TE 0 wave in port 2 when the TE 1 wave is launched into port 1.The transverse step sizes are taken to be ∆x = ∆y = 0.2 µm.The longitudinal step size is taken to be ∆z = 2 µm to make calculation efficient.n 0 is taken to be the effective index of TE 1 mode in port 1.Although the fundamental and the higher order mode are guided in this device, since the relative index difference is very small (∆ ≈ 0.3%) and it is expected that the propagation angle may be small enough, the Fresnel approximation is employed also in this design example.A uniform medium (ρ = 0.3 in the design region) or a linear taper (ρ = 0.8 in core region, 0.3 in cladding region) is given as an initial structure.
Figure 11 shows that the output power of the TE 0 wave increases almost every iteration.The binarized optimized structure and the propagation wave are shown in Fig. 12.The normalized output power of the TE 0 wave is 0.911 in Fig. 12(a), and 0.800 in Fig. 12(b).An optimized structure depends on an initial structure in this design problem.The dependence on an initial structure is mainly determined by the property of the objective function.However, it is difficult to predict the property and a better initial structure.Nevertheless, in weakly guiding waveguides, the dependence is not usually significant if the composition of the objective function is simple, and single-mode and single-wavelength transmission is taken into account in the objective function.Actually, the dependence is not significant in the design of an S-bend and a 1:1 power splitter.Although the better structure is obtained in this design example when a uniform medium is given, it may be better to give the certain initial structure that is based on a known operation principle if the dependence on initial structure is strong.
This device can be used for a TE 0 -to-TE 1 converter by exchanging input and output ports, however, in this conversion, crosstalk may occur in port 1 because port 1 is not a single mode waveguide.In the superior one, when the TE 0 wave is launched into port 2, the normalized output power of the TE 1 wave and the crosstalk are 0.912 and −21.32 dB, respectively.Although the insertion loss is −0.4 dB, a low crosstalk TE 0 -to-TE 1 converter is realized.We investigate the influence of discretization on optimization results because there is possibility that the optimized structure depends on the step size.Although we tried designing with several different step sizes, almost same structure and performance are obtained in this problem, probably because of weakly guiding waveguide.
Finally, we compare the mode converter designed by using our approach with a conventional one.In previous studies on a mode order converter, an asymmetric branching waveguide has been used [22].Using the sizing optimization based on the steepest descent method, we design a TE 1 -to-TE 0 converter based on an asymmetric branching waveguide.Figure 13 is the design  The normalized output power of the TE 0 mode in the optimized structure is 0.846.Thus, it is shown that a mode converter that has higher performance can be obtained by using our design approach than that based on an asymmetric branching waveguide.

Conclusion
An efficient topology optimal design approach using the SVFD-BPM and the AVM has been proposed for the first time.The sensitivity analysis method reported in [16] has been extended for the case that the SVFD-BPM based on the ADIM is employed.S-bend waveguides, Y-branch type power splitters, and TE 1 -to-TE 0 mode converters are designed by using our design approach, then it is shown that the properties of the devices are improved every iteration in the optimization process.Therefore, the validity of our design approach is confirmed in the case of weakly guiding waveguides.
To apply our design approach to strongly guiding waveguides, it may be necessary to develop a constraint condition which prohibits appearance of the structure which is difficult to analyze by using the BPM.The applicability of our design approach to strongly guiding waveguides will be discussed in future works.

Funding
Japan Society for the Promotion of Science (JSPS) KAKENHI Grant Number 15K06009.

Fig. 4 .
Fig. 4. Optimized structure and propagation wave: φ is (a),(b) E x , (c),(d) H y , (e),(f) E y , and (g),(h) H x .The refractive index distribution in the design region is binarized in these structures.Insets (i)-(iii) are magnified views of the optimized structure.

8 Fig. 5 .
Fig.5.Normalized output power as a function of the iteration number in the design of an S-bend waveguide in the case that a uniform medium or a tilted waveguide is given as an initial structure.

Fig. 6 .
Fig. 6.Optimized structure and propagation wave H y : The initial structure is (a),(b) a uniform medium and (c),(d) a tilted waveguide.The refractive index distribution in the design region is binarized in these structures.Insets (i)-(iii) are magnified views of the optimized structure.

8 Fig. 8 .
Fig.8.Normalized output power as a function of the iteration number in the design of a 1:1 power splitter in the case that a uniform medium or a straight Y-branch waveguide is given as an initial structure.

Fig. 9 . 3 Fig. 10 .
Fig. 9. Optimized structure and propagation wave H y : The initial structure is (a),(b) a uniform medium and (c),(d) a straight Y-branch waveguide.The refractive index distribution in the design region is binarized in these structures.

8 Fig. 11 .Fig. 12 .
Fig. 11.Normalized output power of the TE 0 wave as a function of the iteration number in the mode converter.

Fig. 14 .
Fig. 14.Optimized structure of a TE 0 -to-TE 1 converter based on an asymmetric branching waveguide and propagation wave.

Table 1 .
The definitions of D xx and D yy .