## Abstract

In this paper, a novel modeling and simulation method for general linear, time-invariant, passive photonic devices and circuits is proposed. This technique, starting from the scattering parameters of the photonic system under study, builds a baseband equivalent state-space model that splits the optical carrier frequency and operates at baseband, thereby significantly reducing the modeling and simulation complexity without losing accuracy. Indeed, it is possible to analytically reconstruct the port signals of the photonic system under study starting from the time-domain simulation of the corresponding baseband equivalent model. However, such equivalent models are complex-valued systems and, in this scenario, the conventional passivity constraints are not applicable anymore. Hence, the passivity constraints for scattering parameters and state-space models of baseband equivalent systems are presented, which are essential for time-domain simulations. Three suitable examples demonstrate the feasibility, accuracy, and efficiency of the proposed method.

© 2018 Chinese Laser Press

## 1. INTRODUCTION

In recent years, silicon photonic devices and circuits have had rapid development both in complexity and functionality thanks to an increasingly mature manufacturing process. At the same time, several computer-aided design (CAD) tools have emerged both in academic and industrial areas to analyze the behavior of silicon photonic designs.

Time-domain simulation of integrated photonic circuits is an essential part in the design flow, since it gives the most intuitive assessment of system performance. For some basic photonic components, such as waveguides, time-domain simulations can be analytically addressed because of the simple underlying physical principles and equations. However, time-domain simulations cannot be performed analytically when considering more complex devices or effects caused by parasitic elements. In such scenarios, different simulation techniques, such as finite-difference time-domain (FDTD) [1], time-domain traveling wave (TDTW) [2,3], the split-step method (SSM) [4], coupled mode theory (CMT) [5], or convolution-based methods [6], operate on the component or circuit level. However, a trade-off between efficiency and accuracy has to be made when using these techniques [5].

In the electronic field, a similar problem holds for distributed devices, such as nonuniform transmission lines or microstrip filters, since no compact circuit models are readily available for time-domain simulations [7]. A popular solution is based on a frequency-domain system identification technique named the vector fitting (VF) algorithm [8], which is able to build stable and passive rational models of the scattering parameters of the devices under study. Then, these frequency-domain models can be directly converted to an equivalent state-space representation in the time domain. This technique is widely applied for electronic systems, for example in Refs. [8–13].

Since the VF method is developed for linear, time-invariant, passive systems and is based on their transfer function representation (e.g., scattering parameters), it is immediately applicable to passive photonic devices and circuits [14]. However, compared to electronic systems, the frequency range of interest for photonic systems is typically around [187,200] THz, corresponding to a wavelength range of [1.5,1.6] μm. Such a wide range at high frequencies has a direct impact on the modeling and simulation processes, which can become very time and/or memory consuming.

To address this problem, a novel modeling method is proposed in this paper, which is based on baseband equivalent signal and system representation. In particular, the proposed modeling approach computes an accurate baseband equivalent state-space representation, starting from the scattering parameters of the photonic system under study evaluated at optical frequencies. However, such an equivalent state-space model is complex-valued and not physically realizable. Furthermore, the stability and passivity constraints on scattering parameters and state-space models of complex-valued systems, which are fundamental properties for time-domain simulations, appear yet to be missing in the literature. In this paper, we rigorously discuss these conditions for the proposed baseband equivalent systems based on the classic definitions of stability and passivity to validate the feasibility of the proposed time-domain simulation method. The proposed technique offers two main advantages: (1) the modeling process is based on the scattering parameters, which makes it a widely applicable method for generic linear passive photonic components and circuits; (2) the state-space representation is a continuous time-domain model, which can be efficiently simulated in the time domain without involving convolution, fast Fourier transform (FFT), or inverse fast Fourier transform (IFFT), thereby making this method robust and accurate.

The paper is organized as follows. Section 2 presents an overview of the “standard” modeling approach based on the VF algorithm, while Section 3 introduces the baseband equivalent signals and systems and describes the novel proposed modeling framework. The stability and passivity constraints of such systems are discussed in Section 4. A practical guideline for the application of the proposed modeling approach is given in Section 5, while Section 6 validates the proposed method by means of three pertinent numerical examples. Conclusions are drawn in Section 7.

## 2. CONVENTIONAL STATE-SPACE MODELING OF PHOTONIC SYSTEMS

In both electronics and photonics, the scattering matrix is widely used to describe the behaviors of passive devices and circuits

where $s$ is the Laplace variable, $\mathit{a}(s)$ and $\mathit{b}(s)$ are the forward wave and backward wave, respectively, and $\mathit{S}(s)$ is the scattering matrix of the system under study, which can be obtained through simulations or measurements. The aim of the rational modeling is to find a Laplace-domain model of Eq. (1) in a pole-residue form as where $\mathit{D}\in {\mathbb{R}}^{n\times n},{\mathit{R}}_{k}\in {\mathbb{C}}^{n\times n},k=1,\dots ,K$, $n$ and $K$ being the number of ports of the system under study and the number of poles used to approximate the scattering parameters, respectively. Typically, all the elements ${S}_{ij}(s)$ of the scattering matrix representation in Eq. (2) use a common denominator polynomial and pole set $[{p}_{1},{p}_{2},\dots ,{p}_{K}]$, where such poles are either real quantities or complex conjugate pairs [8]. The identification of poles ${p}_{k}$ and residue matrices ${\mathit{R}}_{k}$ can be performed via the VF algorithm [8,15–18], starting from a set of the scattering parameters under study obtained for ${s}_{r}=j2\pi {f}_{r}$ with $r=1,\dots ,R$.However, it is important to note that the sign convention ${e}^{jwt}$ is commonly used in the electronics field to represent the incident and reflected waves in Eq. (1), while ${e}^{-jwt}$ is sometimes adopted in the optics field [19,20]. Hence, the scattering matrix defined with one sign convention is the complex conjugate of the other. The VF algorithm is based on the assumption that the sign convention ${e}^{jwt}$ is adopted, since it has been originally developed for electromagnetic problems. In the case when ${e}^{-jwt}$ is used to define the scattering parameters under study, a simple solution is to apply the VF algorithm to the complex conjugate of the scattering matrix.

Then, the rational model in Eq. (2) can be transformed to state-space form by a simple rearrangement [17,21]

