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.

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, a(s) and b(s) are the forward wave and backward wave, respectively, and 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 D∈Rn×n,Rk∈Cn×n,k=1,…,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 Sij(s) of the scattering matrix representation in Eq. (2) use a common denominator polynomial and pole set [p1,p2,…,pK], where such poles are either real quantities or complex conjugate pairs ^[8]. The identification of poles pk and residue matrices Rk can be performed via the VF algorithm ^[8,15–18], starting from a set of the scattering parameters under study obtained for sr=j2πfr with r=1,…,R.

However, it is important to note that the sign convention ejwt 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 ejwt 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 A∈Cm×m, B∈Rm×n, C∈Cn×m, D∈Rn×n, m=nK, and I is the identity matrix in this paper. In particular, A is a diagonal matrix with all the poles as diagonal elements, while 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 where x(t)∈Rm×1 is the state vector.

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 Δtmax that can be adopted: Δtmax must be smaller than the highest frequency component of the signals considered;

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 Ul(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.

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.

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: u(t)=A(t)cos[2πfct+ϕ(t)],(5)where A(t) is the time-varying amplitude or envelope of the modulated signal, and ϕ(t) is the time-varying phase. In electronics or radio-frequency (RF) applications, both A(t) and ϕ(t) relate to electronic signals, such as voltage, current, or electric field. In photonics, the optical carrier frequency fc is much higher than the bandwidth of A(t) and ϕ(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] ua(t)=u(t)+jH[u(t)]=A(t)ej[2πfct+ϕ(t)],(6)where H[u(t)] is the Hilbert transform of u(t). In the frequency domain, Eq. (6) becomes Ua(f)=2U(f)Step(f),(7)where Ua(f) and U(f) are the Fourier transform of ua(t) and u(t), respectively, and Step(f) is a unit step function defined by Step(f)={1,f>0,12,f=0,0,f<0.(8)

Now, the corresponding baseband equivalent signal of the bandpass signal is defined as ul(t)=ua(t)e−j2πfc=A(t)ejϕ(t),(9)Ul(f)=2U(f+fc)Step(f+fc),(10)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), H[u(t)], and ul(t) in the time and frequency domains are ^[24] u(t)=Re[ul(t)ej2πfct],(11)H[u(t)]=Im[ul(t)ej2πfct],(12)U(f)=12Ul*(−f−fc)+12Ul(f−fc),(13)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 Ua(f) has only a non-zero spectrum in the positive frequencies around the carrier frequency; shifting the spectrum of Ua(f) in the direction of the negative frequencies of fc [or equivalently in the time domain by multiplying ua(t) with e−j2πfct] leads to Ul(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 fc, satisfying fc≫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.

Baseband equivalent signals ul(t) can be easily computed with Eq. (9), where ul(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), ul(t) can be expressed with respect to its in-phase component I(t)=A(t)cosϕ(t) and quadrature component Q(t)=A(t)sinϕ(t) as ul(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 Ul(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 a(t), b(t), and x(t) in the system of ODE from Eq. (4) in the form from Eq. (11), which gives {Redxl(t)ej2πfctdt=ARe[xl(t)ej2πfct]+BRe[al(t)ej2πfct],Re[bl(t)ej2πfct]=CRe[xl(t)ej2πfct]+DRe[al(t)ej2πfct],(14)where al(t), bl(t), and xl(t) are the baseband equivalents of a(t), b(t), and x(t), respectively. Next, using the Hilbert transform and the relation from Eq. (12) to represent a(t), b(t), and x(t) in Eq. (4) leads to {Imdxl(t)ej2πfctdt=AIm[xl(t)ej2πfct]+BIm[al(t)ej2πfct],Im[bl(t)ej2πfct]=CIm[xl(t)ej2πfct]+DIm[al(t)ej2πfct].(15)

Equations (14) and (15) allow us to write {dxl(t)ej2πfctdt=Axl(t)ej2πfct+Bal(t)ej2πfct,bl(t)ej2πfct=Cxl(t)ej2πfct+Dal(t)ej2πfct.(16)

After simple mathematical manipulations, Eq. (16) can be written as {dxl(t)dt=(A−j2πfcI)xl(t)+Bal(t),bl(t)=Cxl(t)+Dal(t),(17)which represents a realization of the baseband equivalent system by means of the state-space matrices A−j2πfcI, B, C, and D; in this contribution we define it as baseband state-space model. It is evident that such a model can be obtained by directly shifting all the poles of the corresponding state-space model from Eq. (4) of bandpass system by j2πfc, considering that A is a diagonal complex-valued matrix with all the poles as diagonal elements, as mentioned in Section 2.

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πfc, 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.

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 S(s) is analytic in Re(s)>0;

The superscriptH stands for the transpose conjugate operator. The first condition is related to causality and stability; the second one is basically a bound for 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.

According to Refs. ^[26,28,29], an n-port electronic system is passive if, for any τ>−∞ and v(t)∈L2n (L2n denotes the space of all vectors whose n components are functions of a real variable t and square integrable over −∞<t<∞), it holds Re∫−∞τvH(t)i(t)dt≥0,(18)where v(t) and 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 a(t) and backward b(t) waves, the passivity definition becomes ^[26,30,31] ∫−∞τaH(t)a(t)−bH(t)b(t)dt≥0,(19)which is more convenient for describing photonic systems.

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 S^l(s) of the baseband equivalent systems: (1)

S^l(s) is analytic in Re(s)>0;

Note that real-valued systems need the extra condition S(s*)=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.

Passivity conditions require both scattering parameters S(s) and their baseband equivalent S^l(s) to be bounded by unity, which implies that all singular values σ of S^l(s) are smaller than unity at all frequencies σi(f)<1,i=1,…,n.(20)

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 M ^[17], defined as M=[A−BL−1DTC−BL−1BTCTQ−1C−AT+CTDL−1BT],(21)where A, B, C, D are the state-space matrices in Eq. (4), while L=DTD−I and Q=DDT−I.

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 M^l for baseband equivalent systems S^l(s) can be derived by following the procedure in Ref. ^[17], leading to Ml=[A^l−B^lLl−1D^lHCl−B^lLl−1B^lHC^lHQl−1C^l−A^lH+C^lHD^lLl−1B^lH],(22)where A^l, B^l, C^l, D^l are the complex-valued baseband equivalent state-space matrices, while L^l=D^lHD^l−I and Q^l=D^lD^lH−I. The derivation of M^l is shown in Appendix C.

One can observe that the only difference between M and M^l is the use of the transpose conjugate operator for the state-space matrices in Ml, while only the transpose operator is required in M. Indeed, state-space models of general electronic or photonic systems satisfy the conjugacy property S*(s)=S(s*); the corresponding scattering parameters do not change if the state-space matrices A, B, C, 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 A^l, B^l, C^l, D^l in Eq. (22) with A^l=A−j2πfcI,B^l=B,C^l=C,D^l=Dgives M^l=M−j2πfcI,(23)where M is the Hamiltonian matrix of the corresponding bandpass system. Then it is easy to derive (see Appendix C) λ^lz=λz−j2πfc,forz=1,…,Z,(24)where Z is the total number of eigenvalues, while λz and λ^lz are the eigenvalues of M and Ml, respectively.

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.

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);

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 fc=193THz 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,∞]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.

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.

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 π/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 fc=193.72THz, 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 [fc−Δ,fc+Δ], where the choice Δ=150GHz 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 −60dB 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.

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 R1 and R2. 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 [fc−Δ,fc+Δ], with fc=195.75THz and Δ=450GHz. The maximum absolute error of both models is less than −65dB. 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.6ps, 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.

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 fc=195.11THz (λ=1.5365μ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 (λ=1.5370μm) or 194.98 THz (λ=1.5375μm), as shown in Fig. 20. Model LB is built for each one of these three situations by adopting a pole shift of fc=195.11THz, 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 −60dB.

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.

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.