**2.2. Mathematical Formulation of Modified Nonlinear Schrödinger Equation (MNLSE)**

In this subsection, we will briefly explain the theoretical analysis of short optical pulses propagation in SOAs. We start from Maxwell's equations (Agrawal, 1995; Yariv, 1991; Sauter, 1996) and reach the propagation equation of short optical pulses in SOAs, which are governed by the wave equation (Agrawal & Olsson, 1989) in the frequency domain.

$$
\nabla^2 \overline{E}(\mathbf{x}, y, z, \alpha) + \frac{\mathcal{E}\_r}{c^2} \alpha^2 \overline{E}(\mathbf{x}, y, z, \alpha) = 0 \tag{1}
$$

where, *Exyz* (,,, ) ω is the electromagnetic field of the pulse in the frequency domain, *c* is the velocity of light in vacuum and *<sup>r</sup>* ε is the nonlinear dielectric constant which is dependent on the electric field in a complex form. By slowly varying the envelope approximation and integrating the transverse dimensions we arrive at the pulse propagation equation in SOAs (Agrawal & Olsson, 1989; Dienes *et al.,* 1996).

$$\frac{\partial V(o\nu, z)}{\partial z} = -i \left\{ \frac{o\nu}{c} \left[ 1 + \mathcal{X}\_{\text{u}}(o\nu) + \Gamma \tilde{\mathcal{X}}(o\nu, N) \right] \Big| \Big{Y} - \beta\_{\text{o}} \right\} V(o\nu, z) \tag{2}$$

where, *V z* ( ,) ω is the Fourier-transform of *Vtz* (, ) representing pulse envelope, ( ) *<sup>m</sup>*χ ω is the background (mode and material) susceptibility, χ( ) ω is the complex susceptibility which represents the contribution of the active medium, *N* is the effective population density, 0 β is the propagation constant. The quantity Γ represents the overlap/ confinement factor of the transverse field distribution of the signal with the active region as defined in (Agrawal & Olsson, 1989).

Using mathematical manipulations (Sauter, 1996; Dienes *et al.,* 1996), including the real part of the instantaneous nonlinear Kerr effect as a single nonlinear index n2 and by adding the TPA term we obtain the MNLSE for the phenomenological model of semiconductor laser and amplifiers (Hong *et al.,* 1996). The following MNLSE (Hong *et al.,* 1996; Das *et al.,* 2000) is used for the simulation of FWM characteristics and optical phase-conjugation characteristics with input pump and probe pulses in SOAs.

$$\begin{split} \left[\frac{\partial}{\partial z} - \frac{i}{2} \mathcal{B}\_2 \frac{\partial^2}{\partial \tau^2} + \frac{\mathcal{Y}}{2} + \left(\frac{\mathcal{Y}\_{2p}}{2} + ib\_2\right) \left| V(\tau, z) \right|^2 \right] V(\tau, z) \\ = \left[\frac{1}{2} g\_N(\tau) \left[\frac{1}{f(\tau)} + i\alpha\_N \right] + \frac{1}{2} \Delta g\_\Gamma(\tau) (1 + i\alpha\_\Gamma) - i \frac{1}{2} \frac{\partial g(\tau, \alpha \theta)}{\partial \alpha} \bigg|\_{\alpha\_0} \frac{\partial}{\partial \tau} - \frac{1}{4} \frac{\partial^2 g(\tau, \alpha \theta)}{\partial \alpha^2} \bigg|\_{\alpha\_0} \frac{\partial^2}{\partial \tau^2} \right] V(\tau, z) \end{split} \tag{3}$$

