Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Extracting causal impulse responses from baseband band-limited S-parameters of passive photonic integrated circuits for time-domain simulations

Open Access Open Access

Abstract

The dispersive characteristics and wavelength-dependent behaviors of passive photonic integrated circuits (PICs) can be well described by S-parameters. However, circuit-level simulations of PICs that commonly consist of both passive and active components have to be conducted in the time domain. Thus, S-parameters need to be converted into a time-domain representation without losing accuracy and violating physical properties (e.g., causality). To address this issue, this paper proposes an approach for extracting causal impulse responses of passive PICs by extrapolating their baseband band-limited S-parameters. The method is efficient and robust for a wide range of passive PICs.

© 2024 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement

1. Introduction

Due to their high refractive index contrast and compatibility with CMOS fabrication technology, silicon photonic integrated circuits (PICs) are experiencing rapid scaling, progressing from integrating thousands to millions of devices [1]. This significant growth is primarily driven by the escalating demand for optical communication, sensing, computing, and various other applications. Similar to electronic integrated circuits, the design of large-scale PICs heavily relies on precise time-domain models and circuit-level simulations [24]. Moreover, PICs are often integrated with electronic integrated circuits to leverage the advantages of both technologies across various applications. This integration calls for the co-design and co-optimization of electronic-photonic integrated circuits, a process that also demands a substantial amount of accurate time-domain evaluations.

Passive PICs are normally characterized by their frequency (or wavelength) response within a limited bandwidth, such as ${[1.5,1.6]}$ ${\mu }$m. It is required to convert this band-limited information into a time-domain model for circuit-level evaluation of PICs. Note that band-limited S-parameters of passive PICs that are simulated or measured in the optical frequency range need to be shifted to baseband using complex envelope representation. This process separates the ultra-high frequency information from the system, enabling efficient time-domain simulations [46]. The transformed set of parameters is referred to as baseband band-limited S-parameters (BBLS). Contrary to conventional S-parameters, BBLS are not Hermitian symmetric.

Given the exceptional nature of BBLS, the work [6] has developed the complex vector fitting (CVF) method, a customized variant of the vector fitting (VF) method [7]. The approach leverages on a reformulation of Sanathanan–Koerner iteration to approximate BBLS using a pole-residue transfer function. Once the optimal poles and residues are identified, the transfer function can be analytically transformed to a continuous-time state-space model for circuit simulations. While this method proves to be accurate and robust for the majority of passive PICs, it may suffer from computational expenses and lead to bulky models when dealing with S-parameters featuring complicated resonances and delays.

Other methods are to implement the frequency dependent behavior of passive PICs by designing digital finite impulse response (FIR) or infinite impulse response (IIR) filters so that their frequency responses approximate that of the passive PICs of interest [8]. These models are usually solved in dynamic data flow simulators that are available in commercial PICs simulation software [9], such as Lumerical Interconnect. However, building accurate FIR/IIR models for the baseband scattering parameters of passive PICs is challenging. In [8], it is mentioned that the accuracy provided by FIR models drops near the edges of modeling frequency range even for the best designed FIR models, which could potentially be attributed to the effects of windowing. The IIR models are also adopted in mainstream commercial simulators, such as Lumerical Interconnect. But, the methods used to build the FIR/IIR models are unknown and the accuracy of simulation results is not always guaranteed.

An alternative method is to convert BBLS into discrete impulse responses by the inverse discrete Fourier transform (IDFT), and then circuit simulations are conducted through convolutions. Applying the IDFT directly on BBLS is not feasible since it could lead to severe causality violations due to limited bandwidth. Since a similar issue also exists for conventional band-limited S-parameters, various methods [1014] have been proposed to extrapolate the conventional band-limited S-parameters to enforce the Kramers-Kronig (K-K) relation, thereby achieving causal impulse responses. However, these methods cannot be directly applied to BBLS since BBLS are not Hermitian symmetric, unlike conventional S-parameters, which exhibit even symmetric magnitude and odd symmetric phase. Another significant limitation of these methods [1014] is that the extrapolated spectrum fails to reach zero at the edge of the extended frequency range, restricting the flexibility of the time steps for the obtained impulse response. When there is a need to adjust the time step to accommodate different requirements or applications, the entire extrapolation process has to be repeated.

Therefore, we propose a novel extrapolation scheme to ensure a causal impulse response while guaranteeing the extrapolated spectra rapidly roll off to zero. This approach directly enforces causality on impulse responses which are represented with linear combinations of S-parameters through the IDFT. The decay-to-zero property is implemented by extrapolating BBLS with a specially constructed polynomial. All these constraints are formulated into an efficiently solvable least squares problem. The resulting causal spectrum can then serve as a valuable model for passive PICs in diverse circuit-level simulations. When an adjustment to the time step in impulse responses is necessitated to satisfy different applications, an impulse response with the desired time step can be obtained promptly by zero-padding the extrapolated spectrum and then performing the IDFT. Time-domain simulations can be efficiently performed by convolving the impulse response with given input signals in any dynamic data flow simulators that are widely adopted by mainstream PIC simulators (e.g., Lumerical Interconnect, VPIphotonics) [8,9]. In contrast, the CVF models are inherently continuous systems of first-order ordinary differential equations and are more suitable to be solved using differential algebraic equation solvers [15], such as SPICE-like simulators. Consequently, time-domain simulations with the proposed model are much more efficient than those with the CVF models especially when the CVF models have many poles, as demonstrated by the numerical examples in the paper.

