Complex Photonic Integrated Resonators Modeled Using Scattering Matrix Methods

1. Introduction

With compact footprint and versatile configuration, integrated photonic resonators have become critical building blocks for photonic integrated circuits (PICs), with applications in a wide range of fields such as optical communications [1,2,3], photonic computing [4,5,6], nonlinear optics [7,8,9], sensing [10,11,12], and optical neural networks [13,14,15].

Modeling the spectral response of integrated photonic resonators is of fundamental importance for their use in different applications. The scattering matrix method (SMM) (also known as the transfer matrix method) [16,17], which is derived based on the Maxwell’s equations for electromagnetic waves, has been widely employed to model the spectral response of integrated photonic resonators with different device configurations [17,18,19,20,21,22].

Although in principle the SMM can be applied to model integrated photonic devices with arbitrary planar configurations [16,23], previous studies mainly focused on modeling devices with simple configurations. A key limitation comes from the fact that the traditional SMM relies on the manual derivation of scattering matrices,, and the obtained scattering matrices need to be multiplied in a specific sequence [24]. For devices with simple configurations, this allows for relatively straightforward calculation of spectral transfer functions. However, for complex integrated photonic resonators, particularly those with bidirectional light propagation, this process becomes much more complicated, which greatly limits its broader applicability.

In this paper, a universal approach based on the SMM is proposed to model integrated photonic resonators with complex structures. The modeling is achieved by dividing a device configuration into basic elements including directional couplers and connecting waveguides, followed by solving a system of linear scattering matrix equations using computational tools. Our approach can be applied to model devices with both unidirectional and bidirectional light propagation, and the simulated spectral response aligns well with measured results for practical devices. The modeling in our approach can also be simplified by dividing the device configuration into independent submodules such as microring resonators (MRRs) and Sagnac interferometers (SIs). Finally, we discuss the limitations of our approach induced by approximations in the modeling, together with strategies for improving modeling accuracy. Our approach provides an effective way for designing and optimizing complex integrated photonic devices across various integrated platforms.

2. Modeling of Devices with Unidirectional Light Propagation

In this section, we use the device configuration shown in Figure 1(a) as an example to show how to model a complex integrated photonic resonator based on the SMM in four steps. The device consists of a microring resonator (MRR) nested within another, with light propagating in only one direction in each waveguide.

First, we divide the device into several directional couplers and connecting waveguides, which are the basic elements in photonic integrated circuits that form integrated photonic resonators. The optical fields at the dividing points between these elements are denoted as Ei (i = 1–16).

Next, we establish a set of scattering matrix equations for the basic elements obtained in the first step, as shown in Table 1. These scattering matrix equations describe the relation between E_i at the input / output ports of the directional couplers and the connecting waveguides. For the directional coupler in Figure 1(b), the field transfer function can be expressed as [16]

$$\left[\begin{array}{c}{E}_{\mathrm{out}-1}\\ E_{\mathrm{out}-2}\end{array}\right]=\left[\begin{array}{cc}t& j\mathsf{\kappa }\\ j\mathsf{\kappa }& t\end{array}\right]\left[\begin{array}{c}{E}_{\mathrm{in}-1}\\ E_{\mathrm{in}-2}\end{array}\right],$$

(1)

where j = $\sqrt{-1}$, t and κ are the self-coupling and cross coupling coefficients, which satisfies the relation t ² + κ² = 1 when assuming lossless coupling, E_in-1, E_in-2, E_out-1, and E_out-2 are the input and output optical fields right before and after the coupling region. For the connecting waveguide in Figure 1(c), the field transfer function can be given by [23]

$$T=a{e}^{-j\mathsf{\phi }},$$

(2)

where a = e^−αL/2 is the round-trip transmission factor, with α and L denoting the power propagation loss factor and the waveguide length, respectively. In Equation (2), φ = 2πn_gL/λ is the round-trip phase shift, with n_g and λ denoting the group index and the wavelength, respectively. For each directional coupler, two equations can be derived, and one equation can be obtained for each connecting waveguide. This results in 2 × 4 + 1 × 6 = 14 equations in total for the device shown in Figure 1(a), which includes four directional couplers and six connecting waveguides.