Here, we introduce the frame of local time τ (=*t* - *z*/ *<sup>g</sup> v* ), which propagates with a group velocity *<sup>g</sup> v* at the center frequency of an optical pulse. A slowly varying envelope approximation is used in (3), where the temporal variation of the complex envelope function is very slow compared with the cycle of the optical field. In (3), *V z* (,) τ is the time domain complex envelope function of an optical pulse, <sup>2</sup> *V z* (,) τ represents the corresponding optical pulse power, and 2 β is the GVD. γ is the linear loss, 2*<sup>p</sup>* γ is the TPA coefficient, 2 *b* (= 0 2 ω *n /cA*) is the instantaneous SPM term due to the instantaneous nonlinear Kerr effect 2 0 *n* , ω *(= 2*π 0 *f )* is the center angular frequency of the pulse, *c* is the velocity of light in vacuum, *A (= wd/*Γ*)* is the effective area (*d* and *w* are the thickness and width of the active region, respectively and Γ is the confinement factor) of the active region.

The saturation of the gain due to the CD is given by (Hong *et al.,* 1996)

$$\log\_N(\tau) = \mathcal{g}\_0 \exp\left(-\frac{1}{W\_s} \int\_{-\infty}^{\tau} e^{-s/\tau\_s} \left| V(s) \right|^2 ds \right) \tag{4}$$

where, ( ) *Ng* τ is the saturated gain due to CD, 0 *g* is the linear gain, *Ws* is the saturation energy, *<sup>s</sup>* τis the carrier lifetime.

The SHB function *f*(τ) is given by (Hong *et al.,* 1996)

88 Optical Communication

equation (MNLSE).

where, *Exyz* (,,, )

where, *V z* ( ,) ω

density, 0 β

(Agrawal & Olsson, 1989).

β

τ

ω

velocity of light in vacuum and *<sup>r</sup>*

(Agrawal & Olsson, 1989; Dienes *et al.,* 1996).

ω

**(MNLSE)** 

modeling/ simulation and the mathematical formulation of modified nonlinear Schrödinger

In this subsection, we will briefly explain the theoretical analysis of short optical pulses propagation in SOAs. We start from Maxwell's equations (Agrawal, 1995; Yariv, 1991; Sauter, 1996) and reach the propagation equation of short optical pulses in SOAs, which are

> <sup>2</sup> (,,, ) (,,, ) 0 *<sup>r</sup> Exyz Exyz c* ε

on the electric field in a complex form. By slowly varying the envelope approximation and integrating the transverse dimensions we arrive at the pulse propagation equation in SOAs

> <sup>2</sup> ( ,) 1 ( ) ( , ) ( ,) *<sup>m</sup> V z <sup>i</sup> N Vz*

which represents the contribution of the active medium, *N* is the effective population

factor of the transverse field distribution of the signal with the active region as defined in

Using mathematical manipulations (Sauter, 1996; Dienes *et al.,* 1996), including the real part of the instantaneous nonlinear Kerr effect as a single nonlinear index n2 and by adding the TPA term we obtain the MNLSE for the phenomenological model of semiconductor laser and amplifiers (Hong *et al.,* 1996). The following MNLSE (Hong *et al.,* 1996; Das *et al.,* 2000) is used for the simulation of FWM characteristics and optical phase-conjugation

 χω

is the Fourier-transform of *Vtz* (, ) representing pulse envelope, ( ) *<sup>m</sup>*

is the propagation constant. The quantity Γ represents the overlap/ confinement

11 1 1 1 (, ) (, ) ( ) ( )(1 ) (,) 2 () 2 <sup>2</sup> <sup>4</sup>

∂ ∂ ∂ ∂ <sup>=</sup> + +Δ + − <sup>−</sup> ∂ ∂ ∂ ∂

τ

velocity *<sup>g</sup> v* at the center frequency of an optical pulse. A slowly varying envelope

*g g g i g ii V z*

τ ω

ω

ω

<sup>∂</sup> =− + +Γ − <sup>∂</sup>

 + = ω

 ω

1

χ( ) ω

is the nonlinear dielectric constant which is dependent

0

ω

0 0

(=*t* - *z*/ *<sup>g</sup> v* ), which propagates with a group

ω  τ 2 2 2 2

τ ω

ω

ω

 τ τ

(3)

(2)

 β

is the electromagnetic field of the pulse in the frequency domain, *c* is the

(1)

χ ωis

is the complex susceptibility

