Main parameters of the SOA device used in the simulation
Semiconductor optical amplifiers (SOAs) have been extensively used in a wealth of telecom and datacom applications as a powerful building block that features large optical gain, all-optical gating function, fast response, and ease of integration with other functional semiconductor devices. As fabrication technologies are steadily maturing toward enhanced yield, SOAs are foreseen to play a pivotal role in complex photonics integrated circuits (PICs) of the near future. From a design standpoint, accurate numerical modeling of SOA devices is required toward optimizing PICs response from a system perspective, while enhanced circuit complexity calls for efficient solvers. In this book chapter, we present established experimentally validated SOA numerical modeling techniques and a gain parameterization procedure applicable to a wide range of SOA devices. Moreover, we describe multigrid concepts and implicit schemes that have been only recently presented to SOA modeling, enabling adaptive time stepping at the SOA output, with dense sampling at transient phenomena during the gain recovery and scarce sampling during the steady-state response. Overall, a holistic simulation methodology approach along with recent research trends are described, aiming to form the basis of further developments in SOA modeling.
- Semiconductor optical amplifier
- numerical modeling
- transfer matrix method
- multigrid techniques
The explosive growth of information traffic and the concomitant increase of bandwidth hungry applications have contributed to a growing demand for communications networks offering greater bandwidth and flexibility at lower cost [1–3]. This has led to a series of technological advances in high-speed backbone networks and telecommunications, where integrated optical systems delivering low-cost, low-power, high-bandwidth transmission, and high-speed switching elements are increasingly required . Toward meaningful and functional subsystems, high-speed all-optical elements such as all-optical signal wavelength conversion , on/off keying modulation [6,7], header recognition , optical buffering , and signal regeneration  have been developed, targeting to perform basic functionalities in the optical domain. Meanwhile, photonic integration technologies have steadily matured over the last decade, achieving improved yield fabrication [11,12], allowing the fabrication of complex single-chip photonic modules with advanced functionality. This has resulted in impressive demonstrations of functional complex circuits and switching architectures [12–17], making firm steps toward medium-scale (MS) photonic integrated chips (PICs).
On the way to develop complex devices with enhanced processing capabilities, semiconductor optical amplifiers (SOAs) have been a powerful building block that provides high-speed all-optical switching operations  by featuring large optical gain, optical gating function, controllable performances by injection current, compactness, fast response, and ease of integration with other functional semiconductor devices . SOAs are utilized in telecom [19,20] and datacom [21,22] PICs either as single SOA travelling waveguides, supporting cross-gain modulation (XGM) phenomena, or arranged in SOA–Mach–Zehnder interferometer (SOA-MZI) configurations, supporting cross-phase modulation (XPM) phenomena. The reason lies in the maturity of the SOA technology to a point where commercial devices are readily available either as bulk single chip elements [23–25], in arrays of certain pitch [26,27], or packaged in discrete fiber pigtailed components  for use in optical communication systems. Moreover, with recent research achievements on mid-board flip-chip bonding of an SOA array to a silicon-on-insulator (SOI) platform  or the development of temperature stable SOA , SOAs are expected to play a pivotal role in future PICs.
In this regime, SOA numerical modeling has attracted a lot of research attention during the last decades, as mathematical models are required to aid in the design of SOAs and to predict their operational characteristics [29–36]. During the development process and prior to fabrication, SOA models that take into account the semiconductor properties of the active material and cross-sectional dimensions are employed when optimizing the amplification gain and output power over a wideband steady state and under various external current injections. In addition, during the system-level performance study of SOA-based circuits, time domain simulations are required for an accurate evaluation of the circuit response with multiple signals propagating bidirectionally along the waveguide. With growing complexity of circuits, system level designs call for experimentally validated yet efficient solvers that support multiple signals at different wavelengths propagating in multiple directions.
This book chapter aims to provide a holistic simulation methodology approach for efficient and validated numerical modeling of SOAs and SOA-based circuits in the time, which may form the basis of further developments in SOA modeling. This chapter is organized as follows: a short review of established SOA modeling approaches in the literature will be presented in Section 2, along with an overview of some recent research efforts. In Section 3, a description of an experimentally validated time domain SOA model relying on the TMM will be presented [37–41], relying on explicit modeling techniques, followed by a gain parameterization procedure for tailoring the TMM model to characterization measurements of other SOA devices in Section 4. Having developed an experimentally validated SOA model, we numerically investigate the impact of the external current injection, the SOA length, and the light propagation direction in the gain recovery using pump–probe techniques to qualitatively study the shortening of the gain recovery time toward achieving high bandwidth switching operations. Section 5 presents a newly introduced efficient solver in the time domain, relying on implicit schemes and multigrid concepts. Implicit schemes allow developing an adaptive step size-controlled solver for the WDM SOA response [42–44], which extends the validated TMM model by lifting the limitations of spatiotemporal grids and releasing adaptive sampling at the SOA output. Toward application in realistic system-level scenarios, the last section presents experimentally validated numerical results for two predominant SOA-based all-optical signal processing circuits, such as a coupled SOA XGM flip-flop arrangement and an XOR gate based on SOA-MZI XPM configuration. The former acts as memory for sequential logic processing circuits, while the latter forms a basic building block of combinational logic signal processing circuits. Conclusions are addressed in the final section.
2. Established SOA numerical models
SOA numerical modeling has progressed on multiple fonts during the last decades, and various simulations approaches have been developed [29–36]. SOA models can be roughly divided in two categories: (i) the material models that target optimization of the emission properties, before and during the development-fabrication process of an SOA device and (ii) the circuit models, which account for the carrier density, the interband, and the intraband phenomena for simulations of light–matter interaction in realistic system-level application scenarios.
The material models have long studied the gain coefficients, emission characteristics, and spectral properties of the active material. Optimization parameters include the peak gain; the molecular fractions of the materials in the semiconductor compounds; the cross-sectional dimensions, e.g., ridge width, core height, etc.; and the structural properties, e.g., number of quantum wells (QW), quantum dots (QD), and so on. Gillner et al.  investigated long-wavelength semiconductor laser amplifiers by means of an experimentally validated SOA model and considering different structural parameters such as thickness of the active layer and amplifier length. The model that takes into account Auger recombination, thermal effects, and spontaneous emission was developed in order to optimize the spectral gain properties of SOA lasers, such as peak gain wavelength shift and width of gain curve. Interestingly, it was shown that there exists an optimum active layer thickness with respect to current density for a certain gain, while increased SOA length allows higher gains with reduced wavelength variation of the peak gain. Similarly, a comprehensive model was presented by Minch et al.  for the calculation of the band edge profile of both the In1–xGaxAsyP1–y and In1–x–yGaxAlyAs quantum-well systems with an arbitrary composition, as typical semiconductor compounds found in SOAs. This model provided in-depth knowledge with mathematical curve fits between the measured net modal gains for both material systems and calculations from the realistic band structure, including valence band mixing effects, allowing to extract the relationship between total current density and carrier density. More recently, research has focused on studying the properties of newly introduced active compound materials, such as the GaInNAs. The analysis of the broadband gain of a GaInNAs single quantum-well (QW) SOA that takes into account the tunability of the gain and incorporates quantum dot (QD) fluctuations due to compositional fluctuations of N within the QW is presented by Xiao et al. . Finally, III/V-on-SOI bonding processes have also been into design consideration recently, with a comprehensive model presented by Cheung et al.  providing valuable insight to the design requirements of the optimization process of the structural and material parameters, e.g., the width, composition, and number of quantum wells of a compressively strained In(1–x–y)Ga(x)Al(y)As quantum-well active region for emission at approximately 1550 nm.
This chapter focuses more on the 2nd category, toward presenting a holistic simulation approach of time domain numerical modeling of circuit-level simulations. An advanced time domain dynamic numeric model that accounts for the ultrafast gain dynamics, including the intraband phenomenon, the gain saturation, and the gain spectral profile of an SOA, has been presented by Toptchiyski et al. , where the model is employed in a system-level simulation to investigate the gain dynamics of the light-propagation direction during pump probe measurement and application in an Sagnac interferometer switch semiconductor laser amplifier in a loop mirror (SLALOM). By exploiting the carrier density fluctuation due to interband phenomena within the SOA, many applications have been demonstrated already. Accurate designs and numerical analysis of various SOA-based devices have resulted in successful demonstrations of optical flip-flops and random access memory (RAM) cells, or SOA–Mach–Zehnder interferometers (MZIs) and cutting edge optical routing devices, profiting from cross-gain modulation (XGM) or cross-phase modulation (XPM) phenomena . Meanwhile, circuit-level modeling approaches have been developed for experimental verification and accuracy.
Traditional time domain SOA models in principle treat the gain dynamics based on the longitudinal carrier distribution along the device, relying on the free carrier rate equation for the electron–photon interaction. To achieve this, a fine longitudinal SOA discretization is usually adopted. Dividing the SOA into cascaded elementary sectors of equally small propagation length and time stepping allows applying wideband steady-state material gain coefficients  across the emission spectrum, using a one-by-one space and time representation of the carrier densities and photon fluxes. This discretization derives from an explicit scheme for solving the associated differential equations, resulting in an equidistant grid with constant time stepping. This technique is the standard approach used in circuit-level designs and has proven capable of delivering quantitative matching with experimental measurements [30–35]. Following this concept, an experimentally validated time domain SOA model relying on the TMM analysis technique was developed by Vagionas et al.  with the simulation results coming in close agreement with the characterization measurements on a commercial SOA device in terms of the emission spectrum, gain profile, and recovery time. This model was later on extended with the plasma effect  to account for the effective refractive index change and hence the phase shift term in order to provide for a comparative study of the performance characteristics of XGM- and XPM-based circuit topologies, indicating that the XPM phenomena are more resilient to high-speed operation.
As adopting explicit schemes for solving the system of coupled differential equations can become sometimes quite large even for a single SOA device, efficiency also needs to be addressed. Especially as integration technology and circuits evolve into complex networks, simulating multiple SOA devices simultaneously in the time domain in complex circuit configurations or with multi-wavelength signals travelling in various directions further exacerbates the numerical requirements. Simplified analytical gain models provide only qualitative conceptual results by suggesting exploring computational efficiency at the cost of accuracy . In this regime, a few SOA models have been recently presented, demonstrating techniques to speed up the simulation time at the cost of some little accuracy. For instance, the use of one auxiliary signal with an effect equivalent to the amplified spontaneous emissions (ASE) was presented by Vujicic et al.  in order to reduce the wavelength channels under study. Almost similarly, a single-state variable depending on the available carrier densities within the SOA, termed as a “reservoir” of excited carriers, was presented by Mathlouthi et al. , resembling the reservoir of excited erbium ions in an erbium-doped fiber amplifier (EDFA). Despite the state variable, the model still relies on an equidistant spatiotemporal grid with constant time sampling at the output of the SOA, even at relatively low bitrates of 1 Gb/s.
3. Accurate Transfer Matrix Method (TMM) SOA model
The TMM numerical analysis technique was initially presented by Davis et al.  for DFB lasers and was initially presented for the extraction of the complete transfer matrix (TM) of large structures with rather simple configuration. The TMM divides the SOA longitudinally into m small sectors of length Δz = L/m, where L is the SOA length along the propagation direction. Considering a short sector length Δz and a short time interval Δt, the light requires only a small time interval Δt = ug∙Δz, where ug is the group velocity, to propagate through each sector. During this small time interval, the structural and material parameters of each sector, such as the carrier density and the material gain, can be considered uniform and constant, while transverse variations are still not allowed. This allows considering a wideband steady-state gain model.
The light propagation through a single sector can then be described by the simple TM, which relates the incoming and outgoing amplitudes of an elementary waveguide sector. For the current analysis and without loss of generality, we assume a signal Es propagating to the right direction and a second signal counterpropagating toward the left side, Ec. By applying the TM of sector i at time t for both the externally injected lights and the spontaneously emitted photons over a wide optical spectrum of 1500–1600 nm, we obtain the following equation:
where the propagation constant γ(λi) = geff/2 + j∙neffω/c describes the amplification (real part) and the phase shift (imaginary part) of the wave amplitude, after propagating through an elementary sector waveguide, as illustrated in Figure 1. Each updated wave amplitude that emerges at the output facet of a single sector can in turn be regarded as incident to the next sector, and the process is repeated until both lights exit the SOA. At each iteration, the waveguide propagation constant is updated and a new operant TM is estimated. This operation implies that a TM of a sector is time dependent since the matrix elements are subject to the net gain coefficient geff and phase shift changes neff. Figure 1(i) depicts the light propagation across the SOA lengths according to the TMM, while Figure 1(ii) depicts the ith sector and the operation of the TM on the incoming and outgoing wave amplitudes.
For the estimation of the real part of the constant, the Connelly gain coefficient has been employed, which was initially presented by Connelly  for the estimation of a wideband steady-state gain approximation of a bulk InP–InGaAsP homogeneous buried ridge stripe SOA, supporting large operating regimes. The gain coefficient is given by
where c is the speed of light in vacuum, na is the refractive index in the active region, τ is the radiative recombination lifetime, v is the optical frequency, me and mhh are the conduction band electron and valence band heavy hole effective masses, respectively, is the Planck’s constant divided by 2π, Eg is the bandgap energy and fc(v), and fv(v) is the Fermi–Dirac distributions in the conduction and valence band, respectively. Subsequently, we incorporate a loss coefficient along the SOA, toward a net gain coefficient given by 
where αα describes the losses in the active region, αc the losses in the cladding layer, and αs the scattering losses m and by αWG; we denote all cumulative losses. Carrier-dependent losses have also been presented in the literature, e.g., as given by α(N)-Κ0+ΓΚ1Ν, with K0 accounting for the overall material and waveguide losses of Eq. (3) and the second term accounting for the intervalence band absorption, which is current dependent. However, the constant term has also been shown to provide very good experimental matching. After the estimation of the updated amplitude, amplified by the net gain coefficient through propagation by one sector, the power emerging at the output facet is given by 
where hv is the photon energy, w is the width, and d is the thickness of the active region of the SOA, while the photon density of each stream is approximated by 
Finally, more rigorous electron statistics and fast nonlinear phenomena of the spectral hole burning and the intraband phenomena are only taken into consideration through the use of the gain compression factor ε, as governed by 
where J is the injection carrier density; c1, c2, and c3 are the current leakage, radiative recombination, and Auger recombination rates; Γ is the confinement factor; and Si is the photon density at the ith sector. The index u = 1, 2, 3 refers to the ASE or the two input streams.
Having updated the carrier density, we can proceed to the estimation of the imaginary part of the propagation constant, which accounts for the phase shifting term of the amplitude. The estimation of the effective refractive index neff in the imaginary part relies on the plasma effect, according to which the change of the refractive index Δnpl is linearly dependent on the change of the free carrier density in the active region ΔN :
where e is the electronic charge, ω is the angular frequency, ε0 is the permittivity of free space, neq is the effective index of the waveguide, and me and mhh are the effective mass of the electron and the heavy hole. An accurate estimation of the real and imaginary part of the propagation constant γ(λi)=geff/2+j∙neffω/c of the ith sector allows describing the forward and backward (counter)propagating amplitudes through a simple time dependent TM given in Eq. (1).
4. Gain parameterization procedure
In this section, we report on a gain parameterization procedure and a methodology toward tailoring the TMM model and producing simulation results close to the experimental measurements obtained from a commercial SOA device. The methodology is in principle generic and compatible with the other SOA devices. The target is to identify a set of material parameters and cross-sectional characteristics by incorporating feedback from experimental measurements such as the bandgap shrinkage to fit the experimental gain spectrum in a spectrum analyzer, the recombination rates to match the gain profile of gain-vs-incoming power measurements, and the gain compression factor for the recovery time at pump probe.
The procedure includes a number of steps clustered into two sections: first, the definition of the material and waveguide properties and, second, the feedback from the experimental characterization of the SOA device under study.
The material and waveguide properties include the following:
To define the length of the active region of the fabricated SOA waveguide, according to foundry.
To estimate of the initial bandgap energy Eg0 of the active region, assuming no external current injection
To define the loss coefficient of the waveguide
To incorporate the waveguide cross-sectional dimensions, e.g., core height and width
To define the external operating current I and the current density injected to the active region.
The cross-sectional dimensions and the length of the SOA waveguide or the active material can be provided by the fabrication foundry. For the definition of the material properties, extensive literature [45–48] has been developed on III/V compounds, especially for common compounds such as the InGaAsP-based devices. The initial bandgap energy Eg0 when no carriers are injected in the active region can be provided by the quadratic approximation , with e being the electronic charge and y the molar fraction of the Arsenide in the active region, and a, b, and c being the coefficient of the quadratic approximation. Similarly, the comprehensive model presented by Cheung et al.  and dealing with the optimization of the structural and material parameters of In(1–x–y)Ga(x)Al(y)As quantum-well (QW) SOAs toward bonding on an SOI platform provide an in-depth insight into bandgap energy.
In addition to the definition of the material properties, feedback from a series of experimental characterization measurements on real SOA devices can be incorporated in the model through the following three steps:
To define the bandgap shrinkage coefficient to fit the experimental gain spectrum of the SOA. The bandgap shrinkage coefficient Kg affects the blue shifting of the gain peak at the ASE output spectrum of the SOA, which occurs when increasing the injection current. An increase in the carrier density is translated into shrinkage of the material bandgap by a term . In order to extract the parameter for the Kg coefficient, static measurements of the ASE spectrum at the output of the SOA under different external current injection and no injected optical power are suggested.
To define the recombination rates for the gain curve versus the incoming optical power. This includes the c1, c2, and c3 parameters that are associated with the linear recombination coefficient at defects (current leakage), the spontaneous recombination, and the Auger recombination rate, respectively. In practice, these three mechanisms exhibit complicated dependences that together account for the overall recombination rate given by . This step requires static power gain measurements using a continuous wavelength light source close to the gain peak wavelength and under varying average optical power in order to estimate net gain of the SOA, and thus defining the small signal gain and the saturation point of the SOA.
To define the gain compression factor to fit the recovery time. The gain dynamics associated with the intraband phenomena are associated with a fast recovery time immediately after the injection of a short pulse, with the gain compression increased under increasing pulse energy and bias current. Although the TMM model focuses on the interband phenomena, the ultrafast intraband phenomena have been considered through the use of a power-dependent gain compression coefficient , which allows tuning the fast recovery time of the intraband phenomena. In this strep, typical pump probe measurements are required with a weak Continuous Wavelength signal with average optical power in the small signal gain fed to an SOA along with a high-power ultrashort return-to-zero pulse with a full width half maximum duration of 5 ps, which sweeps all the available carriers.
Considering the numerous interrelated phenomena that affect the performance of an SOA, determining the optimum parameters that match all system-level parameters, such as the gain–power curve, the spectrum, and the recovery time, implies a certain difficulty. Defining a certain methodology for parameter extraction ensures accuracy of the model for direct applicability in realistic use case scenarios and demonstrations. After developing a certain gain parameterization procedure, the use of multiobjective genetic algorithms supporting an automated iterative procedure can also been adopted  toward a best-fit criterion.
4.1. Quantitative experimentally verified numerical measurements
Following an experimental characterization on a commercial SOA device , the previous gain parameterization procedure was adopted in order to extract the optimum simulation parameters and tailor the TMM model. Figure 2(i) illustrates the experimental setup used for the characterization of the SOA device. Two tunable laser diodes (TLD) operating at wavelengths λ1 and λ2 were employed characterization of the SOA. The CW output at λ1 of the first TLD was amplified through an EDFA and filtered through a band-pass filter (BPF) before being driven to the input of the SOA. The CW output at λ2 of the second TLD was modulated by means of an electroabsorption modulator (EAM), whose operation was driven by an electrical clock signal. In this way, short optical return-to-zero (RZ) pulses were produced amplified through an EDFA and filtered through a BPF at λ2 wavelength. The produced RZ pulse stream was coupled with the λ1 CW and also fed to the SOA. Polarization controllers (PCs) were employed at various stages of the experimental setup to adjust the polarization of the optical signals as discrete off-the-shelf fiber pigtailed components were employed. Using the presented experimental setup, all three measurements of steps 6, 7, and 8 were conducted using variable optical attenuators (VOAs) at the fiber-pigtailed branches before the SOA. VOAs allowed controlling the average power of the optical signals properly. In step 6, no optical input fed to the SOA, and both VOAs provided maximum attenuation, while the SOA output was recorded at an optical spectrum analyzer (OSA). By measuring a peak gain around 1560 nm, at step 7, the CW signal of TLD1 was tuned at λ1 = 1555 nm, and the output of the SOA was measured by a power meter, after filtering through a BPF.
The TMM model was in turn tailored to match the experimental measurements. Figure 2(ii) illustrates the normalized simulated optical field distribution along the SOA using the TMM, with the vertical z-axis being the optical power and the horizontal x-axis being the wavelength range between 1520 and 1580 nm. The plot illustrates how an input simulation signal at 1555 nm wavelength of the SOA is propagated together with the generated ASE noise sector by sector, as marked in the y-horizontal axis, until it emerges at the final sector at the SOA output. Moreover, Figure 2(iii)–(v) demonstrates the simulation-experimental matching: The gain versus input power in Figure 2(iii) for a single λ1-CW signal reveals accurate matching both in the small signal region as well as in the saturate gain regime, with less than 1 dB error for an input power between –30 and 0 dBm. Additionally, the simulated output spectrum fits well the experimentally recorder spectrum with less than 2 dB errors, while the matching in the time domain of the 10%–90% gain–recovery time of the SOA reveals accurate prediction of the response of the SOA in the ps time scale. The simulation parameters extracted by this procedure are presented in Table 1.
|m||Number of longitudinal||20|
|Δt||Temporal sector interval||1 [ps]|
|L||SOA length||1600 [μm]|
|w||Active layer width||1.2 [μm]|
|d||Active layer thickness||0.1 [μm]|
|u g||Group velocity||8.5·108 [m/s]|
|c 1||Linear recombination rate||1·107 [m/s]|
|c 2||Bimolecular recombination rate||50·10–17 [m3/s]|
|c 3||Auger recombination rate||80·10–41 [m6/s]|
|a wg||Waveguide losses||5800 [1/m]|
|R||Reflectivity at facets||0|
|I||External current injection||300 [mA]|
|n ao||Reflective index in the active region||3.22|
|ncl||Refractive index in the cladding||3.1|
|Eg||Bandgap energy||0.7773 [eV]|
4.2. Qualitative results on SOA gain recovery time
To accommodate the increasing demand of data transfer and high-speed optical telecommunication networks with terabit transmission capabilities and high bandwidth switching functionalities, there has been a growing interest in increasing the recovery time of the SOA. High-speed nonlinear SOAs are used to perform either XGM or XPM modulation between two input signals, i.e., a weak CW probe signal and a high-speed data signal with a short pulse width, acting as pump. Increasing the external current injection or the elongating the SOA active region have been shown to shorten the recovery time and thus support higher speed switching operations. After developing an experimentally verified time domain model, we numerically investigate the shortening of the recovery time during pump probe measurements under increasing current or increasing SOA length.
Figure 3(i) illustrates the simulated gain recovery and the normalized output power of the SOA for various external currents. The external currents were tuned from 185 to 350 mA, with the 300 mA being the nominal current operation. The plots reveal a shortening at higher currents, stemming from the increased current density injected to the SOA. The measured 1/e gain recovery times were plotted versus the supplied current in the inset, demonstrating a shortening from 45 ps under 185 mA down to 20 ps under 350 mA.
Similarly, Figure 3(ii) illustrates the simulation results for the normalized SOA power output and the gain recovery after pump probe measurements, while the external current injection was maintained constant at nominal values of 300 mA. The plots reveal a shortening of the recovery time when elongating the SOA from 0.5 to 1.7 mm. The obtained 1/e recovery time was measured and plotted in the inset, indicating a shortening from 63 ps for a 0.5-mm-long SOA down to 20 ps for a 1.7-mm-long SOA.
The effect of the propagation direction of the pump signal was also investigated using the TMM SOA model. Figure 4 depicts the simulated time traces of the normalized output power and the phase shift of the probe light, obtained after pump probe measurements with the pump control pulse being fed in a co- or counterpropagating direction. The CW probe signal featured an average power of 50 μW, which lies in the small signal regime, while the peak of the control pump pulse was tuned, so as to induce a π phase shift at the CW signal for either propagation direction.
The normalized traces are plotted of the XPM operation imprinted on the probe signal are illustrated in Figure 4(i) and the XGM operation in Figure 4(ii). It is obvious that in both cases, a π phase shift was introduced due to XPM operation. However, the probe signal is suppressed more in the case of the copropagating pump pulse, reaching normalized power values down to 0.1, while for the counterpropagating pump pulse the gain suppression is less deep. This owes to the increase pump–probe light interaction along the copropagating direction through the carrier rate equation. In addition, to the recovery time in case of a counterpropagating control pulse is shorter (after the suppression of the gain a local minima), however, the trade-off is a slower response/fall time observed with a delay in the deep (minima of the output power) when the control pulse is injected in the counterpropagating direction. This owes again to the travelling time and the time required for the counterpropagating pump pulse to absorb the available carriers. The findings for XGM effects and XPM effects for π-phase shift of a Gaussian control pulse are shown in Figure 4 and are in full compliance with the ones found in the literature .
5. Efficient multigrid SOA model with adaptive time stepping (TS)
Numerical SOA modeling has so far relied on explicit schemes for solving the associated system of coupled ordinary differential equations (ODE), comprising the spatial discretized carrier density rate equation given in Eq. (6) combined with the material gain coefficient in Eq. (2). In the literature, ODEs are classified in two categories: the stiff and nonstiff ODEs. The latter can be solved efficiently by explicit time stepping schemes. The former requires many steps with explicit schemes as warranted by the smoothness of the solution, as has been the case for the traditional SOA modeling and is schematically illustrated in Figure 2(ii) for the TMM model, where each node of the plot represents a small time step. An interesting alternative option would be to deploy implicit schemes, which could alleviate the problem of the many unneeded, as far as accuracy is concerned, time steps at the cost of having to determine the Jacobian for the set of ODEs and inverting a matrix, at each time step. Implicit schemes have shown to be more efficient for many problems , with multigrid methods known to be among the most efficient solvers for many partial differential equations (PDEs) .
Vagionas and Bos  presented a novel multigrid solver for the dynamic response of the SOA, relying on the wideband steady-state gain coefficient of Eq. (2). Introducing multigrid techniques in SOA modeling enabled extending the accurate time domain modeling of the TMM model, allowing for the development of an efficient solution supporting implicit time discretization schemes. Implicit schemes in turn enable accuracy—instead of stability—restricted time discretization of the signals. This implies a different discretization scheme, where sampling of the signals is optimized for certain accuracy in the solver instead of certain time step restriction. This allows lifting the limitations of an equidistant spatiotemporal grid for the representation of the incoming signals adopted by traditional explicit SOA models, releasing an adaptive step size-controlled solver for the dynamic SOA response with dense time sampling under a rapidly varying SOA signal output and scarce time sampling when negligible changes are observed. Adaptive times stepping adds one more degree of freedom to computational efficiency of accurate SOA modeling, which is of crucial importance when evaluating large input patterns for statistical signal analysis independent of the bitrate or large circuit networks with multiple SOA-based components.
Multigrid methods employ a series of coarser grids to obtain grid independent convergence rates. Drawing from the finest grid of traditional longitudinal division of the SOA into cascaded sectors, coarse and coarser grids with less number of SOA sectors are adopted in order to represent the spatial free carrier density distribution by less and less grid points (carrier density samples). This is schematically illustrated in Figure 5(i), where 4 grid levels have been employed. The 4th level is the finest level, including an SOA longitudinal discretization into 16 sectors, equal to the discretization employed in the TMM model. However, by applying multigrid concepts, the carrier densities of two neighboring sectors can be represented by a single sector in the coarser 3rd-level grid, resulting in half grid points. Equivalently, the 3rd-level can be again restricted to coarser grids, with each transition halving the number of grid points.
In the following the notation, h =Δz and H=2Δz are used for the representation of the SOA sector length in the fine and coarse grid, respectively. In order to focus on the multigrid aspects, the rate equation of Eq. (6) is written in an operator form as follows:
where the operator L, although simple in appearance, is the quite complex, as it involves solving the forward and backward propagation equations for both the signals and noise photon fluxes. By applying the implicit midpoint rule, we obtain the following equation:
where . The transition from a fine toward coarser grid can be achieved through an averaging restrict operator , where a coarse grid point can be considered as an average of two neighboring fine points. Equivalently, the operator , refines the coarse grid functions to a finer grid. In order to generate an intermediate fine point between two coarse grid points, an inverse operator has to be considered compared to the restrict operation, which in our case is the linear approximation. In this way, a multigrid correction V-cycle has been developed, described by the following five steps:
The developed multigrid V-cycle, including the above five steps, is schematically represented in Figure 5(ii), showing four grids of different granularity. Each coarse grid comprises half the grid points compared with the higher/finer grid level, while transitions rely on the restrict and refine operators. The smooth operation exist of second-order distributive Jacobi relaxations has been considered, that is, for each, , an update is calculated as follows:
The solve operation can be implemented recursively, while the solution of the coarse grid problem is used to correct the fine grid approximation. High-frequency components owing to the interpolation of the correction cycle can be removed by the final smooth operations. The presented cycle features a grid-size-independent rate of reducing in the error by . When using the solution of the previous time step as initial approximation to the solution at the current time step, k, O(log(1/h)) cycles are needed to solve the problem to the level of the truncation error. This can be reduced to O(1) by using the coarser levels to generate an initial approximation accurate to within the level of truncation of grid, H .
The proposed solver employs the rate equation (Eq. 6) and the multigrid techniques to solve the carrier density distribution along the SOA in coarser and coarser grids, while the propagation and amplification of the signal is still based on the Connelly material gain coefficient employed in the TMM model . Adopting the previous gain coefficient ensures equivalent steady-state results, such as the optical spectrum and the net gain, and tailoring of the SOA parameters with experimental measurements. On the other hand, incorporating an implicit time discretization scheme and adaptive time sampling suggests that computational efficiency is exploited based on the adaptive time sampling in order to benefit from long bit patterns or small pattern changes at the input bit-streams.
The increase or decrease of the time step is controlled by the implicit midpoint rule with a doubling scheme through the accuracy tolerance parameter ε, as depicted in Figure 6 for various accuracy tolerances of the multigrid solver and single block pulses at different bitrates. The simulation results depict the pulsed NRZ transmission through a single SOA with single pulses of –25 dBm peak power centered at 1550 nm for a bitrate of 1 Gb/s in Figure 6(i) and 5 Gb/s in Figure 6(ii). The accuracy tolerance ranges between 1e–2 for the plot at the top right corner of the graphs and 1e–5 for the bottom left plot, highlighting that a dense time sampling is adopted immediately at the rise and fall time of the pulse, where transient phenomena are observed. On the contrary, the trailing part of a bit pulse before the bit transitions of the output pattern, a steady state is obtained in carrier density distribution resulting in steady-state SOA gain dynamics. Thus, negligible changes are observed at the output of the SOA, revealing than the time step can be adapted to a larger value for enhanced computational efficiency.
The multigrid solver was also employed toward simulating the XGM operation between two signal in pump probe measurements. A pump control bit stream is wavelength converted on a CW probe signal at 10 Gb/s bitrate, with the CW featuring an average power of –25 dBm at 1550 nm and the NRZ bit stream exhibiting a peak power of– 25 dBm at 1555 nm. A custom control bit pattern of 001100 at 10 Gb/s was employed, as this pattern changes at odd bits (1st, 3rd, and 5th bit) with fast changes in the carrier dynamics necessitating dense time sampling. On the contrary, during the even bits, when no changes are observed at the logical values of the input data streams, the response of the SOA is characterized by a constant steady-state response and thus requires only a few time samples. This is characteristically illustrated in Figure 7(i), where the even bits necessitate only very few samples, as highlighted with a yellow marker at the cost of a user-defined accuracy error. Adaptive time stepping can be important for statistical performance evaluation of SOA-based circuits, especially when long patterns or complex structures are employed.
In order to evaluate the performance and computational efficiency of the solver, the operation of wavelength conversion of a control pulsed signal to a CW input signal has been simulated at different bitrates and accuracy levels. The evaluation is based on the number of required samples that describe the output of the SOA, considering an NRZ control signal following a 27-1 pseudo random bit sequence (PRBS) at bitrates of 1 Gb/s, 5 Gb/s, and 10 Gb/s. The results have been summarized in Figure 7(ii), showing the average number of samples per bit required for the overall simulated pattern versus the bitrate. Three different accuracy errors have been considered, namely, 1e–3, 1e–4, and 1e–5. Especially in the case of relatively big tolerance of ε = 1e–3 in the first column, bits are resolved by 4 and 6 samples for the bitrates of 10 Gb/s and 5 Gb/s, respectively. Furthermore, for bitrates of 1 Gb/s, an average of 10 samples per bit is required, suggesting that the 4 samples describe efficiently the first 100 ps (for 10 Gb/s operation), additional two samples are required for the next 100 ps (for 5 Gb/s), while the rest of the four samples out of ten can effectively describe the trailing 800 ps of the bit duration of 1 Gb/s.
6. Circuit-level simulation and applications
Using the developed SOA numeric modeling tools, we numerically evaluated two SOA-based circuits predominant in all-optical signal processing: (i) an all-optical flip-flop architecture that exploits coupled SOA waveguides, operating on XGM phenomena, and (ii) an all-optical XOR gate that exploits an SOA-MZI configuration and XPM phenomena. The numerical results are compared with the results experimentally obtained, showing very good agreement.
6.1. Coupled SOA flip-flop exploiting XGM phenomena
The proposed SR-Flip Flop relies on the bistability between a slightly and a fully saturated regime of two travelling waveguide SOA-XGM switches . Each SOA is powered by a weak CW at λ1 or λ2 wavelength, respectively. Coupled together through a 70/30 coupler as shown in Figure 8(i), they form a simplified version of an SR-Flip Flop. The coupler defines the coupling efficiency between the two SOAs and is also used for inserting the set/reset pulses to the SOAs and driving the SR-Flip Flop states at the outputs. The 70/30 ratio was found to be the optimum for the trade-off between the SOA coupling efficiency and power ratio at the input/output ports. Exploiting XGM phenomena, one SOA at a time acts as master suppressing the other, which consequently acts as slave. Due to symmetrical setup, the role of master and slave can be interchanged and the state of the SR-Flip Flop is determined by the wavelength of the dominating output. A logical ”1” corresponds to SOA1 being dominant and λ1 suppressing the SOA2 output signal, whereas a logical ”0” is obtained when SOA2 dominates the SR-Flip Flop. Switching between the two states requires injecting proper external set or reset pulses at the dominating SOA through the corresponding branches. When a bit of logical ”1” set signal at λ1 is injected into the dominating SOA1 of logical ”1,” its gain is saturated and the transmission of λ1 CW input signal is blocked. This allows SOA2 gain to recover, unblocking the transmission of CW λ2 and switching the SR-Flip Flop state to logical ”0.” The high λ2 CW now is fed into SOA1 serving as the control signal that suppresses its gain even after the set pulse is extinct. As a result, SOA2 acts as master dominating over the slave SOA1, and the SR-Flip Flop will remain in this state until a reset pulse of λ2 wavelength is fed into SOA2. This will saturate its gain and unblock the transmission of λ1 CW, switching the SR-Flip Flop back to its initial state of logical ”1.”
The flip-flop operation was verified numerically and experimentally with set and reset pulse traces, as shown in Figure 8(ii) and (iii), respectively. The simulated output of the flip-flop at 10 Gb/s for both flip-flop outputs employed the same set/reset patterns as the experimentally employed, shown in Figure 8(iii). The experimental demonstration was performed at a low-operation speed as dictated by the 8.5-m-long fiber-pigtailed coupling stage between the two SOAs, which seems to be the main speed-determining factor . The results of Figure 8(iii) were recorded for λ1 and λ2 CW signals of 150 μW before entering the SOAs and set/reset pulses of 1.8 mW prior to reaching the 50/50 couplers and the corresponding branches of the flip-flop. The two SOAs were driven at 300 mA.
Although the experimental data results confirm successful SR-Flip Flop operation, exhibiting an average extinction ratio (ER) of 7 dB and an amplitude modulation (AM) of less than 2 dB, due to the fiber pigtailed implementation and the large fiber lengths, the operation was limited at low operating speeds. In order to confirm that multi-Gb/s operation speed is possible, we have applied in the numerical evaluation of the proposed SR-Flip Flop the same pattern to the external signals, as the ones employed experimentally, but at 10 Gb/s. CW beams of 150 μW at 1548 nm (λ1) and 1550 nm (λ2) are fed into SOA1 and SOA2, respectively, while the set/reset data streams of 1.8 mW average power each follow the corresponding experimentally employed patterns. The evaluation outputs are depicted in Figure 8(ii), assuming 2-mm intermediate coupling length, and they evidently follow the exact same pattern as the experimentally obtained. For SR-Flip Flop operation with noncomplementary control pulses, they exhibit 13-dB ER and 2.4-dB AM, whereas improved performance of 15-dB ER and 1.3-dB AM is demonstrated in case of complementary set and reset pulses. Although both outputs experience some pattern effect and a slight gain overshoot when switching state, it is clearly shown that full switching can be achieved at 10-Gb/s operational speed.
6.2. SOA-MZI XOR gate exploiting XPM phenomena
An SOA-MZI gate operating as an all-optical XOR gate exploiting XPM phenomena has been both experimentally and numerically evaluated, as shown in Figure 9. Figure 9(i) illustrates the arrangement of the SOA-MZI XOR gate. The XOR-gate has been experimentally and numerically evaluated with the results illustrated in Figure 9(ii) and (iii), respectively. The synchronized numerical time traces demonstrated the proof of principle. Control (Ctr) signal 1 at wavelength λ2 and control 2 at λ3, respectively, are fed as control signals to the SOA-MZI. The time trace obtained at λ1 wavelength at the switching output port of SOA-MZI is illustrated at the third row of numerical traces, where it can be seen that a “logical 1” pulse is obtained, when exclusively one of the input control bits bears a “logical 1,” while a “logical 0” pulse is obtained, when the two input control bits are equal. The experimental eye diagram of the XOR output is illustrated in Figure 9(ii), revealing an extinction ratio (ER) of 8 dB and an amplitude modulation (AM) of 1 dB. Equal performance has been obtained for the simulated eye diagram of the XOR output at Figure 9(ii).
A holistic methodology approach on time domain numerical modeling has been demonstrated, targeting to address accuracy and efficiency. Accuracy is addressed through the development of an experimentally validated numerical model and a gain parameterization procedure. Following the development of a validated numerical model relying on the TMM analysis technique, qualitative results are presented so as to investigate the gain dynamics and the recovery time of the SOA during pump–probe measurements. Efficiency is sought through the development of a newly introduced time domain SOA modeling technique based on the multigrid concepts to introduce adaptive time stepping.
This work was supported by FP7 MC-IAPP project COMANDER (contract no. 612257).