Third, the system input is set. For instance, if we assume that there is only a continuous-wave (CW) input from Port 1, then another two equations can be obtained: E₁ = 1 and E₁₀ = 0. Here we set E₁ as 1 because the spectral transfer function at the output port, such as Port 2, is given by f_{Port 2} = E₃ / E₁. By setting E₁ to 1, f_{Port 2} = E₃ / E₁ = E₃, then the transfer function can be directly determined by calculating E₃ in the next step.

Finally, by solving all the linear equations obtained in the second and third steps, one can obtain the spectral transfer functions at the output ports. In Table 1, we summarize all the 16 equations for the device in Figure 1(a). In these equations, there are 16 variables E_i (i = 1–16), with the device’s structural parameters denoted using symbols and treated as constant coefficients. By solving this system of linear equations using computational tools (e.g., symbolic calculation in MATLAB), any of the variables E_i (i = 1–16) can be determined as a function of the device’s structural parameters. For example, the spectral transfer function at Port 2 can be given by f_{Port 2} = E₃ (T_i, t_i, κ_i). Assuming that the power propagation loss factor α and the group index n_g are constant for the device, the spectral transfer function at Port 2 can be expressed as a function of L_i, t_i, and κ_i, i.e., f_{Port 2} = E₃ (L_i, t_i, κ_i).

In Figure 2(a), we show the intensity and phase response spectra at Port 2 for different t_i based on the obtained spectral transfer function f_{Port 2} = E₂ (L_i, t_i, κ_i). Except for the varying parameters, other structural parameters were chosen as follows: L₁ = L₂ = L₃ = L₄ = 41.73 µm, L₅ = L₆ = 62.83 µm, and t₂ = t₄ = 0.98. By using the same method, we obtain the intensity and phase responses at Port 3 for different t_i, as shown in Figure 2(b). In our modeling, we assume that α = 101 m^-1 and group index n_g = 4.335 based on values obtained from our previously fabricated silicon photonic devices [20,25]. Unless otherwise specified, we use the same values for these two parameters in the modeling for all the devices in the following figures.

Our results in Figure 2 show good agreement with the results in Ref. [26], confirming the effectiveness of our method. Since the obtained spectral transfer functions (e.g., f_{Port 2} = E₂ (L_i, t_i, κ_i)) include all the structural parameters, the intensity and phase responses can be easily tailored by adjusting any of these parameters. This provides significant flexibility in designing and optimizing the spectral response of integrated photonic devices with different structural parameters. It should also be noted that our modeling, as mentioned above, is not limited to silicon photonic platforms. By using the corresponding values of α and n_g, it is applicable across various integrated platforms such as silicon nitride, doped silica, and lithium niobate [27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93].

3. Modeling of Devices with Bidirectional Light Propagation

For the device in Figure 1(a), light propagates in only one direction along each waveguide. In more complex integrated photonic resonators with light propagating bidirectionally in the waveguides, the modeling becomes more complicated. In this section, we use the device configuration shown in Figure 3(a) as an example to show how to model a complex integrated photonic resonator with bidirectional light propagation based on the SMM. The device consists of an add-drop (AD) MRR sandwiched between a pair of Sagnac interferometers (SIs), where the SIs introduce back reflection and hence bidirectional light propagation in such device.

To model the device in Figure 3(a), we first divide it into several directional couplers and connecting waveguides ‒ the same as what we did for the device in Figure 1(a). The optical fields at the dividing points between directional couplers and connecting waveguides are denoted as E_i (i = 1–16). Since there are optical fields traveling in two different directions at each dividing point, the optical fields traveling from left to right or in a clockwise direction are defined as “+”, and the opposite direction is defined as “−”.

In the second step, we establish 28 scattering matrix equations describing the relation between E_i⁺ and E_i^- at the input / output ports of the directional couplers and the connecting waveguides, as shown in Table 2. For each directional coupler, four equations can be derived, and two equations can be obtained for each connecting waveguide. Note that the number of equations is doubled due to the bidirectional light propagation in such device.