where $\mathit{A}\in {\mathbb{C}}^{m\times m}$, $\mathit{B}\in {\mathbb{R}}^{m\times n}$, $\mathit{C}\in {\mathbb{C}}^{n\times m}$, $\mathit{D}\in {\mathbb{R}}^{n\times n}$, $m=nK$, and $\mathit{I}$ is the identity matrix in this paper. In particular, $\mathit{A}$ is a diagonal matrix with all the poles as diagonal elements, while $\mathit{C}$ contains all the residues, and they can be always converted to real matrices as long as the poles and residues are real or complex conjugate pairs [17].Now, it is straightforward to convert Eq. (3) to an equivalent state-space representation in the time domain [21] as

Note that fundamental properties for time-domain simulations, such as the stability and passivity of models in Eq. (4), must be assured [13]. While the stability of VF models can be guaranteed by construction by means of suitable pole-flipping schemes [8], their passivity can be checked and, eventually, enforced only after the rational model is computed by adopting passivity enforcement techniques. Indeed, due to the unavoidable numerical approximations, the rational model computed might be non-passive. Several robust passivity enforcement methods have been proposed in the literature; see for example Refs. [16–18]. Now, time-domain simulations can be carried out by solving the first-order system of ordinary differential equations (ODE) in Eq. (4) via suitable numerical techniques [22,23]. These approaches iteratively solve Eq. (4) for a discrete set of values of the time, which are chosen via suitable algorithms (i.e., fixed or adaptive time-step). In particular, the computational cost of solving Eq. (4) depends on three main elements:

- • the bandwidth of the signals considered, which define the maximum time-step $\mathrm{\Delta}{t}_{\text{max}}$ that can be adopted: $\mathrm{\Delta}{t}_{\text{max}}$ must be smaller than the highest frequency component of the signals considered;
- • the numerical technique adopted to solve Eq. (4);
- • the number of poles $K$ and of ports $n$ of the system under study, which directly determine the number of states $m=nK$.

The modeling technique described so far allows one to simulate any generic linear and passive system in the time or frequency domain, and it has found extensive applications in electronic engineering problems [8–13]. However, when it comes to photonic circuits, one substantial difference arises with respect to the electronic domain: the range of frequency of interest is typically around [187, 200] THz, corresponding to a wavelength of ${U}_{l}(f)$, or even higher frequencies for shorter wavelengths. This has a major impact on the modeling and simulation complexity of the approach described so far. Indeed, a large number of poles $K$ can be required to accurately model the scattering parameters in the chosen frequency range, and the passivity enforcement phase can become computationally prohibitive. Furthermore, the corresponding ODE in Eq. (4) will have a large number of equations, and a small time-step (of the order of femtoseconds) must be adopted to solve it.

In order to tackle these issues, a novel approach based on baseband equivalent state-space models is proposed in this contribution.

## 3. BASEBAND EQUIVALENT STATE-SPACE MODELS FOR THE TIME-DOMAIN SIMULATION OF PHOTONIC SYSTEMS

The basic concepts of baseband equivalent signals and systems are first introduced in Section 3.A, given their importance in the definition of the proposed modeling approach, which is described in Section 3.B.

#### A. Baseband Equivalent Signals and Systems

The excitation signal of photonic systems is often an amplitude and/or phase-modulated signal with optical carrier and electronic modulating signals, which can be written in the following form:

where $A(t)$ is the time-varying amplitude or envelope of the modulated signal, and $\varphi (t)$ is the time-varying phase. In electronics or radio-frequency (RF) applications, both $A(t)$ and $\varphi (t)$ relate to electronic signals, such as voltage, current, or electric field. In photonics, the optical carrier frequency ${f}_{c}$ is much higher than the bandwidth of $A(t)$ and $\varphi (t)$, given that the wavelength of light is much smaller than that of RF signals, so the representation in Eq. (5) is often called a bandpass signal.An analytic complex-valued representation of the real-valued signal in Eq. (5), called the analytic signal, is introduced here as [24]

where $\mathcal{H}[u(t)]$ is the Hilbert transform of $u(t)$. In the frequency domain, Eq. (6) becomes where ${U}_{a}(f)$ and $U(f)$ are the Fourier transform of ${u}_{a}(t)$ and $u(t)$, respectively, and $\mathrm{Step}(f)$ is a unit step function defined byNow, the corresponding baseband equivalent signal of the bandpass signal is defined as

which can be considered as the complex envelope optical signal representation and is widely used in photonics and optical fiber communication. The relations between $u(t)$, $\mathcal{H}[u(t)]$, and ${u}_{l}(t)$ in the time and frequency domains are [24] where the superscript * denotes the complex conjugate operator.In the frequency domain, the concepts of analytic signal and baseband equivalent signal are intuitive: $U(f)$ has a symmetric spectrum with respect to the positive and negative frequencies, while ${U}_{a}(f)$ has only a non-zero spectrum in the positive frequencies around the carrier frequency; shifting the spectrum of ${U}_{a}(f)$ in the direction of the negative frequencies of ${f}_{c}$ [or equivalently in the time domain by multiplying ${u}_{a}(t)$ with ${e}^{-j2\pi {f}_{c}t}$] leads to ${U}_{l}(f)$. Such relations are illustrated in Fig. 1.

If a system with impulse response $h(t)$ and frequency response $H(f)$ operates in the bandwidth BW around ${f}_{c}$, satisfying ${f}_{c}\gg \mathrm{BW}$, then it can be considered as a bandpass system. Now, the corresponding baseband equivalent system can be defined by applying the same concepts described for the baseband signals. Thanks to the relations among bandpass signals and systems and their baseband equivalents, it can be proven that the output signal of a bandpass system can be analytically recovered from the output of the corresponding baseband system, as illustrated in Fig. 2. A detailed proof is given in Appendix A.

It is important to remark that performing time-domain simulations of baseband equivalent systems allows one to efficiently recover the corresponding bandpass signals, thus avoiding the expensive time-domain simulations of photonic systems at optical frequencies.

#### B. Realization of Baseband Equivalent Signals and Systems