**2.2. Mathematical Formulation of Modified Nonlinear Schrödinger Equation** 

governed by the wave equation (Agrawal & Olsson, 1989) in the frequency domain.

2 2

χ

 ∇ ω

ε

 ω

*z c*

characteristics with input pump and probe pulses in SOAs.

*N NT T*

τ

 τ

 τα

<sup>2</sup> <sup>2</sup> <sup>2</sup>

*p*

γ γ

*<sup>i</sup> ib V z V z <sup>z</sup>*

(,) (,) 2 22

 α

2 2 2

∂ ∂ − ++ + <sup>∂</sup> <sup>∂</sup>

τ

*f*

Here, we introduce the frame of local time

τ

the background (mode and material) susceptibility,

$$f(\tau) = 1 + \frac{1}{\tau\_{shb} P\_{shb}} \int\_{-\infty}^{+\infty} \mu(\mathbf{s}) e^{-s/\tau\_{shb}} \left| V(\tau - s) \right|^2 ds \tag{5}$$

where, *f*(τ) is the SHB function, *shb P* is the SHB saturation power, *shb* τ is the SHB relaxation time, and α *<sup>N</sup>* and α*<sup>T</sup>* are the line width enhancement factor associated with the gain changes due to the CD and CH.

The resulting gain change due to the CH and TPA is given by (Hong *et al.,* 1996)

$$\begin{split} \Delta g\_{T}(\boldsymbol{\tau}) &= -h\_{1} \int\_{-\infty}^{+\infty} \mu(\boldsymbol{s}) e^{-s/\tau\_{\rm{ch}}} \left( 1 - e^{-s/\tau\_{\rm{sh}}} \right) \left| V(\boldsymbol{\tau} - \boldsymbol{s}) \right|^{2} \, d\boldsymbol{s} \\ &- h\_{2} \int\_{-\infty}^{+\infty} \mu(\boldsymbol{s}) e^{-s/\tau\_{\rm{ch}}} (1 - e^{-s/\tau\_{\rm{sh}}}) \left| V(\boldsymbol{\tau} - \boldsymbol{s}) \right|^{4} \, d\boldsymbol{s} \end{split} \tag{6}$$

where, ( ) *<sup>T</sup>* Δ*g* τ is the resulting gain change due to the CH and TPA, *u*(s) is the unit step function, *ch* τ is the CH relaxation time, 1 *h* is the contribution of stimulated emission and free-carrier absorption to the CH gain reduction and 2 *h* is the contribution of TPA.

The dynamically varying slope and curvature of the gain plays a shaping role for pulses in the sub-picosecond range. The first and second order differential net (saturated) gain terms are (Hong *et al.,* 1996),

$$\left. \frac{\partial \mathcal{g}(\mathbf{r}, a\mathbf{0})}{\partial a} \right|\_{a\mathbf{0}} = A\_1 + B\_1 \left[ \mathcal{g}\_0 - \mathcal{g}(\mathbf{r}, a\mathbf{0}\_0) \right] \tag{7}$$

$$\left.\frac{\partial^2 g(\mathbf{r}, \alpha)}{\partial \alpha^2}\right|\_{\alpha\_0} = A\_2 + B\_2 \left[g\_0 - g(\mathbf{r}, \alpha\_0)\right] \tag{8}$$

$$\mathcal{g}(\tau, a\_0) = \mathcal{g}\_N(\tau, a\_0) / \left[ f(\tau) + \Delta \mathcal{g}\_T(\tau, a\_0) \right] \tag{9}$$

where, *A*1 and *A*2 are the slope and curvature of the linear gain at ω<sup>0</sup> , respectively, while 1 *B* and 2 *B* are constants describing changes in *A*1 and *A*2 with saturation, as given in (7) and (8).

The gain spectrum of an SOA is approximated by the following second-order Taylor expansion in Δω:

$$\log(\tau,\alpha) = \mathcal{g}(\tau,\alpha\_0) + \Delta\alpha \frac{\partial \mathcal{g}(\tau,\alpha)}{\partial \alpha} \bigg|\_{a\_b} + \frac{\left(\Delta\alpha\right)^2}{2} \frac{\partial^2 \mathcal{g}(\tau,\alpha)}{\partial \alpha^2} \bigg|\_{a\_b} \tag{10}$$

The coefficients 0 *g*(, ) ω τ ω ω ∂ ∂ and 0 2 2 *g*(, ) ω τ ω ω ∂ ∂ are related to 11 2 *A* ,*B A*, and <sup>2</sup> *B* by (7) and (8).

Here we assumed the same values of 11 2 *A* , *B A*, and <sup>2</sup> *B* as in (Hong *et al.,* 1996) for an AlGaAs/GaAs bulk SOA.

The time derivative terms in (3) have been replaced by the central-difference approximation in order to simulate this equation by the FD-BPM (Das *et al.,* 2000). In simulation, the parameter of bulk SOAs (AlGaAs/GaAs, double heterostructure) with a wavelength of 0.86 μm (Hong *et al.,* 1996) is used and the SOA length is 350 μm. The input pulse shape is sech2 and is Fourier transform-limited. The detail parameters are listed in Table 1 (Section 3).

**Figure 2.** (a) The gain spectra given by the second-order Taylor expansion about the center frequency of the pump pulse 0 ω . The solid line shows the unsaturated gain spectrum (length: 0 μm), the dotted and the dashed-dotted lines are a saturated gain spectrum at 175 μm and 350 μm, respectively. Here, the input pump pulse pulsewidth is 1 ps and pulse energy is 1 pJ. (b) Saturated gain versus the input pump pulse energy characteristics of the SOA. The saturation energy decreases with decreasing the input pump pulsewidth. The SOA length is 350 μm. The input pulsewidths are 0.2 ps, 0.5 ps, and 1 ps respectively, and a pulse energy of 1 pJ.

The gain spectra of SOAs are very important for obtaining the propagation and wave mixing (FWM and optical phase-conjugation between the input pump and probe pulses) characteristics of short optical pulses. Figure 2(a) shows the gain spectra given by a secondorder Taylor expansion about the pump pulse center frequency ω0 with derivatives of *g(*τ*,*  ω*)* by (7) and (8) (Das *et al.,* 2000). In Fig. 2(a), the solid line represents an unsaturated gain spectrum (length: 0 μm), the dotted line represents a saturated gain spectrum at the center position of the SOA (length: 175 μm), and the dashed–dotted line represents a saturated gain spectrum at the output end of the SOA (length: 350 μm), when the pump pulsewidth is 1 ps and input energy is 1 pJ. These gain spectra were calculated using (1), because, the waveforms of optical pulses depend on the propagation distance (i.e., the SOA length). The spectra of these pulses were obtained by Fourier transformation. The "local" gains at the center frequency at *z* = 0, 175, and 350 μm were obtained from the changes in the pulse intensities at the center frequency at around those positions (Das *et al.,* 2001). The gain at the center frequency in the gain spectrum was approximated by the second-order Taylor expression series. As the pulse propagates in the SOA, the pulse intensity increases due to the gain of the SOA. The increase in pulse intensity reduces the gain, and the center frequency of the gain shifts to lower frequencies. The pump frequency is set to near the gain peak, and linear gain 0 *g* is 92 cm-1 at ω<sup>0</sup> . The probe frequency is set -3 THz from for the calculations of FWM characteristics as described below, and the linear gain 0 *g* is -42 cm-1 at this frequency. Although the probe frequency lies outside the gain bandwidth, we selected a detuning of 3 THz in this simulation because the FWM signal must be spectrally separated from the output of the SOA. As will be shown later, even for this large degree of detuning, the FWM signal pulse and the pump pulse spectrally overlap when the pulsewidths become short (<0.5 ps) (Das *et al.,* 2001). The gain bandwidth is about the same as the measured value for an AlGaAs/GaAs bulk SOA (Seki *et al.,* 1981). If an InGaAsP/InP bulk SOA is used we can expect much wider gain bandwidth (Leuthold *et al.,* 2000). With a decrease in the carrier density, the gain decreases and the peak position is shifted to a lower frequency because of the band-filling effect. Figure 2(b) shows the saturated gain versus input pump pulse energy characteristics of the SOA. When the input pump pulsewidth decreases then the small signal gain decreases due to the spectral limit of the gain bandwidth. For the case, when the input pump pulsewidth is short (very narrow, such as 200 fs or lower), the gain saturates at small input pulse energy (Das *et al.,* 2000). This is due to the CH and SHB with the fast response.