In the third step, we obtain another four equations by setting the system input. If we assume that there is only a CW input from Port 1, then the four equations are: E₁⁺ = 1, E₇⁻ = 0, E₁₀⁺ = 0, and E₁₆⁻ = 0. Here the number of equations is also doubled due to the bidirectional light propagation.

Finally, by solving all the 32 linear equations obtained previously (which include 32 variables, i.e., E_i⁺ and E_i^-, i = 1–16), we can obtain the spectral transfer functions at all the output ports. For example, the transfer function for forward transmission, with output at Port 2, is f_{Port 2} = E₁₆⁺ (L_i, t_i, κ_i), whereas the transfer function for backward reflection is f_{Port 1} = E₁^- (L_i, t_i, κ_i), with output at Port 1. It should be noted that as the number of equations increases or the device configuration becomes more complex, the use of computational software to solve the system of linear equations offers significant advantages as compared to manual derivation of the spectral transfer functions as those in Refs. [25,26].

Based on the obtained spectral transfer functions at Port 2 and Port 1, we plot the corresponding response spectra for different t_i, as shown in Figure 3(b) and Figure 3(c), respectively. Except for t_i, the other structural parameters are kept as constant as: L₁ = L₆ = 129.66 µm, L₂ = L₅ = 77.67 µm, L₃ = L₄ = 63.33 µm, and t₂ = t₃ = 0.95. Our results in Figure 3 also show good agreement with the results in Ref. [25], confirming the effectiveness of our method in modeling complex integrated photonic resonator with bidirectional light propagation.

4. Simplification by Dividing Device into Submodules

In modeling of the device configurations inFigure 1(a) and Figure 3(a), the devices were divided into several basic elements including directional couplers and connecting waveguides. To reduce the number of equations and simplify the calculation of the spectral transfer function, the device configuration can also be divided into several independent submodules such as MRRs and SIs, which are formed by the basic elements.

As illustrated in Figure 4(a), the AD-MRR at the center of the device in Figure 1(a) is treated as an independent module, with the field transfer functions at the through and drop ports denoted as T_RR and D_RR, respectively. This results in 12 equations in total obtained by using our method mentioned previously, as summarized in Table 3. In comparison to the 16 equations obtained for the same device configuration in Table 1, the number of equations is reduced, which helps to simplify the solving of the system of equations. By solving the 12 equations in Table 3, the spectral transfer function at Port 2 can be obtained, which is a function of T_MRR and D_MRR. We compared the response spectra obtained by solving the equations in Table 1 and Table 3 and found them to be identical. This verifies the validity of the new dividing method in Figure 4(a). It is also worth mentioning that when plotting the response spectra based on this dividing method, T_MRR and D_MRR should be further expressed as [94]

$$T_{\mathrm{RR}}={\displaystyle \frac{E_6}{E_5}}={\displaystyle \frac{E_{12}}{E_{11}}}={\displaystyle \frac{t_{\mathrm{RR}-1} - t_{\mathrm{RR}-2}{a}_{\mathrm{RR}}{e}^{-j{\mathsf{\phi }}_{\mathrm{RR}}}}{1 - t_{\mathrm{RR}-1}{t}_{\mathrm{RR}-2}{a}_{\mathrm{RR}}{e}^{-j{\mathsf{\phi }}_{\mathrm{RR}}}}} ,$$

(3)

$$D_{\mathrm{RR}}={\displaystyle \frac{E_{12}}{E_5}}={\displaystyle \frac{E_6}{E_{11}}}={\displaystyle \frac{ - {\mathsf{\kappa }}_{\mathrm{RR}-1}{\mathsf{\kappa }}_{\mathrm{RR}-2}\sqrt{a_{\mathrm{RR}}}{e}^{-j{\mathsf{\phi }}_{\mathrm{RR}}}}{1 - t_{\mathrm{RR}-1}{t}_{\mathrm{RR}-2}{a}_{\mathrm{RR}}{e}^{-j{\mathsf{\phi }}_{\mathrm{RR}}}}} ,$$