Baseband equivalent signals ${u}_{l}(t)$ can be easily computed with Eq. (9), where ${u}_{l}(t)$ can be a real (amplitude modulation) or a complex signal (when both amplitude and phase modulation are applied). For example, in the case of a quadrature amplitude modulation (QAM), ${u}_{l}(t)$ can be expressed with respect to its in-phase component $I(t)=A(t)\mathrm{cos}\text{\hspace{0.17em}}\varphi (t)$ and quadrature component $Q(t)=A(t)\mathrm{sin}\text{\hspace{0.17em}}\varphi (t)$ as ${u}_{l}(t)=I(t)+jQ(t)$.

Note that baseband equivalent signals and systems are widely used in the simulation of communication systems to simplify the modulation, demodulation, and filtering process [24]. In such a scenario, continuous systems and signals are often first sampled and defined as finite discrete sequences, and then convolution, FFT, or IFFT are adopted for the time-domain simulation of the discrete-time representations of such signals and systems, which could lead to inaccurate results [6].

In this section, the goal is to build stable and passive continuous models for baseband equivalent systems in state-space form, whose time-domain simulation can also capture transient behaviors. However, a readily baseband counterpart for Eq. (4) does not exist in the literature. Indeed, baseband systems have an asymmetric frequency response with respect to the positive and negative frequencies [similar to ${U}_{l}(f)$ in Fig. 1] resulting in a non-physical, complex valued system, as described in detail in Appendix A. The VF algorithm [8] is a technique developed for physical systems with a symmetric frequency response, which can be described with real or complex conjugate poles: this situation clearly does not hold for baseband systems, and VF cannot be directly applied to the baseband response of the system under study.

In order to reach our goal, the first step is to express $\mathit{a}(t)$, $\mathit{b}(t)$, and $\mathit{x}(t)$ in the system of ODE from Eq. (4) in the form from Eq. (11), which gives

Equations (14) and (15) allow us to write

After simple mathematical manipulations, Eq. (16) can be written as

It is important to remark upon one difference between the representation in Eq. (17) and the definition of baseband systems: in Eq. (17), the entire frequency response of the system under study is shifted in the baseband, while for the baseband system introduced in Section 3.A, only the frequency response at positive frequencies is shifted in the baseband. However, in Appendix B it is rigorously proven that these two representations are equivalent in terms of time-domain simulations. Hence, in the rest of the contribution, the expression “baseband equivalent system” does not refer to the classic definition given in Section 3.A and Appendix A, but to the new proposed baseband equivalent “shifted” system, where the entire frequency response of the system under study is shifted in the baseband.

A similar realization of baseband equivalent systems in the frequency domain, computed by shifting the poles of the transfer function of the corresponding bandpass system by $j2\pi {f}_{c}$, has been presented in the electronic domain in Refs. [24,25], but the derivation is not given. Note that the time-domain simulation methods in Refs. [24,25] are substantially different from the one presented here. In Ref. [24], once the transfer function of the baseband equivalent system is obtained, it is first sampled and converted to an equivalent discrete system, and then the discrete impulse response is calculated via IFFT. Finally, the time-domain behavior of the baseband equivalent system is simulated by convolution. In Ref. [25], first the inverse Laplace transform is adopted to analytically convert the baseband equivalent transfer function to a continuous impulse response, then a recursive convolution technique is used to perform time-domain simulations. In contrast, the time-domain simulation method presented in this paper directly solves the corresponding ODE, which is more straightforward. However, it is crucial to prove that fundamental properties for time-domain simulations, such as stability and passivity [13], still hold for the proposed baseband equivalent state-space representation.

## 4. PASSIVITY OF THE BASEBAND EQUIVALENT SYSTEM

The poles and residues of rational models of electronic and photonic systems are always real or complex conjugate pairs, as discussed in Section 2. However, for the baseband equivalent state-space model in Eq. (17), the poles do not follow this rule anymore; furthermore, the corresponding frequency response is not symmetrical with respect to positive and negative frequencies, which makes the baseband equivalent system a non-physical, complex-valued system. Finally, the most remarkable difference with respect to bandpass systems is that the impulse response of these baseband equivalent systems is not real, and with a real input, they can generate a complex output.

Then, it is important to verify if such linear, time-invariant complex-valued systems still comply with the passivity conditions of “conventional” real-valued systems, which are listed as follows [26]:

- (1) Each element of $\mathit{S}(s)$ is analytic in $\mathrm{Re}(s)>0$;
- (2) $\mathit{I}-{\mathit{S}}^{H}(s)\mathit{S}(s)$ is a nonnegative-definite matrix for all $s$ such that $\mathrm{Re}(s)>0$;
- (3) ${\mathit{S}}^{*}(s)=\mathit{S}({s}^{*})$.

The ${\mathrm{superscript}\text{\hspace{0.17em}}\text{\hspace{0.17em}}}^{H}$ stands for the transpose conjugate operator. The first condition is related to causality and stability; the second one is basically a bound for $\mathit{S}(s)$; the third ensures that the associated impulse response is real, which requires the system to be real-valued [27]. Evidently, the third condition is not suitable for complex-valued systems anymore. In this section, the passivity constraints for scattering parameters of baseband equivalent systems will be proposed, and a fast assessment of the passivity of the corresponding baseband equivalent state-space model will be presented.

#### A. Passivity Constraints on Scattering Parameters of Baseband Equivalent Systems

According to Refs. [26,28,29], an $n$-port electronic system is passive if, for any $\tau >-\infty $ and $\mathit{v}(t)\in {L}_{2n}$ (${L}_{2n}$ denotes the space of all vectors whose $n$ components are functions of a real variable $t$ and square integrable over $-\infty <t<\infty $), it holds

where $\mathit{v}(t)$ and $\mathit{i}(t)$ are the voltage and current at the system ports, respectivle. It is important to note that this definition is given not only for real signals but also for complex ones. By expressing Eq. (18) in terms of the forward $\mathit{a}(t)$ and backward $\mathit{b}(t)$ waves, the passivity definition becomes [26,30,31]Following the same proof process as Ref. [26], particularly Theorem 2 and Theorem 3, the first and second passivity conditions can be derived from Eq. (19) for the complex-valued systems studied in this paper. Alternatively, the same conclusion can be obtained via the approach in Chapter II of Ref. [30], which gave simpler formal proofs using the theory of distributions. The interested reader may consult Refs. [26,30] for a detailed and comprehensive proof.