90 Optical Communication

expansion in Δω:

The coefficients

the pump pulse 0



Gain, g (c

m-1)

0

50

(a)

100

ω

Probe


Pump (ω0)

Frequency (THz)

respectively, and a pulse energy of 1 pJ.

AlGaAs/GaAs bulk SOA.

00 0 ( , ) ( , )/ ( ) ( , ) *N T g g fg*

and 2 *B* are constants describing changes in *A*1 and *A*2 with saturation, as given in (7) and (8). The gain spectrum of an SOA is approximated by the following second-order Taylor

> (, ) ( ) (, ) (, ) (, ) <sup>2</sup> *g g g g*

∂ ∂ <sup>Δ</sup> = +Δ <sup>+</sup>

0

ω

175 μm 350 μm

<sup>0</sup> Length =

ω

Here we assumed the same values of 11 2 *A* , *B A*, and <sup>2</sup> *B* as in (Hong *et al.,* 1996) for an

The time derivative terms in (3) have been replaced by the central-difference approximation in order to simulate this equation by the FD-BPM (Das *et al.,* 2000). In simulation, the parameter of bulk SOAs (AlGaAs/GaAs, double heterostructure) with a wavelength of 0.86 μm (Hong *et al.,* 1996) is used and the SOA length is 350 μm. The input pulse shape is sech2 and is Fourier transform-limited. The detail parameters are listed in Table 1 (Section 3).

**Figure 2.** (a) The gain spectra given by the second-order Taylor expansion about the center frequency of

The gain spectra of SOAs are very important for obtaining the propagation and wave mixing (FWM and optical phase-conjugation between the input pump and probe pulses)

the dashed-dotted lines are a saturated gain spectrum at 175 μm and 350 μm, respectively. Here, the input pump pulse pulsewidth is 1 ps and pulse energy is 1 pJ. (b) Saturated gain versus the input pump pulse energy characteristics of the SOA. The saturation energy decreases with decreasing the input pump pulsewidth. The SOA length is 350 μm. The input pulsewidths are 0.2 ps, 0.5 ps, and 1 ps

. The solid line shows the unsaturated gain spectrum (length: 0 μm), the dotted and

0

5

10

Saturated Gain (dB)

15

(b)

 ω

2 *g*(, )

τ ω

ω

τ ω  τ

0 2

ω

 τω= + Δ (9)

0 0

2 2

τ ω

ω

ω ω

ω

0.001 0.01 0.1 1 10

Input Energy (pJ)

Pump Pulsewidth = 0.2 ps

0.5 ps 1 ps

<sup>∂</sup> <sup>∂</sup> (10)

are related to 11 2 *A* ,*B A*, and <sup>2</sup> *B* by (7) and (8).

<sup>0</sup> , respectively, while 1 *B*

 τω

τω

where, *A*1 and *A*2 are the slope and curvature of the linear gain at

 τω

2

∂ ∂

τω

*g*(, )

ω

τ ω

∂

0

ω

∂ and

Initially, the MNLSE was used by (Hong *et al.,* 1996) for the analysis of "solitary pulse" propagation in an SOA. We used the same MNLSE for the simulation of FWM and optical phase-conjugation characteristics in SOA using the FD-BPM. Here, we have introduced a complex envelope function V(τ, 0) at the input side of the SOA for taking into account the two (pump and probe) or more (multi-pump or prove) pulses.

#### **2.3. Finite-Difference Beam Propagation Method (FD-BPM)**