(4)

where j = $\sqrt{-1}$, a_RR and φ_RR are the round-trip transmission factor and phase shift along the ring in the AD-MRR, respectively. t_RR-i and κ_RR-i (i = 1, 2) are the self-coupling and cross coupling coefficients of the directional couplers in the AD-MRR, respectively.

Similarly, in Figure 4(b) we divided the device configuration in Figure 3(a) into an add-drop MRR, two SIs, and several connecting waveguides between them. By using this new dividing method, we obtained the same response spectra as those in Figure 3(b) and Figure 3(c). When plotting the response spectra, the forward and backward field transfer functions for the SIs were further expressed as

$$T_{\mathrm{SI}-i}={\displaystyle \frac{E_2^{+}}{E_1^{+}}}={\displaystyle \frac{E_6^{+}}{E_5^{+}}}=\left(t_{\mathrm{SI}-i}^2- {\kappa }_{\mathrm{SI}-i}^2\right)a_{\mathrm{SI}-i}{e}^{-j{\mathsf{\phi }}_{\mathrm{SI}-i}}, (i=1,2)$$

(5)

$$R_{\mathrm{SI}-i}={\displaystyle \frac{E_1^{-}}{E_1^{+}}}={\displaystyle \frac{E_5^{-}}{E_5^{+}}}=2j{t}_{\mathrm{SI}-i}{\kappa }_{\mathrm{SI}-i}{a}_{\mathrm{SI}-i}{e}^{-j{\mathsf{\phi }}_{\mathrm{SI}-i}}, (i=1,2)$$

(6)

where j = $\sqrt{-1}$, a_SI-i and φ_SI-i (i = 1, 2) are the round-trip transmission factor and phase shift along rings in the SIs, respectively. t_SI-i and κ_SI-i (i = 1, 2) are the self-coupling and cross coupling coefficients of the directional couplers in the SIs, respectively.

As shown in Table 4, only 10 equations are established by using the new dividing method, in contrast to 32 equations obtained for the dividing method in Figure 3(a). This indicates that although dividing device into basic elements is a universal approach, dividing it into submodules simplifies the process of solving equations. This is particularly true for devices composed of multiple submodules, where the advantages of the new division method become even more evident.

In Figure 5, we model a more complex device by dividing it into eight SIs and establishing 32 linear equations, as summarized in Table 5. The forward and backward field transfer functions for the SIs are denoted as T_SI-i and R_SI-i (i = 1 ‒ 8), respectively. Dividing the device into basic elements would result in a system of 64 equations, compared to the 32 equations shown in Table 5. By solving the system of 32 equations to obtain E₁₆⁺ and E₁⁻, we obtain the spectral transfer functions at Port 2 and Port 1 (reflection), respectively. The corresponding response spectra are plotted in Figure 5(b) and Figure 5(c), which are consistent with the results in Ref. [95] and further confirms the effectiveness of our method in modeling devices with complex configurations.

When dividing complex devices into submodules, it’s important to ensure that the submodules are independent. Here, "independent" means that the submodules exchange energy with other parts exclusively through connecting waveguides. For the device shown in Figure 6a formed by a self-coupled waveguide, even though there are multiple SIs within the device, it cannot be divided into submodules of SIs. This is because the adjacent SIs are mutually coupled with energy exchange through directional couplers. In this case, the spectral transfer functions can be obtained by using the universal dividing method, which divides the device into basic elements of directional couplers and connecting waveguides.

Table 6 shows a system of 56 equations obtained by using the universal dividing method for the device in Figure 6a. By solving these equations to obtain E₂₈⁺ and E₁⁻, we obtain the spectral transfer functions at Port 1 and Port 2, respectively. The corresponding response spectra for different t_i are shown in Figure 6(b) and Figure 6(c). Except for t_i, all other structural parameters are kept constant as: L₁ = 180 μm, L₂ = 60 μm and L₃ = 120 μm. The response spectra also show a good agreement with the results in Ref. [24], providing additional evidence of the effectiveness of our method.