Therefore, in this paper, we propose the following passivity constraints on the scattering parameters ${\widehat{\mathit{S}}}_{l}(s)$ of the baseband equivalent systems:

- (1) ${\widehat{\mathit{S}}}_{l}(s)$ is analytic in $\mathrm{Re}(s)>0$;
- (2) $\mathit{I}-{\widehat{\mathit{S}}}_{l}^{H}(s){\widehat{\mathit{S}}}_{l}(s)$ is a nonnegative-definite matrix for all $s$ such that $\mathrm{Re}(s)>0$.

Note that real-valued systems need the extra condition $\mathit{S}({s}^{*})={\mathit{S}}^{*}(s)$, which ensures that the impulse response is real, so that a real input results in a real output and makes the system physically realizable. Furthermore, it is clearly mentioned in Section 4 of Ref. [26] that this requirement is independent with respect to the passivity definition in Eqs. (18) and (19). Therefore, this is evidently not required for the passivity of complex-valued systems, which are proposed only for simulation purposes.

#### B. Fast Passivity Assessment of Baseband Equivalent Systems

Passivity conditions require both scattering parameters $\mathit{S}(s)$ and their baseband equivalent ${\widehat{\mathit{S}}}_{l}(s)$ to be bounded by unity, which implies that all singular values $\sigma $ of ${\widehat{\mathit{S}}}_{l}(s)$ are smaller than unity at all frequencies

An efficient and accurate method to assess the passivity properties of state-space models of electronic and photonic systems is based on the Hamiltonian matrix $\mathit{M}$ [17], defined as

A state-space model is passive if its Hamiltonian matrix has no purely imaginary eigenvalues, since any imaginary eigenvalue indicates a crossover frequency where a singular value changes from being smaller to larger than unity, or vice versa. This approach is more reliable and efficient than sweeping the singular values over a set of discrete frequencies, especially for photonic systems, which are defined over a large frequency range.

A similar Hamiltonian matrix ${\widehat{\mathit{M}}}_{l}$ for baseband equivalent systems ${\widehat{\mathit{S}}}_{l}(s)$ can be derived by following the procedure in Ref. [17], leading to

One can observe that the only difference between $\mathit{M}$ and ${\widehat{\mathit{M}}}_{l}$ is the use of the transpose conjugate operator for the state-space matrices in ${\mathit{M}}_{l}$, while only the transpose operator is required in $\mathit{M}$. Indeed, state-space models of general electronic or photonic systems satisfy the conjugacy property ${\mathit{S}}^{*}(s)=\mathit{S}({s}^{*})$; the corresponding scattering parameters do not change if the state-space matrices $\mathit{A}$, $\mathit{B}$, $\mathit{C}$, $\mathit{D}$ are replaced with their conjugate counterparts [17]. Evidently, this is not valid for the baseband equivalent systems.

Note that the eigenvalues of Eq. (22) can be obtained directly from the eigenvalues of the corresponding bandpass system from Eq. (21). According to Eq. (17), replacing ${\widehat{\mathit{A}}}_{l}$, ${\widehat{\mathit{B}}}_{l}$, ${\widehat{\mathit{C}}}_{l}$, ${\widehat{\mathit{D}}}_{l}$ in Eq. (22) with

Hence, the following properties hold:

- • If there are passivity violations in a bandpass state-space model, the corresponding baseband equivalent system from Eq. (17) is not passive either.
- • There is an one-to-one correspondence between the frequencies where passivity violations occur in the state-space models of the bandpass and corresponding baseband equivalent.

## 5. PROPOSED MODELING FRAMEWORK OF THE PHOTONIC SYSTEM FOR TIME-DOMAIN SIMULATIONS

The signals traveling through photonic systems are usually phase- and/or amplitude-modulated signals over a suitable optical carrier. The modulating signals are electronic ones, whose spectrum bandwidth is normally less than a few hundred gigahertz, while the carrier frequency is usually defined in the range [187.5,200] THz, corresponding to a wavelength range of [1.5,1.6] μm.

The proposed modeling approach starts from evaluating the scattering parameters of the photonic system under study in the frequency range of interest. Next, an accurate rational model is computed via the VF algorithm. Stability is enforced during the model-building phase via suitable pole-flipping schemes [8], while the model passivity is checked and, eventually, enforced as a post-processing step via robust passivity enforcement methods [18,32]. A baseband equivalent state-space representation from Eq. (17) can now be obtained with guaranteed passivity by Eq. (24). Such a model can be used to efficiently perform time-domain simulations. The flow chart of the proposed modeling framework is shown in Fig. 3.

In particular, when it comes to building state-space models of photonic systems for time-domain simulations, there are two options:

- (1) modeling the frequency range of interest, e.g., [187.5,200] THz, noted as Model A (covering a large frequency range);
- (2) considering only the frequency range corresponding to the spectrum of the optical input signals under study around the carrier frequency, normally a few hundred gigahertz, noted as Model B (covering a small frequency range).

The corresponding baseband equivalent state-space models are indicated as Model LA and Model LB, respectively. The modeling frequency ranges of these four models are illustrated in Fig. 4 when assuming that ${f}_{c}=193\text{\hspace{0.17em}}\mathrm{THz}$ and the spectrum of the optical input signal of interest is 300 GHz. Note that Model A and Model B can also be used directly to evaluate the behavior of the chosen photonic system in the time domain: such modeling strategies follow the approach outlined in Section 2.

Note that Models A and LA are likely to require more poles as compared to Models B and LB, since they are computed over a larger bandwidth; the modeling complexity is higher and the corresponding system of ODE will be larger. If the scattering parameters under study are very dynamic in the range [187.5,200] THz, the modeling process can become prohibitively expensive, making it practically infeasible to build accurate, stable, and passive models. However, this approach offers more flexibility since the corresponding models can be used for any value of the carrier frequency in the frequency range [187.5,200] THz, while Models B and LB must be constructed anew for each value of the carrier frequency considered.

It is important to note that both Models LA and LB operate at baseband, which means that a relatively large time step can be used to solve the corresponding ODE for time-domain simulation, thereby saving both computational time and memory storage. Table 1 compares the advantages and disadvantages of different approaches considered during the model-building and time-domain simulation process.