To solve a boundary value problem using the finite-differences method, every derivative term appearing in the equation, as well as in the boundary conditions, is replaced by the central differences approximation. Central differences are usually preferred because they

lead to an excellent accuracy (Conte & Boor, 1980). In the modeling, we used the finitedifferences (central differences) to solve the MNLSE for this analysis.

Usually, the fast Fourier transformation beam propagation method (FFT-BPM) (Okamoto, 1992; Brigham, 1988) is used for the analysis of the optical pulse propagation in optical fibers by the successive iterations of the Fourier transformation and the inverse Fourier transformation. In the FFT-BPM, the linear propagation term (GVD term) and phase compensation terms (other than GVD, 1st and 2nd order gain spectrum terms) are separated in the nonlinear Schrödinger equation for the individual consideration of the time and frequency domain for the optical pulse propagation. However, in our model, equation (3) includes the dynamic gain change terms, i.e., the 1st and 2nd order gain spectrum terms which are the last two terms of the right-side in equation (3). Therefore, it is not possible to separate equation (3) into the linear propagation term and phase compensation term and it is quite difficult to calculate equation (3) using the FFT-BPM. For this reason, we used the FD-BPM (Chung & Dagli, 1990; Conte & Boor, 1980; Das *et al.,* 2000; Razagi *et al.,* 2009). If we replace the time derivative terms of equation (3) by the below central-difference approximation, equation (11), and integrate equation (3) with the small propagation step Δz, we obtain the tridiagonal simultaneous matrix equation (12)

$$\frac{\partial}{\partial \tau} V\_k = \frac{V\_{k+1} - V\_{k-1}}{2\Delta \tau}, \quad \frac{\partial^2}{\partial \tau^2} V\_k = \frac{V\_{k+1} - 2V\_k + V\_{k-1}}{\Delta \tau^2} \tag{11}$$

where, ( ) *V V k k* = τ , 1 ( ) *V V k k* τ τ <sup>+</sup> = +Δ , and 1 ( ) *V V k k* τ τ<sup>−</sup> = −Δ

$$\begin{aligned} -a\_k(z + \Delta z) \, V\_{k-1}(z + \Delta z) + \left\{ 1 - b\_k(z + \Delta z) \right\} \, V\_k(z + \Delta z) - c\_k(z + \Delta z) \, V\_{k+1}(z + \Delta z) \\ = a\_k(z) \, V\_{k-1}(z) + \left\{ 1 + b\_k(z) \right\} \, V\_k(z) + c\_k(z) \, V\_{k+1}(z) \end{aligned} \tag{12}$$

where, *k n* = 1, 2, 3, ............, and

$$a\_k(z) = \frac{\Delta z}{2} \left[ \frac{i\beta\_2}{2\Delta\tau^2} + i \frac{1}{4\Delta\tau} \frac{\partial g(\tau, \rho, z)}{\partial \rho \partial} \big|\_{a\_0, \tau\_k} - \frac{1}{4\Delta\tau^2} \frac{\partial^2 g(\tau, \rho, z)}{\partial \rho^2} \big|\_{a\_0, \tau\_k} \right] \tag{13}$$

$$\begin{split} b\_{k}(\mathbf{z}) &= -\frac{\Delta z}{2} \left[ \frac{i\mathcal{B}\_{2}}{\Delta \tau^{2}} + \frac{\mathcal{Y}}{2} + \left( \frac{\mathcal{Y}\_{2p}}{2} + i\mathcal{b}\_{2} \right) \| \, V\_{k}(\mathbf{z}) \|^{2} - \frac{1}{2} \mathcal{g}\_{N}(\tau\_{k}, a\_{0}, \mathbf{z}) (\mathbf{1} + i\alpha\_{N}) \right. \\ &- \frac{1}{2} \Delta \mathcal{g}\_{T}(\tau\_{k}, a\_{0}, \mathbf{z}) (\mathbf{1} + i\alpha\_{T}) - \frac{1}{2\Delta \tau^{2}} \frac{\partial^{2} \mathcal{g}(\tau, a\nu, \mathbf{z})}{\partial a^{2}} \vert\_{a\_{0}, \mathbf{r}\_{k}} \right] \end{split} \tag{14}$$

