## Abstract

For the first time, an integral equation approach for the numerical assessment of Semiconductor Optical Amplifiers (SOAs) is proposed. Performance comparisons between the suggested formulation and the traditional transfer matrix method are carried out in terms of the computation costs in solving the multi-wave mixing process in bidirectional, high-power SOAs. Computation efficiency improvement by more than an order of magnitude was observed with the proposed formulation, achieving better accuracy at equivalent spatial resolution.

© 2006 Optical Society of America

## 1. Introduction

Steady-state, frequency-domain analysis is one of the most powerful approaches in the precise assessment of nonlinear signal interactions (such as four-wave mixing) in Semiconductor Optical Amplifiers (SOAs). To find the solution of the given problem, traditional approaches rely on the coupled differential equations for the SOAs [1, 2, 3]. In the real implementation, the intensity / phase values of the waves are sequentially calculated along the propagation axis, with the corresponding carrier density distributions at each gain segments of the SOA. When a large number of signal waves or gain segments is necessary for the increased spectral resolution / convergence of the signal, much higher degree of calculation efforts are often required, occasionally with repetitive iteration procedures (e.g., for bi-directional SOAs [4]).

Based on coupled differential equations, other types of optical amplifiers (e.g., fiber Raman amplifier, FRA) [5] traditionally have also experienced same fundamental difficulties in their numerical assessment - slow convergence and limited accuracy. Recently, a meaningful advance has been made for the assessment of FRA solution, drastically reducing the computation cost. Based on a novel integral/matrix formulation for the fiber Raman-amplifier problem, more than two orders of reduction in the computation cost was demonstrated; also enabling the successful application to the inverse problem of FRA gain design [6, 7]. In this approach, solutions were sought along the iteration axis - starting with the analytically determined, seed wave values (defined over all the amplifier segments), and then iteratively multiplying two-dimensional transfer matrices (segment and wave, elements determined from lower-order iteration results).

In this paper, we apply the integral equation formalism to the multi-wave, bi-directional high-power SOA problem including nonlinear interactions, and also address the key differences between SOA and FRA in the equations and solution searching process. Results show greater than 30~500 times improvement in computational efficiency, when compared to the traditional differential equation-based approaches for the SOA solution.

## 2. Formulation

We start from the familiar coupled wave equation for SOA [8],

where coefficients *a _{l}*(

*z*) the complex amplitudes of the signal fields at the frequency

*ω*,

_{j}*z*the propagation axis,

*N*the carrier density,

*g*(

_{l}*N*) the modal gain,

*α*the linewidth enhancement factor, and

*γ*is the scattering loss per unit length.

_{sc}*ε*’s are inverse saturation powers from the nonlinearity,

_{m}*β*’s are equivalent linewidth enhancement factors accounting for gain and index modulation at all frequencies, and

_{m}*τ*’s are relaxation times associated with various nonlinear processes, such as carrier population pulsation (

_{m}*m*=

*cpp*), spectral hole burning (

*m*=

*shb*), and carrier heating (

*m*=

*ch*) [8]. Summation on subscripts ‘

*i*,

*j*,

*k*’ It were made to consider the wave component generated at the frequency

*ω*through the four-wave-mixing effect in the SOA (Δ

_{l}*ω*=

_{ij}*ω*-

_{i}*ω*=

_{j}*ω*-

_{l}*ω*).

_{k}The integral equation formulation is now obtained by integrating Eq. (1) over the propagation axis, *z* ; and is given by

To solve the above equation, we assume an adiabatic process and utilize the iteration procedure. After *n* iterations, Eq. (2) becomes

Further, dividing the SOA length with *x* units of elemental segments (of length Δ*z*, Fig. 1), we now successfully convert Eq. (3) into a matrix form (for *y* number signals for each direction),

where *n* is the iteration number, and * A*,

*, and*

**F***are defined as,*

**T**$$\overrightarrow{{A}^{n}}=\left[\begin{array}{ccc}\overrightarrow{{a}_{1}^{n}}\left(0\right)& \cdots & \overrightarrow{{a}_{1}^{n}}\left(x\Delta z\right)\\ \vdots & \ddots & \vdots \\ \overrightarrow{{a}_{y}^{n}}\left(0\right)& \cdots & \overrightarrow{{a}_{y}^{n}}\left(x\Delta z\right)\end{array}\right],\overleftarrow{{A}^{n}}=\left[\begin{array}{ccc}\overleftarrow{{a}_{1}^{n}}\left(0\right)& \cdots & \overleftarrow{{a}_{1}^{n}}\left(x\Delta z\right)\\ \vdots & \ddots & \vdots \\ \overleftarrow{{a}_{y}^{n}}\left(0\right)& \cdots & \overleftarrow{{a}_{y}^{n}}\left(x\Delta z\right)\end{array}\right]$$