Finally, no matter which model is used for the time-domain simulation, the modeling frequency range must be larger than or at least equal to the frequency range of the input signals considered. Indeed, no information on the scattering parameters’ behavior outside such a modeling frequency range is provided to the VF algorithm: the model obtained via the VF approach extrapolates the scattering parameters outside the modeling frequency range. Hence, while the state-space model computed is stable and passive at $[0,\infty ]\text{\hspace{0.17em}}\mathrm{Hz}$, it is not possible to guarantee its accuracy outside the modeling frequency range. Therefore, if the input signal is noisy, the spectrum of the noise should also be considered during the model-building phase.

## 6. NUMERICAL EXAMPLE

This section presents three application examples of the proposed modeling and simulation technique. The scattering parameters of the photonic systems under study are evaluated via Caphe [33], while the time-domain simulations are carried out in MATLAB [34] via the routine lsim on a personal computer with Intel Core i3 processor and 8 GB RAM.

#### A. Mach–Zehnder Interferometer

In this example, the Mach–Zehnder interferometer (MZI) shown in Fig. 5 is studied, which is formed by two identical directional couplers (with coupling coefficient 50/50) and two waveguides with lengths 150 μm (upper one) and 100 μm (lower one). Both waveguides have effective index 2.35 and group index 4.3 at 1.55 μm and a propagation loss of 200 dB/m. The time-domain simulation is carried out with the conventional modeling technique (in Section 2) and the proposed baseband equivalent modeling approach. For comparison, an analytic model for the MZI is also built by considering the loss and dispersion of the waveguides. The directional coupler is assumed to be an ideal signal splitter or combiner, which adds a $\pi /2$ phase delay to the cross-coupled signals. The time-domain simulation of this analytical model is conducted as a benchmark.

The modulating signal is a smooth pulse with amplitude 1 V, a rise/fall time of 5.7 ps, width of 32 ps, initial delay of 18 ps, and a spectrum bandwidth of 100 GHz. An optical carrier of frequency ${f}_{c}=193.72\text{\hspace{0.17em}}\mathrm{THz}$, which is chosen at random in the frequency range [187.5,200] THz, is used to transmit the modulated signal through the MZI. Both the modulated signal at optical frequencies and the smooth pulse are shown in Fig. 6.

Model A is built starting from the MZI scattering parameters in the range [187.5,200] THz, while Model B requires only the scattering parameters in $[{f}_{c}-\mathit{\Delta},{f}_{c}+\mathit{\Delta}]$, where the choice $\mathit{\Delta}=150\text{\hspace{0.17em}}\mathrm{GHz}$ allows one to cover the entire spectrum of the modulated optical signal. In particular, first the frequency samples have been divided in two groups: one to compute the desired rational model (modeling data) and the other to verify its accuracy (validation data). Then Models A and B are built via the VF algorithm with 67 poles and 8 poles, respectively, aiming at a maximum absolute error of less than $-60\text{\hspace{0.17em}}\mathrm{dB}$ between the model and MZI scattering parameters. Finally, Models LA and LB can be derived analytically from Models A and B, as shown in Section 3.B. The accuracy of Models A and LB in the frequency-domain is shown in Figs. 7 and 8, respectively, which show both the magnitude and the phase of the MZI scattering parameters obtained by Caphe and by the corresponding state-space models.

Time-domain simulations are carried out with all four models considered; while for Models A and B a time step of 0.23 fs is adopted, a time step of 0.4 ps can be used for Models LA and LB. Meanwhile, time-domain simulation of the analytic model built according to the underlying physical principle of the MZI is performed in Caphe to validate the accuracy of the other models. The outputs at port P3 of Model A, Model LB, and the analytic model are shown in Fig. 9. According to Section 3, the magnitude of the outputs of Model LB is the envelope of the output of Model A, and this fact is exactly illustrated by Fig. 9. In addition, it is easy to observe that the output of Model LB accurately matches the analytic model prediction.

The time for model building and time-domain simulation for all the different models is presented in Table 2. It clearly shows that modeling only the small frequency range (Models B and LB) rather than the large frequency range (Models A and LA) consumes far less time and results in compact models. Note that the time-domain simulation at baseband with compact models, such as Model LB, is the most efficient, which is consistent with the analysis in Section 5.

Finally, the following test illustrates the importance of choosing the correct modeling frequency range, as mentioned in Section 5. Let us assume an electronic pulse signal with width of 1 ps and spectrum in the range [0,6] THz as the input signal of Models LA and LB of the MZI. The corresponding output at port P3 is shown in Fig. 10: Model LA still gives very accurate results compared to the analytic model, while the output of Model LB is not even close to the benchmark. The reason is that the modeling frequency range (12.5 THz) of Model LA covers the spectrum of the input signal, but this does not hold for Model LB.

#### B. Ring Resonator

In this example, a double ring resonator (RR) is composed of two rings and two waveguides and designed as a narrowband flat-top filter, as shown in Fig. 11. The two rings have different circumferences of 20 μm (lower one) and 20.01 μm (upper one), resulting in slightly different ${R}_{1}$ and ${R}_{2}$. The ring waveguides and bus waveguides have effective index 2.35 and group index 4.3 at wavelength 1.55 μm. The coupling coefficient between waveguides and rings is 0.2, while the same parameter between two rings is 0.03.

First, Model A of the ring resonator is built in the range [187.5,200] THz with 22 poles, while Model B is computed with 6 poles in the range $[{f}_{c}-\mathit{\Delta},{f}_{c}+\mathit{\Delta}]$, with ${f}_{c}=195.75\text{\hspace{0.17em}}\mathrm{THz}$ and $\mathit{\Delta}=450\text{\hspace{0.17em}}\mathrm{GHz}$. The maximum absolute error of both models is less than $-65\text{\hspace{0.17em}}\mathrm{dB}$. Next, Model LB is directly derived by shifting the poles of Model B. Figures 12 and 13 describe the frequency-domain accuracy for Models A and LB, respectively.

In this example, a 4-QAM (quadrature phase-shift keying) modulated input signal is used for time-domain simulations. The in-phase I and quadrature Q parts of the modulating signal are the 4-bit sequences ($-1$, $-1$, 1, 1) and ($-1$, 1, $-1$, 1), respectively, where each bit lasts for 20 ps. As shown in Fig. 14, the modulating signals are realistic analog signals, for example, affected by overshoot and undershoot. As mentioned in Section 3.B, the baseband equivalent of the modulated input signal can be easily calculated, since I and Q are its real and imaginary parts, respectively.

