In this paper, we have implemented and compared two complementary time-domain models that have been widely used for the simulation of SOAs. One of the key differences between them lies in their treatment of the material (gain and refractive index) dispersion. One model named as a spectrum slicing model (SSM) is desirable for the simulation of broadband behaviours of SOAs, but not for the nonlinear effect such as the intermodulation distortion, since the gain dispersion is considered by slicing the entire spontaneous emission spectrum into many stripes. The other model based on effective Bloch equations (EBE’s) is capable of dealing with the SOA nonlinear effects with the material dispersion incorporated explicitly through the susceptibility, but can’t capture the broadband behaviours. Both of them, however, can readily handle the SOA characteristics such as the fibre-to-fibre gain, noise, and crosstalk. Through a direct comparison between them, we have shown that they are in generally good agreement. A discussion on detailed implementations and each model’s salient features is also presented.
©2006 Optical Society of America
Numerical modeling and simulation are essential to predict and characterize the complex behaviours of laser diodes (LDs) and semiconductor optical amplifiers (SOAs) operating under stringent conditions. For the simulation of laser diodes  and vertical-cavity semiconductor optical amplifiers (VCSOAs) , , the classical rate equations and the Fabry-Perot (FP) resonator approaches have been employed. However, the neglect of the longitudinal spatial hole burning (LSHB) effect in the spatial-independent rate equations results in the incorrect predictions. The standing-wave model for the VCSOAs is valid only when the photon lifetime is much shorter than the characteristic time of the laser dynamics . The traveling-wave models for LDs ,  have been employed without making any assumptions about the cavity modes. However, those models may not be applicable or need to be extended further for the simulation of broadband behaviours of SOAs in which the parasitic ASE noise appears over a broad wavelength range.
The full microscopic model  deals with the Maxwell-Bloch and carrier transport equations in a self-consistent manner with the thermal dynamics included. Although it is fully physics-based and complete the most, the calculation of the macroscopic quantities entails a lot of computational efforts. In Ref. , the wave equation and the Bloch equation are so combined into one integro-differential equation as to treat the interaction between the optical field and polarization. Nonuniform gain over a very broad wavelength range is also considered. However, much care is taken to ensure that the least-square fitting for the material gain is accurate within a certain wavelength range. The models in Refs. – consider the material gain dynamics by expanding the linear gain and refractive index in the frequency domain, thereby leading to the second-order time derivative in the traveling-wave equations. The model in Refs. ,  treats the traveling-wave equation and the effective Bloch equation (EBE). This EBE model can handle the gain dynamics and material dispersion, yet requires less computational efforts. The polarization dynamics is extracted from the Lorentzian fit to the microscopic gain/index –. The main peculiarity of this model is an excellent balance between the accuracy and computation efficiency. However, some difficulty is encountered when the multiple channels are spread over a wide wavelength range. In the case where the channel detuning is large, one has to inevitably use an extremely small time step to handle the fast optical field dynamics (a fast oscillatory envelope variation). Moreover, when the optical field dynamics is faster than the system response speed, this model may not be valid. Finally, the model in Refs. – solves a number of traveling-wave equations for the simulation of SOA broadband behaviours. Similarly, the model in Ref  named as the spectrum slicing model (SSM) incorporates the gain dispersion by slicing the entire spontaneous emission spectrum into many strips (=bands) and assigning the traveling-wave equations into each band. Therefore, it is valid only when the signal bandwidth is so narrow that the material gain and refractive index are almost constant (frequency-independent) within each strip.
Those different time-domain models – are all based on a rigorous formulation. Therefore, the propagation and reflection of the optical fields in the cavity, carrier and gain dynamics, and the spatial variation of carrier density and field intensity can all be taken into consideration with the ASE noise contribution further included. Hence, they are applicable to the SOA performance simulation. Of those models, the EBE is desirable for the simulation of SOA nonlinear effects (e.g., channel crosstalk and intermodulation distortion), since the material gain and index dispersions are incorporated explicitly through the susceptibility; whereas the SSM can readily capture the SOA broadband behaviours that are hard to be treated by the EBE, as the entire ASE spectrum is discretized. In this paper, we make a systematic study exclusively on these two complementary models. This will shed light on some of the salient features of these models and provide useful guidelines for the applications of them. The SSM and its implementation are described in full detail in Ref. . However, a detailed computation technique with the EBE, especially in the context of nonlinear phenomena such as the cross-gain saturation (XGS, the origin of channel crosstalk) and nondegenerate four-wave mixing (FWM, the origin of intermodulation distortion), is firstly presented in this paper.
This paper is written with the following objectives. Firstly, we clarify the unique features of each model and its limitations in approach to various problems. Secondly, we make the best use of the EBE to capture the SOA nonlinear effects and discuss a detailed technique for quantifying them. Finally, we make a quantitative comparison between them by simulating typical issues of SOAs.
2.1 Effective Bloch equations (EBE’s)
By substituting the 1D optical field into Maxwell’s wave equations and applying the slow-varying envelope approximation (SVEA), we can obtain the following traveling-wave equations:
where denote the forward and backward coupled polarizations, respectively. ε o,r are the permittivity in free space and the average relative permittivity of the optical waveguide, respectively. ν g is the group velocity, αs is the modal loss, Γ is the confinement factor, β is the propagation constant, and E (f,r) are the forward and backward traveling optical waves, respectively. The material dispersion is considered through the frequency-domain susceptibility χ given as , :
with q defined as the magnitude of the electron charge, h the Planck constant, μ¯ the refractive index of the medium, and m 0 the free-electron mass. M if represents the momentum matrix element that links the transition between the initial and final states. The variables ρc and ρν are the densities of states per unit volume per unit energy for the conduction and valence bands, respectively. f c and f ν are the Fermi factors for electrons and holes, respectively. The carrier-induced refractive index change δ n is computed by 
where P indicates taking the principal value integral.
where χo denotes the frequency-independent (background) susceptibility, δn0 the frequency-independent (background) carrier-induced index change, go the frequency-independent (background) material gain, ωo the reference angular frequency, and ωp the material gain peak position. The mapping is realized by searching for the best fitting between Eq. (2) and Eq. (5a) with A i, δi , and Γi taken as the carrier-dependent Lorentzian fitting parameters. Only one Lorentzian oscillator (T=1) is enough to approximate the material dispersion of SOAs , .
for i=1, 2,… , T. It is noted that in obtaining Eq. (6) by taking the Fourier transformation of Eq. (5), those Lorentzian fitting parameters are assumed to be independent of time. Such an assumption may be valid as long as the time marching step is so small that the carrier density varies much more slowly compared to the inverse gain bandwidth , . We shall return to this point later with the simulation result. The ASE noise with a random phase is treated by a Gaussian-distributed random number generator  satisfying the following condition:
where d z indicates the length of a subsection introduced by the spatial discretization of the active region along the wave propagation direction and γ the coupling coefficient of the spontaneous emission rate R sp into the waveguide. δ() denotes the Dirac-function.
Finally, the carrier rate equation can be written as
where η denotes the current injection efficiency, A the non-radiative (Shockley-Reed-Hall and surface) recombination coefficient, B the bimolecular recombination coefficient, and C the Auger recombination coefficient.
2.2 Spectrum slicing model (SSM)
In the SSM, the broadband behaviors of SOAs are described by slicing the entire spontaneous emission spectrum into N d narrow segments with equal width that cover the whole ASE spectrum (1.2μm~1.65μm). Since it is sliced in the frequency domain, the model is also named as the mixed frequency-time-domain traveling-wave model  written as;
for k=1, 2,… , N d.
The spontaneous emission rate is computed by 
where ΔνL is the subsection width (measured in frequency) sliced in the ASE spectrum. It is shown in Ref.  that the selection of ΔνL is a tradeoff between the accuracy and computation efficiency. However, the total noise contribution to the gain saturation won’t be affected by the selection of ΔνL .
All the optical waves at different wavelength segments in Eq. (9) are coupled through the “carrier-sharing” mechanism, governed by the carrier rate equation:
3. Numerical implementation
To facilitate the simulation, the active layer is spatially divided into M (=61) subsections. The longitudinal variation of carrier density and optical intensity is thus taken into account. The structural and material parameter values used in the simulations are summarized in Table 1.
As discussed in Ref. , the device gain obtained with the SSM is very sensitive to ΔνL for low input signal powers. In fact, the device gain is overestimated for larger ΔνL (coarse discretization), as more spontaneous emission noise is coupled into the low input signal power. Therefore, one needs an adaptive discretization scheme; namely, ΔνL is kept reducing for low input signal powers until a certain convergence is reached.
In the EBE, all of the susceptibility fitting parameters appearing in Eq. (6b) are obtained as a function of carrier density through fitting Eq. (5a) into Eq. (2) using the Levenberg-Marquardt algorithm . Since it is very hard to obtain the complete Lorentzian fit over the entire gain spectrum, we make the best fit around the gain peaks as done in Refs.  and . Though the fit is generally good enough, a numerical error may be inevitably introduced in the EBE. It is noteworthy that it may be difficult to obtain a good fit for laser diodes, as the carrier density may be clamped at a low level. However, this is not the case for SOAs. As also discussed in 2.1, those fitting parameters are assumed to be independent of time. Therefore, an additional numerical error may be parasitic in the EBE unless a time step (dt) is infinitely small. To minimize and possibly ignore such a numerical error, therefore, one has to use a very small time step. As demonstrated in Fig. 1, showing the dependence of the numerical result on a time step, the output power converges into some value as dt decreases and thus the numerical error may become negligibly small as long as dt is small enough (<0.13ps).
As shortly aforementioned, when the multichannel signals are spread afar over a wide wavelength range, the optical field (envelope) dynamics becomes fast as shown in Fig. 2 where the channel spacing between two different signals is changed from 400MHz to 2GHz for the input signal consisting of two identical Gaussian pulses with 10ns 1/e-width. To deal with such fast field dynamics, one also has to use a small time step. Figure 3 shows the average output envelope power as a function of the channel spacing when the average input envelope power [the same envelope shown in Fig. 2(a)] is -13.63dBm. At this input power level, the SOA provides the device gain as high as 20dB (also seen in Fig. 8). With a larger time step, the device gain increases sharply (up to 28dB) as the channel spacing increases. With a smaller time step, however, the numerical result is less sensitive to the channel spacing. In short, the numerical accuracy of the SSM and EBE depends sensitively on ΔνL and dt, respectively.
4. Results and discussions
4.1 Unique feature of SSM: broadband behaviours
In the simulation of SOA broadband behaviours, the broadband gain dynamics needs to be considered. It can readily be incorporated in the SSM, since the entire spontaneous emission spectrum is discretized and a number of traveling-wave equations are solved accordingly. To demonstrate it, we have computed the ASE spectra of a SOA and gain-clamped (GC-) SOA under the bias current of 140mA and presented the results in Fig. 4. In the GC-SOA, the lasing oscillation is generated at 1514nm by DBR gratings (200μm in length) with the grating coupling coefficient of 6cm-1. The ASE noise is shown to be considerably suppressed by the internal lasing mode. In addition, the ASE peak is red-shifted. In typical GC-SOAs, the lasing mode is normally positioned far away (more than 20nm blue-shift) from the signal band (near gain peak). In such a case, it is not feasible to capture the interaction between the lasing mode and signals with the EBE, unless additional Bloch equations are introduced at the lasing wavelength.
4.2 Salient features of EBE: pulse propagation and intermodulation distortion
As discussed earlier, the gain and refractive index within one sliced band are constant (frequency independent) in the SSM, while the material dispersion in the EBE is considered explicitly through the susceptibility in Eq. (2). For this reason, we can carry out the simulation of short-pulse propagation only with the EBE. Figure 5 shows the pulse distortion through the SOA depending on the pulse width (τ0 ). It is evident that pulses with a width shorter than the carrier lifetime (≈1/A) are distorted the most .
Due to the discrete nature of the SSM in the frequency domain, it can’t be applicable for the case where the signal bandwidth is wider than the channel spacing so that there is an overlap in the output spectrum. In this situation, only the EBE can correctly capture the beatings between the signals. This is demonstrated in Fig. 6, showing the input and output signal spectra of the SOA calculated with the EBE when the input signal envelope consists of two identical Gaussian pulses with 4.5ns 1/e-width in different channels. In each output signal, it is apparent that the nonlinear effects such as the intermodulation distortion and crosstalk are mixed up.
The nonlinear effect such as the intermodulation distortion (IMD) that originates from the nondegenerate four-wave mixing (FWM)  can also be treated only with the EBE for the same reason stated above. To simulate it, two Gaussian pulses with 10ns 1/e-width in different channels (dual-channel) are launched into the SOA. Figure 7(a) is the spectra of the input and output envelopes for the input power of -15dBm in each channel (channel spacing=400MHz). Due to the nondegenerate FWM of SOAs, intermodulation distortion products appear in the output spectrum. For a quantitative study, we calculate the relative third-order IMD defined in Ref. . Shown in Fig. 7(b) is the computed IMD as a function of the Ch2 input power for the Ch1 input power of -20dBm. As generally known, the relative IMD is getting severe (increases very fast) as the total input power increases.
4.3 Comparison: gain and noise figure
The SOA performances such as the device gain and noise figure are readily computed with both models. Therefore, it is of theoretical interest to make a direct comparison between them. Figure 8 shows the fibre-to-fibre gain and noise figure with respect to the input signal powers. A good agreement has been obtained except for a little discrepancy in the noise figure. This discrepancy may be due to an overestimation of the noise contribution in the SSM. In fact, after the spectrum slicing in the SSM, the beating among noises in each strip is considered coherently. In the EBE, however, the noise contribution from the entire spectrum is added incoherently through the pseudo-random phase.
4.4 Comparison: crosstalk
The crosstalk that originates from the XGS , ,  can be treated by both models. However, the way of calculation is slightly different. In the SSM, the input power of Ch1 is kept constant while the input power of Ch2 is modulated at a certain frequency. Therefore, the crosstalk is defined as the ratio of the spurious modulation appearing in the time-domain output envelope in Ch1 to the intentionally imposed modulation in Ch2 . Since only one oscillatory signal envelope is excited in the EBE, the time-domain output envelope in each channel can be recovered only by taking inverse Fourier transform of the associated output spectrum. For instance, the rectangular window marked in Fig. 9 is inversely Fourier transformed to recover the time-domain output pulse in Ch1. Noticing the fact that the Gibb’s phenomenon would be so pronounced in the output spectrum when the constant input power is excited in Ch1, we rather inject a Gaussian pulse. This enables us to recover the time-domain output pulse (still in Gaussian shape) with more accuracy. At the same time, we inject the constant input power or modulated Gaussian-pulse stream in Ch2, which is identical with the one in Ch1. The average input power of the Gaussian pulse is equal to the constant input power. The crosstalk from Ch2 to Ch1 is therefore defined as the relative output power difference in Ch1 when the input power in Ch2 is modulated and constant at the same power level.
In this subsection, the nonlinear gain due to the carrier heating (CH) and spectral hole burning (SHB) is taken into consideration to identify its effect on the crosstalk. To this end, we simply multiply the material gain in Eq. (3) by (1-ε 1 S)  where ε1 is the gain compression coefficient and S the total photon density inside the cavity. Shown in Fig. 10(a) is the crosstalk versus the Ch2 input power at 1535.7nm when the Ch1 input power at 1538.9nm is -23dBm and I=140mA. The crosstalk increases with increasing Ch2 input power and reaches saturation at about -5dBm. A small discrepancy between the models may come from the different way of treating the material dispersion or slightly different way of calculating the crosstalk. The calculation of the crosstalk with the SSM is done directly from the time-domain output envelope, while its calculation with the EBE requires one more step, an inverse Fourier transform, to get the envelope. Another finding is that, although the effect of the nonlinear gain is not so remarkable, it indeed causes more crosstalk. Figure 10(b) is the crosstalk as a function of the detuning frequency (channel spacing) when the input power at each channel is -10dBm. Both models show that the crosstalk keeps unchanged as the detuning frequency increases, indicating that the channel crosstalk is attributed to the cross-gain saturation, a nonlinear phenomenon that doesn’t depend on the beating frequency (i.e., the channel spacing) explicitly .
Of the existing time-domain models, two complementary models (the spectrum slicing model and the effective Bloch equation model) have been implemented and compared. The motivation is to clarify the advantages and shortcomings of each model and utilize both of them in approach to various problems. The key difference between them lies in their treatment of the material dispersion. The SSM considers the gain dispersion by slicing the entire ASE spectrum into many stripes. It can therefore treat the broadband behaviours and the channel crosstalk of SOAs, but not the IMD. With the material dispersion included explicitly through the frequency-dependent susceptibility, the EBE is capable of dealing with both IMD and crosstalk, but not the broadband behaviours. The EBE can also handle the short-pulse propagation and the situation where the signal bandwidth is wider than the channel spacing. However, it requires a smaller time step to minimize a numerical error when the frequency detuning among signal channels is getting larger. We have made a direct comparison between them in terms of the device gain, noise figure, and crosstalk, and demonstrated that they are in generally good agreement. Both models have also shown that the nonlinear gain indeed increases the crosstalk, but its effect is not so pronounced.
The authors would like to thank Dr. Xun Li of McMaster University, Hamilton, Canada, for his valuable discussions.
References and links
1. T. J. Menne, “Analysis of the uniform rate equation model of laser dynamics,” IEEE J. Quantum Electron. 2, 38–44 (1966). [CrossRef]
2. C. Tombling, T. Saitoh, and T. Mukai, “Performance predictions for vertical-cavity semiconductor laser amplifiers,” IEEE J. Quantum Electron. 30, 2491–2499 (1994). [CrossRef]
3. J. Piprek, S. BjÖrlin, and J. E. Bowers, “Design and analysis of vertical-cavity semiconductor optical amplifiers,” IEEE J. Quantum Electron. 37, 127–134 (2001). [CrossRef]
4. W. Li, W.-P. Huang, X. Li, and J. Hong, “Multiwavelength gain-coupled DFB laser cascade: design modeling and simulation,” IEEE J. Quantum Electron. 36, 1110–1116 (2000). [CrossRef]
5. L. M. Zhang, S. F. Yu, M. Nowell, D. D. Marcenac, and J. E. Carroll, “Dynamic analysis of radiation and side mode suppression in second-order DFB lasers using time-domain large signal traveling wave model,” IEEE J. Quantum Electron. 30, 1389–1395 (1994). [CrossRef]
6. A. J. Lowery, “New dynamic semiconductor laser model based on the transmission line modeling method,” IEE Proc. J. 134, 281–289 (1987).
7. E. Gehrig, O. Hess, and R. Wallenstein, “Modeling of the performance of high-power diode amplifier systems with an optothermal microscopic spatio-temporal theory,” IEEE J. Quantum Electron. 35, 320–331 (1999). [CrossRef]
8. M. Kolesik and J. V. Moloney, “A spatial digital filter method for broadband simulation of semiconductor lasers,” IEEE J. Quantum Electron. 37, 936–944 (2001). [CrossRef]
9. G. P. Agrawal, “Effect of gain dispersion on ultrashort pulse amplification in semiconductor laser amplifiers,” IEEE J. Quantum Electron. 27, 1843–1849 (1991). [CrossRef]
11. M. Homar, J. V. Moloney, and M. San Miguel, “Traveling wave model of a multimode Fabry-Perot laser in free running and external cavity configurations,” IEEE J. Quantum Electron. 32, 553–566 (1996). [CrossRef]
12. G. C. Dente and M. L. Tilton, “Modeling multiple-longitudinal-mode dynamics in semiconductor lasers,” IEEE J. Quantum Electron. 34, 325–335 (1998). [CrossRef]
13. T. Durhuus, B. Mikkelsen, and K. E. Stubkjaer, “Detailed dynamic model for semiconductor optical amplifiers and their crosstalk and inter-modulation distortion,” J. Lightwave Technol. 10, 1056–1065 (1992). [CrossRef]
14. A. Mecozzi and J. Mork, “Saturation effects in nondegenerate four-wave mixing between short optical pulses in semiconductor laser amplifiers,” IEEE J. Sel. Top. Quantum Electron. 3, 1190–1207 (1997). [CrossRef]
15. C. Z. Ning, R. A. Indik, and J. V. Moloney, “Effective Bloch equations for semiconductor lasers and amplifiers,” IEEE J. Quantum Electron. 33, 1543–1550 (1997). [CrossRef]
16. C. Z. Ning, J. V. Moloney, A. Egan, and R. A. Indik, “A first-principle fully space-time resolved model of a semiconductor laser,” Quantum Semiclassic. Opt. 9, 681–691 (1997). [CrossRef]
17. U. Bandelow, M. Radziunas, J. Sieber, and M. Wolfrum, “Impact of gain dispersion on the Spatio-temperal dynamics of multisection lasers,” IEEE J. Quantum Electron. 37, 183–188 (2001). [CrossRef]
18. M. Bahl, H. Rao, N. C. Panoiu, and R. M. Osgood, Jr, “Simulation of mode-locked surface-emitting lasers through a finite-difference time-domain algorithm,” Opt. Lett. 29, 1689–1691 (2004). [CrossRef] [PubMed]
19. M. A. Summerfield and R. S. Tucker, “Frequency-domain model of multiwave mixing in bulk semiconductor optical amplifiers,” IEEE J. Sel. Top. Quantum Electron. 5, 839–850 (1999). [CrossRef]
20. M. J. Connelly, “Wideband semiconductor optical amplifier steady-state numerical model,” IEEE J. Quantum Electron. 37, 439–1103 (2001). [CrossRef]
21. M. J. Connelly, “Wideband dynamic numerical model of a tapered buried ridge stripe semiconductor optical amplifier gate,” IEE Proc.: Circuits Devices Syst. 149, 173–178 (2002). [CrossRef]
22. J. W. Park, X. Li, and W. P. Huang, “Comparative study on mixed frequency-time-domain models of semiconductor laser optical amplifiers,” IEE Proc.: Optoelectron. 152, 151–159 (2005). [CrossRef]
23. G. P. Agrawal and N. K. Dutta, “Semiconductor Lasers,” (Van Nostrand Reinhold, New York, 1993).
24. C. H. Henry, R. A. Logan, and K. A. Bertness, “Spectral dependence of the change in refractive index due to carrier injection in GaAs lasers,” J. Appl. Phys. 52, 4457–4461 (1981). [CrossRef]
25. W. H. Press, B. P. Flannery, S. A. Teukolssy, and W. T. Vetterling, “Numerical Recipes: The art of Scientific Computing,” (Cambridge Univ. Press, Cambridge, MA, 1986).
26. G. P. Agrawal, “Fiber-optic communication systems,” 3rd edition, (Wiley-Interscience, 2002). [CrossRef]
27. J. Sun, G. Morthier, and R. Baets, “Numerical and theoretical study of the crosstalk in gain clamped semiconductor optical amplifiers,” IEEE J. Sel. Top. Quantum Electron. 3, 1162–1167 (1997). [CrossRef]
28. H. E. Lassen, P. B. Hansen, and K. E. Stubkjaer, “Crosstalk in 1.5μm InGaAsP optical amplifiers,” J. Lightwave Technol. 6, 1559–1565 (1988). [CrossRef]