$$\overrightarrow{f\left(z\right)}=\left(\frac{1}{2}{g}_{l}\left(N\left(z\right)\right)\left[\left(1-\mathit{j\alpha}\right)\overrightarrow{{a}_{l}\left(z\right)}-\sum _{m}\left(\sum _{\Delta {\omega}_{\mathit{ij}}}\frac{\left(1-j{\beta}_{m}\right){\epsilon}_{m}}{1+j\Delta {\omega}_{\mathit{ij}}{\tau}_{m}}\overrightarrow{{a}_{i}^{*}\left(z\right)}\overrightarrow{{a}_{j}\left(z\right)}\overrightarrow{{a}_{k}\left(z\right)}\right)\right]-\frac{{\gamma}_{\mathit{sc}}\overrightarrow{{a}_{l}\left(z\right)}}{2}\right)\times \mathrm{\Delta z}$$

$$\overleftarrow{f\left(z\right)}=\left(\frac{1}{2}{g}_{l}\left(N\left(z\right)\right)\left[\left(1-\mathit{j\alpha}\right)\overleftarrow{{a}_{l}\left(z\right)}-\sum _{m}\left(\sum _{\Delta {\omega}_{\mathit{ij}}}\frac{\left(1-j{\beta}_{m}\right){\epsilon}_{m}}{1+j\Delta {\omega}_{\mathit{ij}}{\tau}_{m}}\overleftarrow{{a}_{i}^{*}\left(z\right)}\overleftarrow{{a}_{j}\left(z\right)}\overleftarrow{{a}_{k}\left(z\right)}\right)\right]-\frac{{\gamma}_{\mathit{sc}}\overleftarrow{{a}_{l}\left(z\right)}}{2}\right)\times \mathrm{\Delta z}$$

Based on above construction, the numerical solution of the SOA problem can be obtained from the following process [also refer Fig. 1(b)].

Step 0. Initialize **A**
^{0} for all segments, assuming transparent SOA.

Step I. Calculate **F**
^{0} using **A**
^{0}

Step II. Obtain **A**
^{1}=**F**
^{0}
**T**

Step III. Repeat step I & II for higher numbers of n, with Eq. (6) and **A**
^{n}=**F**
^{(n-1)}
**T**