The rest of the paper is organized as follows. Section 2 introduces baseband signal and BBLS of PICs. Section 3 provides a detailed explanation of our proposed method. In Section 4, we present two examples of passive PICs to showcase the effectiveness of our approach. Finally, Section 5 offers a summary of the paper.

2. Baseband representations for passive photonic integrated circuits

Lightwaves in PICs are normally modeled with their electric field. However, this electric field can oscillate at frequencies of hundreds of terahertz. For optical interconnect applications, the frequencies could fall within a range, for example [187.5, 200] THz, corresponding to wavelengths of [1.5, 1.6] $\mu$m. In reality, the fast oscillating electric field is hard to measure, and it is not directly utilized in practical applications. Instead, it is the modulating radio frequency (RF) signals on optical carriers that are the valuable information of interest. Therefore, the modulated optical signals

$$E(t)=A(t) \cos (2\pi f_c t+\phi(t))$$
are represented with the complex baseband equivalent by removing the carrier
$$E_b(t)=A(t) e^{j\phi(t)},$$
where $A(t)$ and $\phi (t)$ are amplitude and phase modulating RF signals, respectively, and $f_c$ is the optical carrier frequency. Note that, the phase signal $\phi (t)$ wraps around $2\pi$, resulting in multiple solutions and discontinuities. Thus, baseband signals are often represented with their real and imaginary parts rather than amplitude and phase in PICs simulations. In the scenario of quadrature amplitude modulation, the real and imaginary parts correspond to the in-phase and quadrature components, respectively.

Accordingly, to accommodate the baseband signal, the baseband equivalent of passive PICs can also be derived. Assuming that the S-parameters of a passive PIC under study is simulated or measured in the frequency range of interest, its baseband representation is obtained by shifting the S-parameters downwards by carrier frequency. The red solid lines in Fig. 1 depict the magnitude and phase of BBLS while the blue solid lines represent that of the measured or simulated S-parameters in optical frequency range. It is evident that BBLS do not comply with the Hermitian symmetric constraint, distinguishing them from conventional S-parameters. It is important to note that baseband S-parameters are a complex envelope representation that are widely used in communication system to simplify the modulation and demodulation process [16].

 figure: Fig. 1.

Fig. 1. Magnitude and phase of the baseband band-limited S-parameters (red solid lines) and the extrapolated spectrum (red dash lines) of passive PICs.

Download Full Size | PDF

To extract causal impulse responses from BBLS, it is necessary to extrapolate both the lower (negative) and higher (positive) sides of the spectrum. The extrapolation should be performed to satisfy the Kramers-Kronig (K-K) relation while not enforcing the Hermitian symmetric constraint. This requirement immediately renders ineffective the extrapolation approaches [1014] that have been proposed for electronic circuits.

The CVF technique [6] has been proposed to handle the unique characteristics of BBLS and it constructs continuous-time state-space models in the following form for passive PICs