5. Deviations Induced by Approximations

As evidenced by the results in previous sections, our method proves effective for modeling integrated photonic resonators and shows advantages for modeling those with complex configurations. Despite this, the method still has limitations that could cause deviations between the simulation results and measurements for practical devices. In this section, we discuss the limitations of our method to model integrated photonic resonators and methods for improving modeling accuracy.

In our modeling, we assume the self-coupling coefficient t and cross-coupling coefficient κ of the directional couplers to be wavelength-independent constant, while this holds true only within a certain wavelength range. According to Ref. 34 95 32], the field cross-coupling coefficient κ of a directional coupler with a coupling length of L_c can be given by

$$\mathsf{\kappa }=\sin \left({\displaystyle \frac{\mathsf{\pi }}{2}}·{\displaystyle \frac{L_c}{L_x}}\right),$$

(7)

where L_x is the cross-over length, defined as the minimum distance at which optical power completely transfers from one waveguide to the other. The L_x can be further expressed as [95]

$$L_x={\displaystyle \frac{\lambda }{2(n_{\mathrm{eff}, \mathrm{even}}-n_{\mathrm{eff}, \mathrm{odd}})}}.$$

(8)

where λ is the light wavelength, n_{eff, even} and n_{eff, odd} are the effective indices of the two fundamental eigenmodes of the coupled waveguides, respectively. From Eqs. (7) and (8), it can be seen that the coupling strength of the directional coupler is actually wavelength dependent due to the existence of dispersion (which leads to wavelength-dependent values of n_{eff, even} and n_{eff, odd}). Therefore, the variation in the coupling strength of the directional couplers with wavelength can no longer be ignored when modeling the device’s response spectra over a broad bandwidth (typically > 30 nm).

Another potential limitation arises from the fact that we use a constant waveguide group index n_g instead of wavelength-dependent waveguide effective index n_eff when calculating the phase shift along waveguides. In fact, the relation between n_g and n_eff can be expressed as [94]:

$$n_g=n_{\mathrm{eff}}\left(\mathsf{\lambda }\right)- \mathsf{\lambda }{\displaystyle \frac{d{n}_{\mathrm{eff}}\left(\mathsf{\lambda }\right)}{d\mathsf{\lambda }}},$$

(9)

where n_eff (λ) is the effective index as a function of light wavelength λ. The group index n_g is widely used for calculating the free spectral ranges (FSRs) of integrated photonic resonators. For example, the FSR for a MRR can be approximately given by [94]

$${\mathsf{\lambda }}_{F\mathrm{SR}} \text{≈} {\displaystyle \frac{{\mathsf{\lambda }}_0^2}{n_{g}L}}$$

(10)

where λ₀ is the resonance wavelength and L is the ring circumference.

Similar to the cross-coupling coefficient κ of a directional coupler, the group index n_g can be regarded as a constant coefficient only within a specific wavelength range. In contrast to introducing n_eff at different wavelengths, using a constant n_g can greatly simplify the plotting of the response spectra based on the spectral transfer function. However, this approximation only holds for modeling spectral response over a relatively narrow wavelength range (e.g., < 30 nm). When modeling the spectral response over a broad bandwidth or for nonlinear optical devices that are sensitive to dispersion-induced phase mismatch [8,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128], the wavelength-dependent n_eff should be used to achieve more accurate simulation results.

6. Conclusion

In summary, a universal approach for modeling integrated photonic resonators with complex structures is proposed based on the SMM. By dividing the device configuration into basic elements such as directional couplers and connecting waveguides, our approach shows effectiveness in modeling integrated photonic resonators with both unidirectional and bidirectional light propagation, with the simulated spectral response agreeing well with experimental results. By dividing the device configuration into independent submodules MRRs and SIs, the modeling in our approach can also be simplified. Finally, we discuss the limitations arising from approximations in our modeling and the corresponding strategies for improving modeling accuracy. Our approach offers an efficient way for designing and optimizing complex integrated photonic devices, which is applicable across a wide range of integrated platforms.