After a sufficient number of iterations, the final solution **A**^{n} can be obtained, converging to the exact value. It is important to note here that, for the conventional coupled differential equation based approach, solutions are obtained for each segment in a sequential manner, while their values are repeatedly updated as the iteration procedure continues [Fig. 1(a)]. In contrast, the proposed formulation uses the full set of wave information (at every segment and wavelength) from the previous iterations (with **F**^{(n-1)} - enabling a simpler and much faster calculation (note that the information on the SOA segments are updated with simple multiplication of two-dimensional matrices-without the need of solving differential equations for each segment).

The practical implementation of the above process to the SOAs is different from the integral equation for the FRAs, and here we note - 1: the convolution process in Eq. (3) can be treated as a vector product in the Fourier domain (with an additional inverse Fourier transform step) for faster calculation. 2: the carrier density N (absent for FRAs) needs to be calculated for each segment, and 3: for sufficiently small signal, variations in the signal phase have to be ignored (when amplitude of the wave is too small, the phase becomes too sensitive to the amplitude change caused by interactions with other waves, and leads to oscillation).

## 3. Results

To verify the accuracy of the algorithm, we compared our result against the previous report on the SOA multi-wave-mixing process [8]. SOA parameters were taken from the same reference. As can be seen in Fig. 2, perfect match in the result was obtained between the previous approach and proposed integral equation method, over the whole wavelength range. Worth to note, with the Runge-Kutta method, it took approximately 1.5 seconds with a 2-GHz desktop processor to obtain the one-way result. This result is consistent with the calculation time reported in the previous art [8], when considering the difference in the processor speed.

Still, for SOAs operating at higher driving currents and larger input signal power, different response curves were observed for those two methods, when tested with the same number (10) of SOA segments (Fig. 3 : forward/backward input power = 10 / 8 dBm at 1549.2/1550.8 nm. 0.04 nm resolution. SOA input/output coupling loss = 4 dB). Signal power/carrier density distributions in SOA at 500 mA of driving current is also shown in Fig. 4. With the finite difference method, more than 2,000 SOA segments were necessary to reach the final solution. However, with the proposed integral equation method, only 10 SOA segments were necessary to obtain the final, converged solution with much less computation time.

To investigate the convergence behavior, we plot in Fig. 5 the signal power distributions in the SOA at the different stages of the iteration process. The solution obtained by using 2,000 SOA segments with the finite difference method agreed very well with the solution of the integral equation method (10 segments, 20 iterations), and therefore we took it as an exact, reference datum. In contrast to the proposed integral equation method for which the signal power distribution converges to the exact value uniformly over the whole SOA length, the differential equation method exhibited a slow asymptotic convergence to the exact solution at the input section (or output, depending on the direction of solution search) of the SOA. Roughly stated, this is because, with our approach, the update on the wave information at each local segments is made utilizing the previous states (signal power, phase, and carrier density in the transfer matrix) of the SOA, for the whole section of SOA at the same time -rather than sequentially updating the wave along the propagation direction (case of differential equation method), from the states of previous SOA segments. In this sense, the proposed iteration process better mimics the real SOA dynamics, and thus considered to exhibit much efficient convergence with reduced error.

Finally, to detail the relationship between the convergence/computation cost and the number of SOA segments used in the simulation, we plot in Fig. 6 the output power convergence error and the required computation time as a function of SOA segmentation (The SOA operating condition is identical to that of Figs. 3, 4, and 5. output power = 14.8 dBm/13.8 dBm for forward/backward signal). As clearly can be seen, with the integral method, the output power of the SOA converges to its limiting value much faster (0.7 seconds, error < 0.0012dB), with a much less number of SOA segments (6), when compared to the previous finite difference method (352 seconds, error < 0.0012dB, with 2,000 segments). More than five hundred times faster, and precise numerical assessment of the result was achieved. It is worth to note that, even when compared to the differential equation solving based on higher orders of approximation/coding effort (Runge-Kutta with 4^{th} order), ~30 times faster assessment of the result was possible when using the proposed approach.

## 4. Conclusion

For the first time, we demonstrated that it is possible to apply the integral equation method to the analysis of multi-wave, high power, bi-directional SOAs. Solutions were sought along the iteration axis from the analytically-determined seed wave and subsequential multiplication of transfer matrices. Comparisons are carried out in terms of computation expenses between the suggested formulation and traditional differential-equation approach. Results show greater than 30~500 times improvement in the computation efficiency for the proposed method, when compared to the traditional/improved Runge-Kutta coupled differential equation methods. The application of the proposed formulation can be found in the faster and precise characterization of SOA response curves, for example, calculating the L-I curve of high-power SOAs, which involve wide range (numbers) of operating condition/driving currents.

## References and links

**1. **C. Y. J Chu and H. Ghafouri-Shiraz, “Analysis of gain and saturation characteristics of a semiconductor laser optical amplifier using transfer matrices,” J. Lightwave Technol. **12**, 1378–1386 (1994). [CrossRef]

**2. **S. L. Zhang and J. J. O’Reilly, “Modelling of four-wave-mixing wavelength conversion in a semiconductor laser amplifier,” Physical Modelling of Semiconductor Devices, IEE Colloquium on, 5/1-5/6 (1995).

**3. **H. Lee, H. Yoon, Y. Kim, and J. Jeong, “Theoretical study of frequency chirping and extinction ratio of wavelength-converted optical signals by XGM and XPM using SOA’s,” IEEE J. Quantum Electron. **35**, 1213–1219 (1999). [CrossRef]

**4. **M. J. Connelly, “Wideband semiconductor optical amplifier steady-state numerical model,” IEEE J. Quantum Electron. **37**, 439–447 (2001). [CrossRef]

**5. **B. Min, W. J. Lee, and N. Park, “Efficient formulation of Raman amplifier propagation equations with average power analysis,” IEEE Photonics Technol. Lett. **12**, 1486–1488 (2000). [CrossRef]

**6. **N. Park, P. Kim, J. Park, and L. K. Choi, “Integral form expansion of fiber Raman amplifier problem,” Opt. Fiber Technol., **11**, 111–130 (2005). [CrossRef]

**7. **J. Park, P. Kim, J. Park, H. Lee, and N. Park, “Closed integral form expansion of Raman equation for efficient gain optimization process,” IEEE Photonics Technol. Lett. **16**, 1649–1651 (2004). [CrossRef]

**8. **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]