$$c\_{k}(\mathbf{z}) = \frac{\Delta \mathbf{z}}{2} \left| \frac{i\beta\_{2}}{2\Delta \tau^{2}} - i \frac{1}{4\Delta \tau} \frac{\partial \mathbf{g}(\tau, \rho oz)}{\partial \rho o} \vert\_{a\_{0}, \tau\_{k}} - \frac{1}{4\Delta \tau^{2}} \frac{\partial^{2} \mathbf{g}(\tau, \rho oz)}{\partial \rho o^{2}} \vert\_{a\_{0}, \tau\_{k}} \right| \tag{15}$$

where, Δτ is the sampling time and *n* is the number of sampling. If we know ( ) *V z <sup>k</sup>* , ( *k n* = 1, 2, 3, .........., ) at the position *z* , we can calculate ( ) *Vz z <sup>k</sup>* + Δ at the position of *z z* + Δ

which is the propagation of a step Δ*z* from position *z* , by using equation (12). It is not possible to directly calculate equation (12) because it is necessary to calculate the left-side terms ( ) *<sup>k</sup> az z* + Δ , ( ) *<sup>k</sup> bz z* + Δ , and ( ) *<sup>k</sup> cz z* + Δ of equation (12) from the unknown ( ) *Vz z <sup>k</sup>* + Δ . Therefore, we initially defined ( ) ( ) *k k az z az* +Δ ≡ , ( ) ( ) *k k bz z bz* +Δ ≡ , and ( ) ( ) *k k cz z cz* +Δ ≡ and obtained (0)( ) *Vzz <sup>k</sup>* + Δ , as the zero-th order approximation of ( ) *Vz z <sup>k</sup>* + Δ by using equation (12). We then substituted (0)( ) *Vzz <sup>k</sup>* + Δ in equation (12) and obtained (1)( ) *Vzz <sup>k</sup>* + Δ as the first order approximation of ( ) *Vz z <sup>k</sup>* + Δ and finally obtained the accurate simulation results by the iteration as used in (Brigham, 1988; Chung & Dagli, 1990; Das *et al.,* 2000; Razagi *et al.,* 2009).

92 Optical Communication

where, ( ) *V V k k* =

where, Δ

τ

τ

where, *k n* = 1, 2, 3, ............, and

lead to an excellent accuracy (Conte & Boor, 1980). In the modeling, we used the finite-

Usually, the fast Fourier transformation beam propagation method (FFT-BPM) (Okamoto, 1992; Brigham, 1988) is used for the analysis of the optical pulse propagation in optical fibers by the successive iterations of the Fourier transformation and the inverse Fourier transformation. In the FFT-BPM, the linear propagation term (GVD term) and phase compensation terms (other than GVD, 1st and 2nd order gain spectrum terms) are separated in the nonlinear Schrödinger equation for the individual consideration of the time and frequency domain for the optical pulse propagation. However, in our model, equation (3) includes the dynamic gain change terms, i.e., the 1st and 2nd order gain spectrum terms which are the last two terms of the right-side in equation (3). Therefore, it is not possible to separate equation (3) into the linear propagation term and phase compensation term and it is quite difficult to calculate equation (3) using the FFT-BPM. For this reason, we used the FD-BPM (Chung & Dagli, 1990; Conte & Boor, 1980; Das *et al.,* 2000; Razagi *et al.,* 2009). If we replace the time derivative terms of equation (3) by the below central-difference approximation, equation (11), and integrate equation (3) with the small propagation step Δz,

2

τ

*k k*

<sup>+</sup> = +Δ , and 1 ( ) *V V k k*

{ }

*kk k k kk*

*V V*

 τ

 τ

<sup>2</sup> , <sup>2</sup> *k k k kk*

1 1

− +

( ) ( )1 ( ) ( ) ( ) ( )

*a z zV z z b z z V z z c z zV z z*