After conducting the proposed time-domain simulation, the outputs of Model LB are complex, and their magnitude are the envelopes of the outputs of Model A as shown in Fig. 15. Note that the outputs of Model A can be analytically recovered from the outputs of Model LB, according to Eq. (B3). Hence, Fig. 16 shows a side-by-side comparison of the output of Model A at port P4 and the corresponding value recovered from Model LB. For a better observation of the accuracy of the recovered signal, Fig. 17 shows a zoom of Fig. 16 around $t=45.6\text{\hspace{0.17em}}\mathrm{ps}$, which demonstrates an excellent agreement.

As far as the computational time is concerned, building the Models A and LB required 0.28 s and 0.04 s, respectively, while their time-domain simulations took 9.29 s and 0.05 s, respectively, which clearly demonstrates the superior efficiency of the proposed technique when dealing with amplitude- and phase-modulated signals.

#### C. Lattice Filter

A fifth-order filter with a Chebyshev window, designed by using a discrete finite impulse response (FIR) filter design method [35], is realized via an MZI lattice filter (LF) [36]. As illustrated in Fig. 18, it is formed by six directional couplers with power coupling coefficients of 0.008, 0.067, 0.175, 0.175, 0.067, and 0.008, and waveguides with a length difference of 179 μm between the upper and lower ones, whose effective and group indices are 2.30 and 4.18, respectively. In practice, due to process variations, when manufacturing photonic devices, geometrical or optical parameters can vary in a relatively small range around their nominal value [37], which in turn can lead to variations in the device frequency response, such as frequency shifts. In this example, we study the time-domain influence of frequency shifts in the response of the lattice filter via an eye diagram analysis.

For eye diagram analysis, the input signal and time-domain simulation should last a relatively long period of time (long bits sequence), which could make the time-domain simulation of Models A and B unfeasible. In this example, a pseudo-random sequence of 1000 bits with a bit time of 30 ps and a Gaussian jitter having a standard deviation of 1.5 ps is used as modulating signal $A(t)$. The amplitude of the signal up to 1 ns is shown in Fig. 19. The total number of time steps required for time-domain simulations of Models A and B with such an input signal is 60 million (30 ns/0.5 fs), while this number reduces to only 30000 time-steps (30 ns/1 ps) for Models LA and LB.

The scattering matrices of the lattice filter are computed in the range [187.5,200] THz. However, due to the dynamic behavior of the filter frequency response in such a wide bandwidth, the modeling complexity of Model A (LA) is very high. Considering that the efficiency and accuracy of Model LB have been already demonstrated in Sections 6.A and 6.B, only the time-domain simulation of Model LB is performed.

The sequence signal is modulated on ${f}_{c}=195.11\text{\hspace{0.17em}}\mathrm{THz}$ ($\lambda =1.5365\text{\hspace{0.17em}}\mathrm{\mu m}$), which is chosen as the filter passband center frequency during the design phase. Due to manufacturing tolerances, let us assume that the center frequency can shift to 195.05 THz ($\lambda =1.5370\text{\hspace{0.17em}}\mathrm{\mu m}$) or 194.98 THz ($\lambda =1.5375\text{\hspace{0.17em}}\mathrm{\mu m}$), as shown in Fig. 20. Model LB is built for each one of these three situations by adopting a pole shift of ${f}_{c}=195.11\text{\hspace{0.17em}}\mathrm{THz}$, since the excitation signal is modulated on this frequency. In particular, the models for the three wavelengths considered are built with 36 poles, achieving a maximum absolute error of $-60\text{\hspace{0.17em}}\mathrm{dB}$.

Then, the time-domain simulations can be easily carried out at the baseband with the pseudo-random sequence of 1000 bits. Figure 21 shows the eye diagram of the power of the complex output signals at port P4 of the three baseband equivalent systems over a two-bit span resulting from the entire 1000-bit input stream. It is evident that the signal is completely distorted when the center frequency shifts from 195.11 to 194.98 THz. The computational time of the time-domain simulation for generating each eye diagram is 1.09 s, while building each model took 1.67 s, which is very efficient. This example shows that expensive time-domain simulations can be efficiently performed via the proposed technique without a loss in accuracy.

## 7. CONCLUSION

A novel modeling and simulation technique for passive photonic components and circuits that is flexible, efficient, accurate, and robust has been proposed in this paper,. Photonics systems can be characterized by the proposed baseband equivalent state-space models via the robust VF algorithm, which allows for the time-domain simulations to be conducted at the baseband rather than at the optical carrier frequency. The outputs of photonic systems can be immediately recovered from the outputs of the corresponding baseband equivalent models, thereby significantly decreasing the simulation time and memory usage. The passivity conditions of the proposed baseband equivalent systems are rigorously discussed, and a fast passivity assessment method for the corresponding state-space models is presented in this paper. The accuracy and efficiency of the proposed approach are verified by three suitable numerical examples.

## APPENDIX A: TIME-DOMAIN SIMULATION OF BASEBAND EQUIVALENT SIGNALS AND SYSTEMS

If a system with impulse response $h(t)$ and frequency response $H(f)$ operates in the bandwidth BW around ${f}_{c}$, satisfying ${f}_{c}\gg \mathrm{BW}$, then it can be considered as a bandpass system. Now, in a similar manner as with the baseband equivalent signal, a baseband equivalent system with impulse response ${h}_{l}(t)$ and frequency response ${H}_{l}(f)$ can be defined as [24]

where ${h}_{a}(t)$ is the analytic signal of $h(t)$ and is defined in the same way as Eq. (6).Compared with the definition of baseband equivalent signals, a factor 1/2 is introduced into the definition of baseband equivalent systems [24]. Again, the relations between $h(t)$, $\mathcal{H}[h(t)]$, and ${h}_{l}(t)$ in the time and frequency domains are

It is important to note that baseband equivalent signals and systems are not physical, but constitute a mathematical representation developed only for simplifying analysis and simulation of bandpass signals and systems, as discussed in the following.