$$\left\{ \begin{aligned} \frac{\text{d}{{\boldsymbol x}}(t)}{\text{d}t} & = {\boldsymbol A}{\boldsymbol x}(t)+{\boldsymbol B}{\boldsymbol a}(t)\\ {\boldsymbol b}(t) & = {\boldsymbol C}{\boldsymbol x}(t)+{\boldsymbol D}{\boldsymbol a}(t), \end{aligned} \right.$$
where ${\boldsymbol a}(t)$ and ${\boldsymbol b}(t)$ are the baseband incident and reflected waves, respectively, and contain magnitude and phase information, thus they are complex-valued. ${\boldsymbol x}(t)$ represents a state vector, and ${\boldsymbol A}$, ${\boldsymbol B}$, ${\boldsymbol C}$, ${\boldsymbol D}$ are constant and complex state-space parameters. The CVF method is inherently good at modeling resonances in S-parameters. Nevertheless, when modeling frequency responses that include intricate resonances and time delays, it could suffer from high computational cost and lead to cumbersome models [1214].

Our objective in this study is to devise a novel approach to extrapolate BBLS at both lower and higher frequencies, satisfying the K-K relation, while ensuring that the extrapolated spectra rapidly decrease to zero at the edges, as indicated by the red dash lines in Fig. 1. It is important to note that the requirement for the extrapolated spectra to decay to zero at the edges also enforces periodicity of the overall spectrum. The proposed method will be compared to the CVF technique in Section 4.

3. Extrapolation algorithm

To enforce causality, band-limited S-parameters must be extrapolated so that the real and imaginary parts of the full spectrum comply with the K-K relation. It is often implemented in the way that the real (or imaginary) part is first extrapolated up to a certain frequency and then the imaginary (or real) part is calculated via the Hilbert transform (HT)

$$\begin{aligned} S_{\Re}(f) = \frac 1\pi p.v. \int_{-\infty}^{\infty} {\frac{S_{\Im}(f')}{f-f'}} df', \\ S_{\Im}(f) ={-}\frac 1\pi p.v. \int_{-\infty}^{\infty} {\frac{S_{\Re}(f')}{f-f'}} df', \end{aligned}$$
where $p.v.$ is the principal value of the integral, and $S_{\Re }(f)$ and $S_{\Im }(f)$ are the real and imaginary parts of $S(f)$. The real part can be extrapolated with Gaussian functions [10] or polynomials [1113]. Then parameters in the extrapolation function are optimized iteratively to minimize the error between $S_{\Im }(f)$ and the imaginary part computed from HT of the extrapolated real part [10,11]. The work [12,13] takes a further step by integrating the real part extrapolation and imaginary part error minimization processes into a least square scheme, thereby solving the problem efficiently. However, this approach cannot take full control of the extrapolated spectrum in the extended frequency range to force it to drop to zeros since the extrapolation poses constraints on either the real or imaginary part, rather than on both simultaneously. In [14], discrete HT is utilized to extrapolate the real and imaginary parts at the same time without any extrapolation functions. But the work incorporates computationally expensive optimization methods to find an accurate solution and it does not guarantee a declining extrapolation spectrum either. In the following, we introduce a new extrapolation approach that is particularly developed for BBLS. The method enforces causality directly on impulse responses via IDFT, without relying on the K-K relation.

Assume that the BBLS {$f_i$, $S(f_i)$} under study is obtained in the baseband frequency range $f_i\in$[$f_{min}$, $f_{max}$], where $f_{min}<0<f_{max}$. The lower and higher bands of the BBLS will be extrapolated to $f_{L}$ and $f_{R}$, respectively. And the corresponding extrapolated spectrum in [$f_{L}$, $f_{min}$) and ($f_{max}$, $f_{R}$] are denoted as $S_{L}$ and $S_{R}$. $S_{full}$ represents the entire spectrum composed of $S_{L}$, $S$ and $S_{R}$.

To ensure that $S_{L}$ smoothly decays to zeros at $f_{L}$, a high-order polynomial is adopted to generate $S_{L}$

$$S_{L}(f_i)=\sum_{m=4}^{M}a_m (f_i-f_L)^m,{f_L\leq f_i < f_{min}}.$$

Similarly, $S_{R}$ can be expressed as

$$S_{R}(f_i)=\sum_{k=4}^{K}b_k (f_i-f_{R})^k,{f_{max}<f_i\leq f_{R}} ,$$
where $a_m$ and $b_k$ are the polynomial coefficients and will be determined later. Note that, these extrapolation functions are constructed so that $S_{L}$ and $S_{R}$ equal zero at $f_{L}$ and $f_{R}$, respectively. The locations of $f_{L}$, $f_{min}$, $f_{max}$ and $f_{R}$, as well as the modeling strategy, can be better observed in Fig. 1. Moreover, the choice of starting the polynomial order from four is intended to ensure that its higher-order derivatives are also zero. This constraint enforces a decay of $S_{L}$ and $S_{R}$ towards zero with a zero angle to the frequency axis. As a result, padding zeros to $S_{L}$ and $S_{R}$ creates minimal discontinuity to the spectrum, thus only limited causality violation will be introduced. This approach allows for generating impulse responses with smaller time steps (or higher sampling frequencies), catering to various requirements or applications, without the need for repeating the extrapolation process.

The discrete impulse response of $S_{full}$ can be calculated with the IDFT

$$h(n)=\sum_{p={-}P/2}^{P/2-1} S_{full} (p) \varphi_{np} ,$$
where $\varphi _{np}=\text {exp}({-j2\pi {np}/{P}})$. For the sake of simplifying the explanation of the proposed method, an integer (e.g. $p$) is adopted to index $S_{full}$, consistently used throughout the rest of the paper, in place of the notation $f_i$. $P$ is the total sample number of $S_{full}$ and is assumed to be even.

Mathematically, the K-K relation is equivalent to that impulse responses are vanishing for $t<0$ [17]. For the $h(n)$ in Eq. (7), causality indicates that

$$\sum_{p={-}P/2}^{P/2-1} S_{full} (p) \varphi_{np}=0, \quad -P/2 \leq n <0.$$

Since $S_{full}$ is composed of $S_L$, $S$, and $S_R$, there is

$$\begin{aligned} \sum_{p={-}P/2}^{{-}P/2+L_N} & S_{L} (p) \varphi_{np}+\sum_{p={-}P/2+L_N+1}^{{-}P/2+L_N+N} S(p) \varphi_{np}+ \sum_{p={-}P/2+L_N+N+1}^{P/2-1} S_{R} (p) \varphi_{np}=0, \, -P/2 \leq n<0 \end{aligned}$$
where $L_N+1$ and $N$ are sample numbers of $S_L$ and $S$, respectively. Sample number of $S_R$ is denoted as $R_N$, satisfying $L_N+1+N+R_N=P$. Representing Eq. (9) in matrix form and moving the known term $S$ to the right hand side lead to
$${\boldsymbol A}_{L} {\boldsymbol S}_{L}+{\boldsymbol A}_{R} {\boldsymbol S}_{R}={-}{\boldsymbol AS},$$
where
$$\begin{aligned} {\boldsymbol A}_{L}= & \begin{bmatrix} \vec{A}_{{-}P/2} & \cdots & \vec{A}_{{-}P/2+L_N} \end{bmatrix}\in \mathbb{C}^{(P/2)\times(L_N+1)} \\ {\boldsymbol A}= & \begin{bmatrix} \vec{A}_{{-}P/2+L_N+1} & \cdots & \vec{A}_{{-}P/2+L_N+N} \end{bmatrix}\in \mathbb{C}^{(P/2)\times N} \\ {\boldsymbol A_{R}}= & \begin{bmatrix} \vec{A}_{P/2-R_N} & \cdots & \vec{A}_{P/2-1} \end{bmatrix}\in \mathbb{C}^{(P/2)\times R_N} \\ \vec{A}_i= & \begin{bmatrix} \varphi_{{-}P/2, i} & \cdots & \varphi_{{-}1, i} \end{bmatrix}^T\in \mathbb{C}^{(P/2)\times 1}\\ \end{aligned}$$

The matrices ${\boldsymbol A}_{L}$, ${\boldsymbol A}_{R}$, and ${\boldsymbol A}$ can be immediately computed via $\varphi _{np}=\text {exp}({-j2\pi {np}/{P}})$.

The extrapolation vectors $S_L$ in Eq. (5) and $S_R$ in Eq. (6) can also be formulated into matrix form as follows

$$\begin{aligned} {\boldsymbol S_{L}}={\boldsymbol F} \cdot {\boldsymbol a}, \\ {\boldsymbol S_{R}}={\boldsymbol G} \cdot {\boldsymbol b}, \end{aligned}$$
where
$$\begin{aligned} & {\boldsymbol F}= \begin{bmatrix} (f_{{-}P/2}-f_{L})^{4} & \cdots & (f_{{-}P/2}-f_{L})^{M}\\ \vdots & \vdots & \vdots \\ (f_{-\frac{P}{2}+L_N}-f_{L})^{4} & \cdots & (f_{-\frac{P}{2}+L_N}-f_{L})^{M} \end{bmatrix} \in \mathbb{R}^{(L_N+1)\times (M-3)}\\ & {\boldsymbol G}= \begin{bmatrix} (f_{P/2-R_N}-f_{R})^{4} & \cdots & (f_{P/2-R_N}-f_{R})^{K}\\ \vdots & \vdots & \vdots\\ (f_{P/2-1}-f_{R})^{4} & \cdots & (f_{P/2-1}-f_{R})^{K} \end{bmatrix} \in \mathbb{R}^{R_N\times (K-3)}\\ & {\boldsymbol a}= \begin{bmatrix} a_{4} & \cdots & a_{M} \end{bmatrix}^T \in \mathbb{C}^{(M-3)\times 1}\\ & {\boldsymbol b}= \begin{bmatrix} b_{4} & \cdots & b_{K} \end{bmatrix}^T \in \mathbb{C}^{(K-3)\times 1}\\ \end{aligned}$$

It is important to note that $f_{-P/2}=f_{L}$ and $f_{P/2-1}=f_{R}$, which ensures extrapolation polynomials decay to zero at these frequency points.

Combining Eq. (10) with Eq. (12) can derive the least square problem

$$\begin{bmatrix} {\boldsymbol A}_{L} {\boldsymbol F} & {\boldsymbol A}_{R} {\boldsymbol G} \end{bmatrix} \begin{bmatrix} \boldsymbol{a}\\ \boldsymbol{b} \end{bmatrix} ={-}{\boldsymbol A S}.$$

To preserve continuity at the joints ($f_{min}$ and $f_{max}$) of the BBLS and the extrapolated spectra, several extra constraints are added to the least square problem

$$\begin{bmatrix} {\boldsymbol A}_{L} {\boldsymbol F} & {\boldsymbol A}_{R} {\boldsymbol G} \\ \boldsymbol{F_c} & \boldsymbol{0}\\ \boldsymbol{0} & \boldsymbol{G_c} \end{bmatrix} \begin{bmatrix} \boldsymbol{a}\\ \boldsymbol{b} \end{bmatrix} = \begin{bmatrix} -{\boldsymbol AS}\\ \boldsymbol{S_c} \end{bmatrix},$$
where
$$\begin{aligned} {\boldsymbol F}_c= & \begin{bmatrix} (f_{-\frac{P}{2}+L_N+1}-f_{L})^{4} & \cdots & (f_{-\frac{P}{2}+L_N+1}-f_{L})^{M}\\ \vdots & \vdots & \vdots\\ (f_{-\frac{P}{2}+L_N+4}-f_{L})^{4} & \cdots & (f_{-\frac{P}{2}+L_N+4}-f_{L})^{M} \end{bmatrix} \in \mathbb{R}^{4\times (M-3)}\\ {\boldsymbol G}_c= & \begin{bmatrix} (f_{P/2-R_N-4}-f_{R})^{4} & \cdots & (f_{P/2-R_N-4}-f_{R})^{K}\\ \vdots & \vdots & \vdots\\ (f_{P/2-R_N-1}-f_{R})^{4} & \cdots & (f_{P/2-R_N-1}-f_{R})^{K} \end{bmatrix} \in \mathbb{R}^{4\times (M-3)}\\ \boldsymbol{S_c}= & \begin{bmatrix} S({-{P}/{2}+L_N+1}) \cdots S({-{P}/{2}+L_N+4})\quad S({P/2-R_N-4}) \cdots S({P/2-R_N-1}) \end{bmatrix}^{T} \end{aligned}$$

This requirement enforces that four points, located at the joints between $S_L$ and $S$, as well as between $S$ and $S_R$, must be equal to each other. The use of four points is recommended because fewer points would not be sufficient to ensure smoothness at the joints, whereas using more points might introduce additional errors when solving the least square problem.

By solving the least squares problem in Eq. (15), we can determine the polynomial coefficients in ${\boldsymbol a}$ and ${\boldsymbol b}.$ These coefficients are then used to calculate the extrapolation spectra $S_L$ and $S_R$ via Eq. (12). Note that the maximum polynomial order $M$ and $K$, as well as the number of extrapolated points $L_N$ and $R_N$, may have an impact on the accuracy and passivity of the recovered spectrum. To select these integer parameters effectively, a similar method proposed in [12,13] can be adopted.

The impulse response $h(n)$ can be immediately computed with $S_{full}=[S_L,S,S_R]$ via the inverse fast Fourier transform (IFFT) algorithm. It is worth mentioning that solving the overdetermined system in Eq. (15) with the least square methods will inevitably introduce errors, resulting in a non-zero impulse response at $n<0$. In this work, the non-zero impulse response at $n<0$ is extremely small (on the order of 1e-5) due to the well-conditioned nature of Eq. (15). Subsequently, the non-zero portion before $n=0$ can be set to zero, and the signal is converted back to the frequency domain by fast Fourier transform (FFT), leading to a causal spectrum $\bar {S}_{full}$. It is demonstrated in Section 4 that the error induced by this operation has minimal effects on the spectrum. Note that, this issue exists in all methods that include solving a least square problem. For instance, in [12,13], the imaginary parts calculated from the HT of the real parts of S-parameters are used to replace the original imaginary parts.

4. Numerical examples

4.1 Directional coupler

In this example, the proposed extrapolation method is applied to a silicon-on-insulator based directional coupler (DC). The coupling length is 10 $\mu$m, the gap between coupling region is 0.2 $\mu$m, and the waveguide has a thickness of 220 nm and a width of 500 nm. The structure is simple and not shown here. The transmission $S_{13}$ are obtained via electromagnetic simulations in Lumerical FDTD at 101 frequency points in the range [187.5, 200] THz, corresponding to the wavelength [1.5, 1.6] $\mu$m. By downshifting the S-parameters by $f_c$=193.75 THz, we get the BBLS $S$, which are plotted in Fig. 2(a) with a blue solid line. Please note that the magnitude of the BBLS exhibits a gradual variation, transitioning from 0.7 to 0.9, resembling a nearly linear slope. The time-domain output can be approximated by simply multiplying the input with $S_{13}(f_c)$ only when the bandwidth of the input signal for $S_{13}$ is extremely narrow. However, in cases where the input signal has a broader bandwidth, the frequency-dependent effects must be considered.

 figure: Fig. 2.

Fig. 2. (a) The full spectrum $S_{full}$ of the DC after extrapolation. (b) The impulse responses computed via $S_{full}$ and via directly padding zeros to the BBLS.

Download Full Size | PDF

By applying the proposed technique, the extrapolated spectrum $S_L$ and $S_R$ are computed with the highest polynomial order $M=K=10$ and an extended frequency $-f_L=f_R=12.5$ THz. The full spectrum is very smooth at the joints and the extrapolated portions approach zero with a zero angle to the axis, as demonstrated in Fig. 2(a). Despite the magnitude of the spectrum $S$ being close to one (at $f_{min}$) and prone to violating passivity, no such issue occurs in the extrapolated spectrum.

Then, the impulse response is calculated from $S_{full}$. For comparison, we also extrapolate the BBLS by directly padding zeros to extend it to the same frequency range and calculate its impulse response, which demonstrates considerable causality violations. Two impulse responses are illustrated in Fig. 2(b). Please note that the impulse responses of baseband systems are complex numbers, thus the amplitude (or absolute values) is shown in the plot. The proposed impulse response before $n=0$ is not exactly zero but rather extremely small, less than 1e-5 in this case, owing to the least square error. Therefore, this part can be directly zeroed to compute a completely causal spectrum $\bar {S}_{full}$ without really degrading the accuracy. Indeed, Fig. 3(a) and (b) presents the error of the magnitude and phase between $\bar {S}_{full}$ and $S$. We can observe that the maximum absolute error is below -80 dB. Considering that the existing extrapolation methods are not applicable to BBLS, we compare the proposed technique to the CVF modeling method. Thus, a CVF model is built for the BBLS with 8 poles after 10 iterations, achieving a maximum absolute error less than -60 dB, as depicted side by side with the proposed method in Fig. 3(c) and (d). The results clearly demonstrate the superior accuracy of the proposed method over the CVF method.

 figure: Fig. 3.

Fig. 3. Modeling accuracy for the magnitude and phase of the BBLS under analysis. (a)-(b) Magnitude and phase of the full causal spectrum $\bar {S}_{full}$ compared to that of the BBLS. (c)-(d) Magnitude and phase of the CVF model compared to that of the BBLS.

Download Full Size | PDF

Time-domain simulations are carried out for both models using an input signal generated from the sinc function. The input signal shares the same $\Delta t=1/(2f_R)=0.04$ ps as the computed impulse response and exhibits a rectangular-shaped frequency spectrum with a bandwidth of 2.5 THz. The causal impulse response associated with $\bar {S}_{full}$ is convolved with the input signal in SIMULINK, resulting in the generation of the output signal. In order to establish a benchmark, the CVF model is also simulated by solving ordinary differential equations within the same SIMULINK environment. It’s important to highlight that SIMULINK is chosen for its similarity to widely-used commercial simulators, such as Lumerical Interconnect and VPIphotonics, due to its dynamic data flow simulation nature. Note that, at the time of this research, these commercial simulators have not provided interfaces that we can import the proposed models and CVF models into them to perform time-domain simulations. There is also a lack of standard modeling languages for PICs that is supported across different commercial PIC simulators.

The input and resulting outputs from the two simulations are displayed in Fig. 4. Remarkably, it is observed that the outputs from two models exhibit a high level of agreement, with an root mean square (RMS) error of 7.7e-4. As far as the efficiency is concerned, computing the proposed model is 3.2 times faster than building the CVF model in this case. The time-domain simulations for both methods consume comparable amounts of time because the input signal is very short in duration, and the CVF model has a limited number of poles.

 figure: Fig. 4.

Fig. 4. Time-domain input and outputs for the proposed method and the CVF model.

Download Full Size | PDF

When the same DC model is reused in other applications where the input signals have much higher bandwidth and smaller time steps, $\Delta t$ for the impulse response should also become smaller to align with that of input signals. In this example, we assume the bandwidth of the input signal increases to 5 THz and time steps decrease to 0.02 ps to ensure the smoothness of the sinc curve. Now, the time steps for the impulse response should also be reduced to 0.02 ps for conducting convolutions. Thanks to that the extrapolated spectrum falls to zero at the edges of the extended frequency range, we can directly padding zeros at the beginning and the end of ${S}_{full}$, thereby increasing the sampling frequency. Then, a new impulse response with smaller $\Delta t$ can be computed via IFFT, as shown in Fig. 5(a). Note that, the new calculated impulse response also exhibit weak causality violations. We can enforce causality by zeroing out the non-zero impulse response values for $n < 0$, which introduces only minimal errors in the frequency domain. This is demonstrated in Fig. 5(b) where the corresponding frequency response is calculated from the new causal impulse response and compared to the BBLS. The maximum absolute error is still below -68 dB. The new causal impulse response is then utilized for convolution with the input with smaller $\Delta t$. The time-domain results are once again compared to the CVF simulation results in Fig. 6, demonstrating an RMS error of 0.0014 and reaffirming an exceptional agreement between the two.

 figure: Fig. 5.

Fig. 5. (a) The impulse response calculated by padding zeros to $S_{full}$. (b) The frequency spectrum of the new impulse response, compared to the BBLS.

Download Full Size | PDF

 figure: Fig. 6.

Fig. 6. The time-domain simulation results of the CVF model and the impulse response with $\Delta t=0.02$ ps.

Download Full Size | PDF

4.2 MZI-based demultiplexer

A Mach–Zehnder interferometer (MZI) based eight-way demultiplexer, as shown in Fig. 7, is designed and simulated in IPKISS. The power couplings of all DCs at 1.55 $\mu$m are indicated in the figure. The transmission of one channel at 251 frequency points within the range of [193.22, 193.57] THz is downshifted to the baseband by 193.4 THz. Upon implementing the proposed technique, we calculate the extrapolated spectra $S_L$ and $S_R$ using the maximum polynomial order of $M=K=9$ and an extended frequency range of $-f_L=f_R=0.31$ THz. The full spectrum ${S}_{full}$ and its smoothness is presented in Fig. 8 (a). Figure 8 (b) shows the impulse response directly calculated from ${S}_{full}$, with a maximum amplitude of 4.3e-5 for the non-causal part. By zeroing out the non-causal part, we obtain the causal spectrum $\bar {S}_{full}$ using the FFT, which is then compared to $S$ in Fig. 9(a) and (b). Meanwhile, the CVF model is also built with 48 poles after a total iteration of 16. As demonstrated in Fig. 9, the proposed model and the CVF model exhibit maximum absolute errors of -75 dB and -60 dB in magnitude, 0.037 rad and 0.09 rad in phase, respectively, when compared to the data $S$.

 figure: Fig. 7.

Fig. 7. The structure of the MZI-based demultiplexer.

Download Full Size | PDF

 figure: Fig. 8.

Fig. 8. (a) The full spectrum $S_{full}$ of the demultiplexer after extrapolation. (b) The impulse response computed from $S_{full}$.

Download Full Size | PDF

 figure: Fig. 9.

Fig. 9. Modeling accuracy for the magnitude and phase of the BBLS of the demultiplexer. (a)-(b) Magnitude and phase of the full causal spectrum $\bar {S}_{full}$ compared to that of the BBLS. (c)-(d) Magnitude and phase of the CVF model compared to that of the BBLS.

Download Full Size | PDF

To evaluate the time-domain performance, we consider a raised cosine pulse-shaped PAM4 signal as the input with a symbol rate of 80 Gbaud/s and a total length of 4000 symbols. To achieve accurate transient time simulation, the time steps should be much smaller than the highest frequency components of the input signal, and the frequency range of the simulated or measured BBLS should cover the spectrum of the input signal. In this example, the frequency range of the BBLS is 350 GHz and is larger than that of the PAM4 signal. However, the time step can be calculated as $\Delta t=1/(2f_R)=1.61$ ps and is too large to obtain accurate transient simulation results. Thanks to the proposed technique, zeros can be directly padded at the beginning and the end of ${S}_{full}$ to expand the maximum frequency range, for example to 3.3 THz, leading to a new impulse response with much smaller time steps $\Delta t=0.3$ ps. It is important to remark that smaller time steps result in more accurate results but lead to longer simulation times. For the time-domain simulation, the convolution is performed in SIMULINK with the impulse response and the input sequence. Meanwhile, the CVF model is also simulated with the same input. Time-domain simulation results from both methods are displayed in Fig. 10, showing an RMS error of 0.003 with $\Delta t=0.3$ ps. In contrast, an RMS error of 0.016 is obtained in the simulation with $\Delta t=1.61$ ps. Note that, only the time window from 0 ps to 600 ps is presented for clarity. To better observe that full output signal and the filtering effect, Fig. 11 illustrates the PAM4 eye diagrams of the input and computed output.

 figure: Fig. 10.

Fig. 10. The time-domain simulation results of the CVF model and the impulse response ($\Delta t=0.3$ ps) with a PAM4 input signal.

Download Full Size | PDF

 figure: Fig. 11.

Fig. 11. Eye diagrams of the PAM4 input signal (left) and the output (right) after filtering.

Download Full Size | PDF

In term of efficiency, it is worth mentioning that building the proposed model is 10.2 times faster compared to the CVF model. The time-domain simulation for the proposed method only takes 1.52 s with $\Delta t=0.3$ ps and 0.12 s with $\Delta t=1.61$ ps, whereas solving the CVF model consumes 116 s. This stark contrast illustrates that the proposed method is significantly more efficient than the CVF technique when simulating filters with long time sequences.

5. Conclusion

We have presented an accurate extrapolation scheme for extracting causal impulse responses from the baseband band-limited S-parameters of passive PICs. This method ensures a rapid and smooth decay of the extrapolated spectrum to zeros in the extended frequency range, allowing for zero-padding to generate new causal impulse responses. As a result, time steps in impulse responses can be customized to meet the needs of various applications. We have provided two examples to validate our proposed method. The technique is applicable to a variety of passive PICs and can be potentially integrated into PICs simulation tools for design automation.

Funding

National Natural Science Foundation of China (62165003, 62205075); Natural Science Research Project of Guizhou Province (ZK[2022]-136).

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. S. Shekhar, W. Bogaerts, L. Chrostowski, et al., “Silicon photonics – roadmapping the next generation,” arXiv, arXiv:2305.15820 (2023). [CrossRef]  

2. A. Narayan, Y. Thonnart, P. Vivet, et al., “System-level evaluation of chip-scale silicon photonic networks for emerging data-intensive applications,” in Proc. Design, Autom. Test Eur. Conf. Exhib. (DATE), (2020), pp. 1444–1449.

3. W. Bogaerts and L. Chrostowski, “Silicon photonics circuit design: methods, tools and challenges,” Laser Photonics Rev. 12(4), 1700237 (2018). [CrossRef]  

4. M. J. Shawon and V. Saxena, “Rapid simulation of photonic integrated circuits using Verilog-A compact models,” IEEE Trans. Circuits Syst. I 67(10), 3331–3341 (2020). [CrossRef]  

5. C. Sorace-Agaskar, J. Leu, M. R. Watts, et al., “Electro-optical co-simulation for integrated CMOS photonic circuits with VerilogA,” Opt. Express 23(21), 27180 (2015). [CrossRef]  

6. Y. Ye, D. Spina, D. Deschrijver, et al., “Time-domain compact macromodeling of linear photonic circuits via complex vector fitting,” Photonics Res. 7(7), 771–782 (2019). [CrossRef]  

7. B. Gustavsen and A. Semlyen, “Rational approximation of frequency domain responses by vector fitting,” IEEE Trans. Power Delivery 14(3), 1052–1061 (1999). [CrossRef]  

8. S. Mingaleev, A. Richter, E. Sokolov, et al., “Towards an automated design framework for large-scale photonic integrated circuits,” Proc. SPIE 9516, 951602 (2015). [CrossRef]  

9. J. Pond, C. Cone, L. Chrostowski, et al., “A complete design flow for silicon photonics,” Proc. SPIE 9133, 913310 (2014). [CrossRef]  

10. C. Yang, H.-D. Brüns, P. Liu, et al., “Impulse response optimization of band-limited frequency data for hybrid field-circuit simulation of large-scale energy-selective diode grids,” IEEE Trans. Electromagn. Compat. 58(4), 1072–1080 (2016). [CrossRef]  

11. F. Rao, “Optimization of spectrum extrapolation for causal impulse response calculation using the hilbert transform,” U.S. Patent No. 7 962 541 B2 (2011).

12. J. Cho, K. Hwang, S. Jeung, et al., “An efficient extrapolation method of band-limited S-parameters for extracting causal impulse responses,” IEEE Trans. Comput.-Aided Des. Integr. Circuits Syst. 38(11), 2086–2098 (2019). [CrossRef]  

13. J. Cho, J. Ahn, J. Kim, et al., “Low- and high-frequency extrapolation of band-limited frequency responses to extract delay causal time responses,” IEEE Trans. Electromagn. Compat. 63(3), 888–901 (2021). [CrossRef]  

14. J. Becerra, F. Vega, and F. Rachidi, “Extrapolation of a truncated spectrum with Hilbert transform for obtaining causal impulse responses,” IEEE Trans. Electromagn. Compat. 59(2), 454–460 (2017). [CrossRef]  

15. Y. Ye, T. Ullrick, W. Bogaerts, et al., “SPICE-compatible equivalent circuit models for accurate time-domain simulations of passive photonic integrated circuits,” J. Lightwave Technol. 40(24), 7856–7868 (2022). [CrossRef]  

16. M. C. Jeruchim, P. Balaban, and K. S. Shanmugan, Simulation of communication systems: modeling, methodology and techniques (Springer, New York, NY, USA, 2006), 2nd ed.

17. L. L. Barannyk, H. A. Aboutaleb, A. Elshabini, et al., “Spectrally accurate causality enforcement using SVD-based Fourier continuations for high-speed digital interconnects,” IEEE Trans. Compon., Packag., Manuf. Technol. 5(7), 991–1005 (2015). [CrossRef]  

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (11)

Fig. 1.
Fig. 1. Magnitude and phase of the baseband band-limited S-parameters (red solid lines) and the extrapolated spectrum (red dash lines) of passive PICs.
Fig. 2.
Fig. 2. (a) The full spectrum $S_{full}$ of the DC after extrapolation. (b) The impulse responses computed via $S_{full}$ and via directly padding zeros to the BBLS.
Fig. 3.
Fig. 3. Modeling accuracy for the magnitude and phase of the BBLS under analysis. (a)-(b) Magnitude and phase of the full causal spectrum $\bar {S}_{full}$ compared to that of the BBLS. (c)-(d) Magnitude and phase of the CVF model compared to that of the BBLS.
Fig. 4.
Fig. 4. Time-domain input and outputs for the proposed method and the CVF model.
Fig. 5.
Fig. 5. (a) The impulse response calculated by padding zeros to $S_{full}$. (b) The frequency spectrum of the new impulse response, compared to the BBLS.
Fig. 6.
Fig. 6. The time-domain simulation results of the CVF model and the impulse response with $\Delta t=0.02$ ps.
Fig. 7.
Fig. 7. The structure of the MZI-based demultiplexer.
Fig. 8.
Fig. 8. (a) The full spectrum $S_{full}$ of the demultiplexer after extrapolation. (b) The impulse response computed from $S_{full}$.
Fig. 9.
Fig. 9. Modeling accuracy for the magnitude and phase of the BBLS of the demultiplexer. (a)-(b) Magnitude and phase of the full causal spectrum $\bar {S}_{full}$ compared to that of the BBLS. (c)-(d) Magnitude and phase of the CVF model compared to that of the BBLS.
Fig. 10.
Fig. 10. The time-domain simulation results of the CVF model and the impulse response ($\Delta t=0.3$ ps) with a PAM4 input signal.
Fig. 11.
Fig. 11. Eye diagrams of the PAM4 input signal (left) and the output (right) after filtering.

Equations (16)

Equations on this page are rendered with MathJax. Learn more.

E ( t ) = A ( t ) cos ( 2 π f c t + ϕ ( t ) )
E b ( t ) = A ( t ) e j ϕ ( t ) ,
{ d x ( t ) d t = A x ( t ) + B a ( t ) b ( t ) = C x ( t ) + D a ( t ) ,
S ( f ) = 1 π p . v . S ( f ) f f d f , S ( f ) = 1 π p . v . S ( f ) f f d f ,
S L ( f i ) = m = 4 M a m ( f i f L ) m , f L f i < f m i n .
S R ( f i ) = k = 4 K b k ( f i f R ) k , f m a x < f i f R ,
h ( n ) = p = P / 2 P / 2 1 S f u l l ( p ) φ n p ,
p = P / 2 P / 2 1 S f u l l ( p ) φ n p = 0 , P / 2 n < 0.
p = P / 2 P / 2 + L N S L ( p ) φ n p + p = P / 2 + L N + 1 P / 2 + L N + N S ( p ) φ n p + p = P / 2 + L N + N + 1 P / 2 1 S R ( p ) φ n p = 0 , P / 2 n < 0
A L S L + A R S R = A S ,
A L = [ A P / 2 A P / 2 + L N ] C ( P / 2 ) × ( L N + 1 ) A = [ A P / 2 + L N + 1 A P / 2 + L N + N ] C ( P / 2 ) × N A R = [ A P / 2 R N A P / 2 1 ] C ( P / 2 ) × R N A i = [ φ P / 2 , i φ 1 , i ] T C ( P / 2 ) × 1
S L = F a , S R = G b ,
F = [ ( f P / 2 f L ) 4 ( f P / 2 f L ) M ( f P 2 + L N f L ) 4 ( f P 2 + L N f L ) M ] R ( L N + 1 ) × ( M 3 ) G = [ ( f P / 2 R N f R ) 4 ( f P / 2 R N f R ) K ( f P / 2 1 f R ) 4 ( f P / 2 1 f R ) K ] R R N × ( K 3 ) a = [ a 4 a M ] T C ( M 3 ) × 1 b = [ b 4 b K ] T C ( K 3 ) × 1
[ A L F A R G ] [ a b ] = A S .
[ A L F A R G F c 0 0 G c ] [ a b ] = [ A S S c ] ,
F c = [ ( f P 2 + L N + 1 f L ) 4 ( f P 2 + L N + 1 f L ) M ( f P 2 + L N + 4 f L ) 4 ( f P 2 + L N + 4 f L ) M ] R 4 × ( M 3 ) G c = [ ( f P / 2 R N 4 f R ) 4 ( f P / 2 R N 4 f R ) K ( f P / 2 R N 1 f R ) 4 ( f P / 2 R N 1 f R ) K ] R 4 × ( M 3 ) S c = [ S ( P / 2 + L N + 1 ) S ( P / 2 + L N + 4 ) S ( P / 2 R N 4 ) S ( P / 2 R N 1 ) ] T
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.