1 1 (, ,) (, ,) ( ) <sup>|</sup> <sup>|</sup> 2 4 2 4 *k k <sup>k</sup> <sup>z</sup> <sup>i</sup> gz gz a z <sup>i</sup>*

2

*<sup>z</sup> <sup>i</sup> b z ib V z g z i*

2 0

<sup>1</sup> ( ) | ( )| ( , , )(1 ) 2 22 <sup>2</sup>

*k k Nk N*

=− + + + <sup>−</sup> <sup>+</sup> Δ

0 , 2 2

τ

ω τ

<sup>Δ</sup> ∂ ∂ =− − <sup>Δ</sup> Δ ∂ Δ ∂

τ ω

2 2 2

*p*

1 1 (, ,) ( , , )(1 ) <sup>|</sup> 2 2 *<sup>k</sup>*

> τ ω

> > ω

 α

*g z g zi*

−Δ + − Δ ∂

1 1 (, ,) (, ,) ( ) <sup>|</sup> <sup>|</sup> 2 4 2 4 *k k <sup>k</sup> <sup>z</sup> <sup>i</sup> gz gz c z <sup>i</sup>*

( *k n* = 1, 2, 3, .........., ) at the position *z* , we can calculate ( ) *Vz z <sup>k</sup>* + Δ at the position of *z z* + Δ

 ω

τ

γ γ

τ

*T k T*

Δ

− +Δ +Δ + − +Δ +Δ − +Δ +Δ

*V V V VV*

1 1 1 1 2 2

τ

2 , , 2 2

2

2 , , 2 2

is the sampling time and *n* is the number of sampling. If we know ( ) *V z <sup>k</sup>* ,

τ ω

∂

 ω

> τ

ω τ

Δ ∂ ∂ =+ − <sup>Δ</sup> Δ ∂ Δ ∂

<sup>−</sup> = −Δ

 τ

{ }

*kk k k kk*

1 1

τ ω

 ω

0

τ ω

ω

τω

ω τ

0 0 2

ω τ

> α

ω τ (12)

(13)

(14)

(15)

− +

() () 1 () () () ()

*a zV z b z V z c zV z*

0 0 2

τ

+ − + − ∂ ∂ − − <sup>+</sup> = = ∂ Δ ∂ Δ (11)

 τ

= ++ +

differences (central differences) to solve the MNLSE for this analysis.

we obtain the tridiagonal simultaneous matrix equation (12)

τ

 , 1 ( ) *V V k k* τ

2

τ

β

τ

τω

2

τ

β

β

**Figure 3.** A simple schematic diagram of FD-BPM in the time domain, where, ( /) *<sup>g</sup>* τ = −*t zv* is the local time, which propagates with the group velocity *<sup>g</sup> v* at the center frequency of an optical pulse and Δτ is the sampling time, and *z* is the propagation direction and Δz is the propagation step.

Figure 3 shows a simple schematic diagram of the FD-BPM in time domain. Here, ( /) *<sup>g</sup>* τ = −*t zv* is the local time, which propagates with the group velocity *<sup>g</sup> v* at the center frequency of an optical pulse and Δτ is the sampling time. z is the propagation direction and Δz is the propagation step. With this procedure, we used up to 3-rd time iteration for more accuracy of the simulations.

The FD-BPM (Conte & Boor, 1980; Chung & Dagli, 1990; Das *et al.,* 2000; Razagi *et al.,* 2009a & 2009b) is used for the simulation of several important characteristics, namely, (1) single pulse propagation in SOAs (Das *et al.,* 2008; Razaghi *et al.,* 2009a & 2009b), (2) two input pulses propagating in SOAs (Das *et al.,* 2000; Connelly *et al.,* 2008), (3) Optical DEMUX characteristics of multi-probe or pump input pulses based on FWM in SOAs (Das *et al.,* 2001), (4) Optical phase-conjugation characteristics of picosecond FWM signal in SOAs (Das

*et al.,* 2001), and (5) FWM conversion efficiency with optimum time-delays between the input pump and probe pulses (Das *et al.,* 2007).