Let us assume that the bandpass input signal, system, and output are $u(t)$, $h(t)$, and $r(t)$, respectively, while their corresponding Fourier transforms are indicated as $U(f)$, $H(f)$, and $R(f)$. Then, the following relations hold:

where $\otimes $ represents the convolution operator. Now, the corresponding baseband equivalents of the input signal and system are ${u}_{l}(t)$, ${h}_{l}(t)$, ${U}_{l}(f)$, ${H}_{l}(f)$. Hence, the output signal of the baseband equivalent system can be defined asIn the following, it is proven that the output of the baseband equivalent system ${r}_{l}(t)$, ${R}_{l}(f)$ and the output of the bandpass system $r(t)$, $R(f)$ have the same relations as the baseband equivalent and bandpass signals [see Eqs. (11) and (13)]. Indeed, starting from Eqs. (A6) and (A7), the following relations can be derived [24]:

## APPENDIX B: BASEBAND EQUIVALENT “SHIFTED” SYSTEM

In the following, we prove that the baseband equivalent “shifted” system represented by Eq. (17) is equivalent to the based equivalent system ${h}_{l}(t)$ in Eq. (A1), in the sense of time-domain simulations.

According to Section 3.B, the transfer function ${\widehat{H}}_{l}(f)$ and impulse response ${\widehat{h}}_{l}(t)$ of the proposed baseband equivalent state-space model from Eq. (17) can be described as

since it is obtained by shifting all the poles of the corresponding state-space model of a bandpass system by $j2\pi {f}_{c}$, considering that $\mathit{A}$ is a diagonal complex-valued matrix with all the poles as diagonal elements.By comparing the results obtained in Eqs. (B1) and (B2) to the baseband equivalent system definition given in Eqs. (A1) and (A2), one difference is clear: only the frequency response of $H(f)$ at positive frequencies is shifted by ${f}_{c}$ in the definitions Eqs. (A1) and (A2), while in Eqs. (B1) and (B2) the entire frequency response of the bandpass system considered is shifted. This difference is illustrated in Fig. 22.

Then it is proven that the relation in Eq. (A8) still holds for baseband equivalent “shifted” systems calculated by means of Eqs. (B1) and (B2). Indeed, the output signals of the bandpass system in the frequency-domain can be written as

Note that Eq. (B5) holds because ${\widehat{H}}_{l}(f-{f}_{c})$ and ${\widehat{H}}_{l}^{*}(-f-{f}_{c})$ have a non-zero frequency response at both positive and negative frequencies, while ${U}_{l}(f-{f}_{c})$ and ${U}_{l}^{*}(-f-{f}_{c})$ have a non-zero frequency response only at positive and negative frequencies, respectively.

Finally, Eq. (B3) demonstrates that the state-space representation from Eq. (17) of the baseband equivalent “shifted” system can effectively be used to replace the expensive time-domain simulations of the bandpass system.

## APPENDIX C: HAMILTONIAN MATRIX OF BANDPASS EQUIVALENT SYSTEM

Following the procedure in Ref. [17] and assuming that ${\widehat{\mathit{S}}}_{l}(s)$ is the scattering matrix of a baseband equivalent system, such a system is switching from a non-passive to a passive state (or the other way around) at the frequencies where $\mathit{I}-{\widehat{\mathit{S}}}_{l}^{H}(s){\widehat{\mathit{S}}}_{l}(s)=0$. To identify these frequencies, with input $|{\mathit{u}}_{l}|\ne 0$, we write

Let us assume that ${\widehat{\mathit{S}}}_{l}(s)$ has state-space parameters ${\widehat{\mathit{A}}}_{l}$, ${\widehat{\mathit{B}}}_{l}$, ${\widehat{\mathit{C}}}_{l}$, ${\widehat{\mathit{D}}}_{l}$; while ${\widehat{\mathit{A}}}_{l}^{H}$, ${\widehat{\mathit{C}}}_{l}^{H}$, ${\widehat{\mathit{B}}}_{l}^{H}$, ${\widehat{\mathit{D}}}_{l}^{H}$ are the state-space parameters of ${\widehat{\mathit{S}}}_{l}^{H}(s)$. Then, Eqs. (C2) and (C3) can be written in the form

Combining Eqs. (C9) and (C10) leads to

Thus, the Hamiltonian matrix of the baseband equivalent systems is

Finally, by indicating the eigenvalues of ${\widehat{\mathit{M}}}_{l}$ with the symbol ${\widehat{\lambda}}_{l}^{z}$, the following equation holds:

Now, assuming that a matrix $\mathit{M}$ exists with eigenvalues ${\lambda}^{z}$ satisfying ${\hat{\mathit{M}}}_{l}=\mathit{M}-j2\pi {f}_{c}\mathit{I}$ leads to

## Funding

Fonds Wetenschappelijk Onderzoek (FWO) (G013815N); Flanders Innovation & Entrepreneurship (VLAIO) and Luceda Photonics through the MEPIC project.

## REFERENCES

**1. **A. F. Oskooi, D. Roundy, M. Ibanescu, P. Bermel, J. Joannopoulos, and S. G. Johnson, “Meep: a flexible free-software package for electromagnetic simulations by the FDTD method,” Comput. Phys. Commun. **181**, 687–702 (2010). [CrossRef]

**2. **L. M. Zhang, S. F. Yu, M. C. Nowell, D. D. Marcenac, J. E. Carroll, and R. G. S. Plumb, “Dynamic analysis of radiation and side-mode suppression in a second-order DFB laser using time-domain large-signal traveling wave model,” IEEE J. Quantum Electron. **30**, 1389–1395 (1994). [CrossRef]

**3. **H. Bahrami, H. Sepehrian, C. S. Park, L. A. Rusch, and W. Shi, “Time-domain large-signal modeling of traveling-wave modulators on SOI,” J. Lightwave Technol. **34**, 2812–2823 (2016). [CrossRef]

**4. **S. Balac and F. Mahe, “An embedded split-step method for solving the nonlinear Schrodinger equation in optics,” J. Comput. Phys. **280**, 295–305 (2015). [CrossRef]

**5. **M. Fiers, T. V. Vaerenbergh, K. Caluwaerts, D. V. Ginste, B. Schrauwen, J. Dambre, and P. Bienstman, “Time-domain and frequency-domain modeling of nonlinear optical components at the circuit-level using a node-based approach,” J. Opt. Soc. Am. B **29**, 896–900 (2012). [CrossRef]

**6. **S. Mingaleev, A. Richter, E. Sokolov, C. Arellano, and I. Koltchanov, “Towards an automated design framework for large-scale photonic integrated circuits,” Proc. SPIE **9516**, 951602 (2015). [CrossRef]

**7. **R. Achar and M. S. Nakhla, “Simulation of high-speed interconnects,” Proc. IEEE **89**, 693–728 (2001). [CrossRef]

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

**9. **G. Antonini, “SPICE equivalent circuits of frequency-domain responses,” IEEE Trans. Electromagn. Compat. **45**, 502–512 (2003). [CrossRef]

**10. **A. Chinea, S. Grivet-Talocia, H. Hu, P. Triverio, D. Kaller, C. Siviero, and M. Kindscher, “Signal integrity verification of multichip links using passive channel macromodels,” IEEE Trans. Compon. Packag. Manuf. Technol. **1**, 920–933 (2011). [CrossRef]

**11. **D. Spina, F. Ferranti, G. Antonini, T. Dhaene, L. Knockaert, and D. Vande Ginste, “Time-domain Green’s function-based parametric sensitivity analysis of multiconductor transmission lines,” IEEE Trans. Compon. Packag. Manuf. Technol. **2**, 1510–1517 (2012). [CrossRef]

**12. **M. Sahouli and A. Dounavis, “An instrumental-variable QR decomposition vector-fitting method for modeling multiport networks characterized by noisy frequency data,” IEEE Microwave Wireless Compon. Lett. **26**, 645–647 (2016). [CrossRef]

**13. **S. Grivet-Talocia and B. Gustavsen, “Black-box macromodeling and its EMC applications,” IEEE Electromagn. Compat. Mag. **5**, 71–78 (2016). [CrossRef]

**14. **L. Chrostowski and M. Hochberg, *Silicon Photonics Design: From Devices to Systems* (Cambridge University, 2015).

**15. **D. Saraswat, R. Achar, and M. S. Nakhla, “A fast algorithm and practical considerations for passive macromodeling of measured/simulated data,” IEEE Trans. Adv. Packag. **27**, 57–70 (2004). [CrossRef]

**16. **D. Deschrijver, M. Mrozowski, T. Dhaene, and D. De Zutter, “Macromodeling of multiport systems using a fast implementation of the vector fitting method,” IEEE Microwave Wireless Compon. Lett. **18**, 383–385 (2008). [CrossRef]

**17. **B. Gustavsen and A. Semlyen, “Fast passivity assessment for S-parameter rational models via a half-size test matrix,” IEEE Trans. Microwave Theory Tech. **56**, 2701–2708 (2008). [CrossRef]

**18. **D. Deschrijver and T. Dhaene, “Fast passivity enforcement of S-parameter macromodels by pole perturbation,” IEEE Trans. Microwave Theory Tech. **57**, 620–626 (2009). [CrossRef]

**19. **R. H. Muller, “Definitions and conventions in ellipsometry,” Surf. Sci. **16**, 14–33 (1969). [CrossRef]

**20. **R. Atkinson and P. H. Lissberger, “Sign conventions in magneto-optical calculations and measurements,” Appl. Opt. **31**, 6076–6081 (1992). [CrossRef]

**21. **S. Grivet-Talocia and A. Ubolli, “On the generation of large passive macromodels for complex interconnect structures,” IEEE Trans. Adv. Packag. **29**, 39–54 (2006). [CrossRef]

**22. **D. G. Schultz and J. L. Melsa, *State Functions and Linear Control Systems* (McGraw-Hill, 1967).

**23. **J. Butcher, *Numerical Methods for Ordinary Differential Equations* (Wiley, 2003).

**24. **M. C. Jeruchim, P. Balaban, and K. S. Shanmugan, *Simulation of Communication Systems: Modeling, Methodology and Techniques* (Springer, 2006).

**25. **J. B. King and T. J. Brazil, “Time-domain simulation of passband S-parameter networks using complex baseband vector fitting,” in *Integrated Nonlinear Microwave and Millimetre-Wave Circuits Workshop (INMMiC)* (2017), pp. 1–4.

**26. **D. Youla, L. Castriota, and H. Carlin, “Bounded real scattering matrices and the foundations of linear passive network theory,” IRE Trans. Circuit Theory **6**, 102–124 (1959). [CrossRef]

**27. **P. Triverio, S. Grivet-Talocia, M. S. Nakhla, F. G. Canavero, and R. Achar, “Stability, causality, and passivity in electrical interconnect models,” IEEE Trans. Adv. Packag. **30**, 795–808 (2007). [CrossRef]

**28. **R. W. Newcomb, “On the energy in passive systems,” Proc. IEEE **53**, 1651–1652 (1965). [CrossRef]

**29. **R. Nedunuri, “On the definition of passivity,” IEEE Trans. Circuit Theory **19**, 72 (1972). [CrossRef]

**30. **M. R. Wohlers, *Lumped and Distributed Passive Networks* (Academic, 1969).

**31. **S. Boyd and L. O. Chua, “On the passivity criterion for LTI N-ports,” Int. J. Circuit Theory Appl. **10**, 323–333 (1982). [CrossRef]

**32. **B. Gustavsen, “Fast passivity enforcement for S-parameter models by perturbation of residue matrix eigenvalues,” IEEE Trans. Adv. Packag. **33**, 257–265 (2010). [CrossRef]

**33. **https://www.lucedaphotonics.com.

**34. **https://www.mathworks.com.

**35. **K. Jinguji and M. Kawachi, “Synthesis of coherent two-port lattice-form optical delay-line circuit,” J. Lightwave Technol. **13**, 73–82 (1995). [CrossRef]

**36. **S. Suzuki, Y. Inoue, and T. Kominato, “High-density integrated 1 × 16 optical FDM multi/demultiplexer,” in *IEEE Lasers and Electro-Optics Society Annual (LEOS) Meeting* (1994), Vol. 2, pp. 263–264.

**37. **Y. Xing, D. Spina, A. Li, T. Dhaene, and W. Bogaerts, “Stochastic collocation for device-level variability analysis in integrated photonics,” Photon. Res. **4**, 93–100 (2016). [CrossRef]