## Abstract

A novel predictor-corrector method for the coupled equations for ultrabroad-band Raman amplifiers with multiple pumps is proposed and derived, for the first time, based on the Adams formula. The proposed algorithm is effective in solving Raman amplifier equations that include pumps, signals, noises, and their backscattering waves. The detail procedure is given, and proves the excellence of our algorithm. Simulation results show that, in designing the Raman amplifier, our multistep method can effectively improve the accuracy and stability compared with the one-step method and explicit multistep method. The numerical results show that the power of backscattering pumps and signals is lower by ~30 dB and 20 dB than their original power, respectively, and the power of forward and backward noises is less than that of input signals by ~30 dB under our simulation conditions.

©2003 Optical Society of America

## 1. Introduction

The advent of the optical amplifier changes the optical fiber communication systems and invites the rapid development of wavelength-division-multiplexed (WDM) transmission systems and networks. Although the entire band of erbium-doped fiber amplifier (EDFA) is fully utilized, the need for more optical channels along installed fibers, increasing channel count and optical bandwidth, pushes EDFA technology beyond its performance limit [1–15]. To increase WDM data capacities and realize the ultra-long-haul system reach, the noise performance of amplifiers becomes critical [8–10, 15–19]. Although the capacities can be upgraded by means of reducing the WDM channel spacing, the nonlinear effects in the transmission fiber such as cross-phase modulation (XPM) and four-wave mixing (FWM) starts to impair the optical waveform significantly [11, 16, 20, 21]. And although these nonlinear effects can be suppressed by reducing the optical signal power launched into the fiber, in turn, the transmission distance is limited due to optical signal-noise-ratio (OSNR) degradation [22–26]. Fortunately, distributed Raman amplifier (DRA) makes the signal-spontaneous beat noise performance significantly improve, longer transmission distances be realized and existing systems be upgraded. What is more, DRA can mitigate fiber nonlinear effects, and improve the OSNR without the need to increase the input optical signal power [12, 19]. Two critical merits of DRA are the low noise and the arbitrary gain band [22–30]. Experiments show that 2.5 Gb/s system could be upgraded to 10 Gb/s by only adding Raman amplifier [24].

A fiber Raman amplifier employs the intrinsic properties of silica fiber to obtain the amplification. Thus, the transmission fiber can be used as the amplification medium, where the gain is created along the transmission [14, 31–33]. The amplification is realized by stimulated Raman scattering (SRS), which occurs when a sufficiently powerful pump is within the same fiber as the signal [23–29]. Due to the random density fluctuations of the optical fiber during the manufacture, small inhomogeneities or microscopic variations in the refractive index exist in the fiber and they can cause Rayleigh scattering [34]. One important issue for DRA is the effect of double-Rayleigh scattering, which makes the maximum DRA gain restricted to around 20 dB [30].

The fiber Raman amplifier is usually designed employing the power analysis [25–28, 35–40], which includes such physical effects as SRS, Rayleigh scattering, spontaneous Raman emission, temperature dependences, arbitrary interaction between the channels, high-order Stokes generation, etc. In the design of the ultra-broad bandwidth Raman amplifier with multiple pumps, however, it takes exhaustive computing time to achieve well-behaved results by using direct integration of coupled differential equations [39, 40]. The required simulation time may even make the design be impractical if the bandwidth is large and the transmission fiber length is long. Fortunately, many practical methods are proposed [36, 39, 40]. For example, Ref. [36] can predict the gain tilt for *N* channels with the equal power and spacing, based on the so-called triangular approximation. Min et al., obtained an efficient method with the average power analysis, which reduces the computing time of over two orders of magnitude compared with the direct integration approach based on ordinary coupled differential equations [40]. Although one of us has reported a multistep method based on the four-step Adams-Bashforth method in [39], which can reduce the computing time over three orders of magnitude with excellent accuracy promises in comparison with the direct integration method for the Raman amplifier propagation equations, it is not perfect and its stability is limited.

In this paper, a predictor-corrector method (PCM) based on the Adams-Bashforth-Moulton formula for the Raman amplifier propagation equations is proposed for the first time to the authors’ knowledge. This proposed method is compared with the four-step method in [39] and the one-step method in [40]. The results show that it is a fast and stable numerical method, which can be applied to real system performance evaluations. The numerical results also exhibit that the power of backscattering pumps and signals is lower than their original power by ~30 dB and 20 dB, respectively, under our simulation conditions.

## 2. Mathematical model

In the DRA WDM systems with the bi-directional pump/signal flow, the analysis of the bi-directional signal propagation in fibers is essential. The noises in this system include amplified spontaneous emission (ASE), thermal noise, double Rayleigh backscattering, cross gain saturation noise, nonlinear impairment, pump-to-pump SRS and signal-to-signal SRS, and so on [1–42]. In the Raman amplifiers, SRS can lead to the energy exchange between forward- and backward-propagating waves. The spontaneous Raman scattering of unidirectional signals occurs in both directions with a frequency shift, and the scattered waves are further amplified by SRS. Therefore, the linear and nonlinear effects couple signals/pumps propagating in the same and opposite directions, and to solve for their propagation becomes a two-point boundary problem. For the DRA WDM systems, the iterative method is usually used to solve such a problem [41].

The analysis of DRA is based on a set of coupled steady-state equations that include spontaneous Raman emission and its temperature dependence, Rayleigh scattering including multiple reflections, ASE, SRS, high-order Stokes generation, and arbitrary interaction between an unlimited number of pumps and signals. The forward and backward power evolution of pumps, signals and ASE can be expressed in terms of the following equations [25–29]

$$\pm {P}^{\pm}(z,{v}_{i})\sum _{m=1}^{i-1}\frac{{g}_{R}\left({v}_{m}-{v}_{i}\right)}{\Gamma {A}_{\mathit{eff}}}\left[{P}^{\pm}(z,{v}_{i})+{P}^{\mp}(z,{v}_{i})\right]$$

$$\pm h{v}_{i}\sum _{m=1}^{i-1}\frac{{g}_{R}\left({v}_{m}-{v}_{i}\right)}{\Gamma {A}_{\mathit{eff}}}\left[{P}^{\pm}(z,{v}_{i})+{P}^{\mp}(z,{v}_{i})\right]\left[1+{\left({e}^{\frac{h\left({v}_{m}-{v}_{i}\right)}{kT}}-1\right)}^{-1}\right]\Delta v$$

$$\mp {P}^{\pm}(z,{v}_{i})\sum _{m=i+1}^{n}\frac{{v}_{i}}{{v}_{m}}\frac{{g}_{R}\left({v}_{i}-{v}_{m}\right)}{\Gamma {A}_{\mathit{eff}}}\left[{P}^{\pm}(z,{v}_{i})+{P}^{\mp}(z,{v}_{i})\right]$$

$$\mp 2h{v}_{i}{P}^{\pm}(z,{v}_{i})\sum _{m=i+1}^{n}\frac{{v}_{i}}{{v}_{m}}\frac{{g}_{R}\left({v}_{i}-{v}_{m}\right)}{\Gamma {A}_{\mathit{eff}}}\left[1+{\left({e}^{\frac{h\left({v}_{i}-{v}_{m}\right)}{kT}}-1\right)}^{-1}\right]\Delta \mu $$

where *P*
^{+}(*z, v _{i}*) and

*P*

^{-}(

*z, v*) are optical power of forward- and backward-propagating waves within infinitesimal bandwidth around

_{i}*v*, respectively;

_{i}*α*,

*η*,

*h*,

*k*and

*T*are attenuation coefficient, Rayleigh-backscattering coefficient, Planck’s constant, Boltzmann constant and temperature, respectively;

*A*is effective area of optical fiber at frequency

_{eff}*v*;

_{m}*g*(

_{R}*v*) is Raman gain parameter at frequency

_{i}-v_{m}*v*due to pump at frequency

_{i}*v*[41]; the factor of Γ accounts for polarization randomization effects, whose value lies between 1 and 2. Because spontaneous emission and thermal noise are not correlated with signal, a factor of two in the noise terms [i.e., the sixth terms in Eq. (1)] is put. The frequency ratio

_{m}*v*/

_{i}*v*describes vibrational losses. Terms from

_{m}*m*=1 to

*m*=

*i*-1 and from

*m*=

*i*+1 to

*m*=

*n*cause amplification and attenuation of the channel at frequency

*v*, respectively. Δ

_{i}*v*and Δ

*µ*are the spectral noise interval for the noise increase and loss (in this paper, they are assumed to be equal; see Fig. 1) in the calculation. The first two terms in the right hand of Eq. (1) denote the fiber loss and Rayleigh back-scattering, the third term does the Raman gain due to shorter wavelength, the forth term does the ASE noise with thermal factor, the fifth term does the pump depletion due to longer wavelength, and the sixth term does the loss due to noise emission.

## 3. Numerical method

#### 3.1. Predictor-corrector method

For the mathematical convenience, Eq. (1) is rewritten as

$$\mp \sum _{m=i+1}^{n}\frac{{v}_{i}}{{v}_{m}}\frac{{g}_{R}\left({v}_{i}-{v}_{m}\right)}{\Gamma {A}_{\mathit{eff}}}\left[{P}^{\pm}(z,{v}_{i})+{P}^{\mp}(z,{v}_{i})+2h{v}_{i}\left(1+{\left({e}^{\frac{h\left({v}_{i}-{v}_{m}\right)}{kT}}-1\right)}^{-1}\right)\Delta \mu \right]$$

$$\pm \sum _{m=1}^{i-1}\frac{{g}_{R}\left({v}_{m}-{v}_{i}\right)}{\Gamma {A}_{\mathit{eff}}}\left[{P}^{\pm}(z,{v}_{i})+{P}^{\mp}(z,{v}_{i})\right]\left[1+\frac{h{v}_{i}}{{P}^{\pm}(z,{v}_{i})}\left(1+{\left({e}^{\frac{h\left({v}_{m}-{v}_{i}\right)}{kT}}-1\right)}^{-1}\right)\Delta v\right]$$

Because the light amplification or attenuation resembles the exponential pattern with the propagation distance, the exponential change of variable is advantageous in amplifier simulations [43]. In Ref. [40], the average power analysis technique was applied. This technique is one-step method, which involves the information from only one of the previous mesh points *z _{j}*, although it could greatly decrease the computing time. According to this method, Eq. (2) can be solved as

where Δ*z* is the step size for each elemental amplifier section.

Although this method might use functional evaluation information at points between *z _{j}* and

*z*

_{j+1}, they do not retain that information for direct use in future approximations. Since the approximate solution is available at each of the mesh points z

_{0}, z

_{1}… z

*before the approximation at z*

_{j}_{j+1}is obtained, and because the multistep method can improve the accuracy of solution at each step [44–47], it can effectively decrease the computing time under the same accuracy [39]. The Adams methods are the most elementary multistep methods, since they use interpolation on the immediate previous approximation. The traditional Adams-Bashforth formula is an explicit formula with the order three (a four-step continuation formula). In the previous report, one of us proposed a multistep method based on the four-step Adams-Bashforth formula [39], which greatly shortened the computing time.

However, the implicit Adams-Moulton method gives considerably better results than the explicit Adams-Bashforth method of the same order [44–47]. Although this is generally the case, the implicit methods have the inherent weakness of first having to convert the method algebraically to an explicit representation. Therefore, the combination of an explicit and implicit technique is usually used, which is called a PCM. The explicit method (i.e., four-step Adams-Bashforth method) predicts an approximation, and the implicit method (i.e., three-step Adams-Moulton method) corrects this prediction.

Firstly, one can use the four-step method based on the Adams-Bashforth formula as predictor. Then, Eq. (2) can be solved as [39]

Secondly, this approximating is improved by use of the three-step method based on the Adams-Moulton formula as corrector, i.e.,

where
$\overline{F({z}_{j+1},v)}$
is calculated from Eqs. (4) and (2b). *F*(z* _{j}*,

*v*),

*F*(z

_{j-1},

*v*),

*F*(z

_{j-2},

*v*) and

*F*(z

_{j-3},

*v*), which correspond the previous mesh points z

*, z*

_{j}_{j-1}, z

_{j-2}and z

_{j-3}, are obtained from Eq.(2b) after suitable substitutions.

Comparing Eqs. (4), (5) with (3), we see that our multistep method obtains the solution of point z_{j+1} from the previous mesh points z_{j}, z_{j-1}, z_{j-2} and z_{j-3}, instead of only one point z_{j} as in Eq. (3). For the starting values of *F*(z_{1}, *v*), *F*(z_{2}, *v*) and *F*(z_{3}, *v*), one can obtain them from one-step, two-step and three-step methods based on the Adams-Bashforth-Moulton formula, respectively. They are

for the one-step method, and

for the two-step method, and

for the three-step method.
$\overline{F({z}_{2},v)}$
and
$\overline{F({z}_{3},v)}$
in Eqs. (7b) and (8b) are from
$\overline{{P}^{\pm}({z}_{2},v)}$
and
$\overline{{P}^{\pm}({z}_{3},v)}$
in Eqs. (7a) and (8a), respectively. Then, *F*(z_{0}, *v*), *F*(z_{1}, *v*), *F*(z_{2}, *v*) and *F*(z_{3}, *v*) can be obtained after substituting *P*(z_{0}, *v*), *P*(z_{1}, *v*), *P*(z_{2}, *v*) and *P*(z_{3}, *v*) into Eq.(2b), where *P*(z_{0}, *v*) is the initial value.

#### 3.2. Simulation algorithm

Equations (1) or (2) describing the signal propagation in the forward direction contains amplitudes of the backward-propagating signals and vice versa. In the general calculation, firstly, the reasonable value for *P*
^{-}(*z*) along the fiber length is assumed. Then, Eq. (1) is solved in the forward direction, using boundary conditions for *P*
^{+}(*z*) at *z*=0 and the initial guess for *P*
^{-}(*z*). The computed values give an approximate solution, which can be used to integrate equations for the backward-propagating waves, with boundary conditions for *P*
^{-}(*z*) at *z*=*L*, in the backward direction. A new and improved approximation for *P*
^{-}(*z*) is given, which can be used in the next iterations. The signal values at both fiber ends *P*
^{-}(0) and *P*
^{+}(*L*) obtained in the successive iterations are checked for the convergence.

However, the simulation results show that, for the backward multiple pumps with large power (Note that: this value is related to the fiber length, the number of backward-propagating pumps, the demand of accuracy, and so on), the iterations may not converge if the boundary conditions are defined at *z*=0 according to the algorithm explained in the previous paragraph. Fortunately, the iterations converge under the common conditions if one defines the boundary conditions for the signal channels at *z*=*L* instead of *z*=0 and skip the noise effects [38, 41]. The processing detail is similar to the previous paragraph algorithm except the starting boundary conditions at *z*=*L* instead of *z*=0. Therefore, we divide the computational method for multiple pumps DRA into two steps. In the first step, the noises are neglected and the shooting method is employed [41]. Then, the boundary conditions for the backward-propagating pumps at *z*=0 can be solved. In the second step, with the results of the first step and reasonable initial values ^{1)}, one can calculate Eq. (2) including the noise effects by use of our proposed multistep method. The simulation results show that this two-step algorithm can be convergent.

#### 3.3. Test and results

To compare the accuracy of our multistep method with that of one-step method in [40] and of four-step method in [39], an ordinary differential equation is tested, i.e.,

Substituting y(*z*) and (-*y*(*z*)+*z*+1)/*y*(*z*) of Eq. (9) into *P*
^{±}(*z*, *v*) and *F*(*z*, *v*) of Eqs.(2a) and (2b), we can obtain Fig. 2
^{2}). Lines with square, triangular and circle symbols in Fig. 2 are the absolute errors ^{3)}, which are obtained from the one-step method, the four-step method based on Adams-Bashforth formula in [39], and the PCM in this paper, respectively. Details of the procedure are shown in the Appendix, where such symbols as *y*, *z*, *F*, *dz* are same with Eq.9. Fig. 2 shows that, under the same step size Δ*z*, the absolute error for the PCM is lower about three and one order of magnitude compared with the one-step method [40] and the four-step method [39]. Additionally, the simulation results also imply that, under the same accuracy, the step size Δ*z* of four-step method based on the Adams-Bashforth formula can be more 20 times than that of one-step method ^{4)} [39].

To compare the stability between the one-step method in [40], the four-step method in [39], and the PCM in this paper, we employ the same test function of Eq. (4) in [39], which is

The exact solution of Eq. (10) is

Lines with none, square, circle and triangular symbols in Fig. 3 are obtained from the exact solution, one-step method, four-step Adams-Bashforth method in [39], and PCM in this paper, respectively. Figure 3 shows that, at the same step size of Δ*z*=0.1, ① the one-step and four-step Adams-Bashforth methods do diverge when *z*>1.2 and *z*>1.5, respectively; ② the PCM is stable during the entire regime.

Although the PCM costs the computing time as about two times as the four-step Adams-Bashforth method [39], its accuracy (see Fig. 2) and stability (see Fig. 3) is evident, and then the step size can be lengthened. Hence, the computing time in the PCM can be reduced with a little sacrifice of the accuracy. In fact, numerical results show that the PCM is more accurate and stable than the four-step Adams-Bashforth method under the same computing time in calculating Eq. (2). This method can be applied to real system performance evaluations.

## 4. Numerical simulation and results

In the numerical calculations, we assume that: *L*=80 km, *T*=300, *η*=7×10^{-8} m^{-1}, *A _{eff}*=80×10

^{-12}m

^{2},

*α*=0.20 and 0.25 dB/km for signals and pumps, respectively. The amplifier is designed to have a total input power of 50 mW for 100 signal channels spaced Δ

*v*=100 GHz/channel, from 185 THz to 195 THz. The pumps are chosen to realize an 80-nm bandwidth with the gain variation of less than 1 dB, and their powers are 20.6, 20.8, 21.0, 19.4, 18.6, 18.1, 16.9 and 21.9 dBm corresponding to wavelengths of 1432, 1442, 1452, 1462, 1472, 1485, 1499 and 1510 nm, respectively. In this work, we select

*v*

_{0}=211.95 THz and

*v*=182.05 THz, having 300 channels for noises in all (see Fig. 1). From Eq. (2) and Fig. 1, it is easily seen that the number of equations for this numerical simulation is 2×(300+100+8)

_{n}^{5)}.

Figures 4 and 5 demonstrate the use of Raman gain to compensate attenuation within the transmission fiber, the values of which come from the numerical simulation for Eq. (2) (with 816 equations). The PCM is employed with the fiber parts of 29, and the relative error is less than 10^{-3}.

Figure 4 shows spectra of the fiber input, output and backscattering, where (a) for input pump (backward) and signal (forward); (b) for output pump, backscattering signal and backward noise at the input-end of the fiber (i.e., *z*=0) due to SRS and single-Rayleigh scattering; (c) for output signal, forward noise and backscattering pump at the output-end of the fiber (i.e., *z*=*L*) due to SRS and double-Rayleigh scattering. The inset figure in Fig. 4(c) shows that the WDM channels obtain nearly equal gain after the careful arrangements of pumps. By comparison Fig. 4(a) with Fig. 4(b), the ‘bluest’ pumps suffer the greatest depletion and the ‘redder’ pump can cause the depletion of the ‘bluer’ pump. Thus, the power and wavelength of pump should be suitably designed because of strong interaction between them. It can be seen, from Figs. 4(a)–(c), that SRS, the single and double Rayleigh scattering can cause noise and idler waves. Therefore, the filter and isolator are needed for the systems and networks in order to improve OSNR [37].

Figure 5 exhibits the internal powers of signals, pumps, noises and their backscattering waves within the fiber, where the backwards pumps deplete towards the input-end of the fiber [see Fig. 5(a)].In Fig. 5(a), the dashes in the plane of the pump power with the wavelength are the projection of eight curves. From the plot of the pump power along the transmission fiber (backward direction) in Fig. 5(a), we see that the interaction of different pump waves is significant, as causes amplification of some pumps before they are depleted by the signal. Figure 5(c) shows that signals decrease when *z*<55 km, then increase up to 80 km. In the former, the fiber attenuation dominates, but the Raman gain dominates continuously in the latter. Both Figs. 4 and 5 show that: ① there is the strong interaction between the signals due to SRS; ② the Rayleigh scattering creates intense backscattering pumps and signals; ③ backward and forward noises are produced by the SRS, single- and double-Rayleigh scattering; ④ the power of backscattering pumps and signals is lower by about 30 dB and 20 dB than that of pumps and signals; ⑤ the power of forward and backward noises is less ~30 dB than that of input signals. By the suitable design, the bandwidth of >80 nm and the fluctuation of < 1 dB can be implemented for the length of >80 km by the Raman amplifiers.

## 5. Discussion

In Section 3.3, the numerical calculations show that our proposed PCM is superior to one-step and four-step method in the accuracy and stability. To further comprehend the superiority of PCM, we compare the convergence of this method with that of four-step method in [39]. The results are demonstrated in Fig. 6, where (a) for PCM and (b) for four-step method. All parameters in this figure are same with those in Fig. 5 except that the single-pump power is 1.3 W at the wavelength of 1470 nm and there are only fifteen signals, from 192.15 THz to 194.95 THz. The traditional iteration algorithm based on the shooting method is employed for the two-point boundary condition problem in the numerical calculation. It can be easily seen that PCM is convergent [Fig. 6(a)], but that four-step method is divergent [Fig. 6(b)] under above simulation conditions. Therefore, the proposed PCM is more stable than the reported four-step method in [39].

In [41], we proposed another effective shooting algorithm based on the Runge-Kutta formula, which can be availably employed into fiber amplifiers. However, when the pump power exceeds 1.1 W in Fig. 6, the method in [41] is divergent. By comparing the proposed method in this paper with that in [41], it is easily seen that the former is excellent to the latter.

To compare the starting boundary conditions at *z*=*L* with those at *z*=0, the detail calculation procedure of our proposed algorithm is given in Fig. 7, where shows the evolution of all signals along the fiber in each iteration. In the simulations, all parameters are same with Fig.5 and the noises are neglected (correspond to the first step of two-step procedure for the simulation algorithm in Section 3.2). Figures 7(a), (b), (c), and (d) are the first-, second-, third-, and fourth-iteration for our algorithm, and their corresponding maximum of relative error of all signals are 0.3388, 0.0591, 0.0019, 9.844×10^{-4}, respectively. The initially assumed value for each signal at *z*=*L* (i.e., 80 km) is 0.15 mW (i.e., -8.239 dBm), which are shown in Fig. 7(a). After only 4 iterations, the maximum of relative error of all signals reaches the required value of <10^{-3} [see Fig. 7(d)]. However, if it is calculated from *z*=0 to *L* (correspond to the starting boundary conditions at *z*=0), there need 29 iterations to obtain the same error (the detail procedure is not given here, due to too many iterations). Therefore, under the same conditions, the procedure of calculating from *z*=*L* to 0 can great exceed to that of calculating from *z*=0 to *L*. The comparisons of two cases show the asymmetric convergence speed of our proposed approach on Raman amplifier problems.

Although we employ the four-step PCM to simulate the Raman amplifier propagation equations, there are other order PCM, which can be derived from other order Adams-Bashforth-Moulton formula similarly. For example, Eqs. (7a)–(7b) and (8a)–(8b) are two-step and three-step PCM, respectively. The five-step PCM is as

$$\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}+2616\xb7F({z}_{j-2},v)-1274\xb7F({z}_{j-3},v)+251\xb7F({z}_{j-4},v)\left)\frac{\Delta z}{720}\right],$$

$$\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}-264\xb7F({z}_{j-1},v)+106\xb7F({z}_{j-2},v)-19\xb7F({z}_{j-3},v))\frac{\Delta z}{720}].$$

Theoretically, the accuracy can increase with the order of multistep. On the other hand, the more the order of multistep is, the more the computing time costs under the same number of steps. In general, four-step Adams-Bashforth-Moulton method can availably improve the accuracy and decrease the computing time with larger step sizes [44–48]. In fact, we find out that four-step PCM is better than other order PCM at given consideration to the accuracy and computing time.

## 6. Conclusions

The theoretical model for ultrabroad bandwidth Raman amplifiers with the multi-wavelength and different power level pumping is presented in detail. A fast and accurate numerical method for this model is proposed and derived, for the first time to the authors’ knowledge, based on the four-step Adams-Bashforth and three-step Adams-Moulton method. The proposed effective algorithm for the two-point boundary problem can numerically solve the Raman amplifier equations, which includes pumps, signals, noises, and their backscattering waves. The numerical simulation shows that only 4 iterations can reach the required relative error when calculating from *z*=*L* to 0, but it need 29 iterations to obtain the same error under the same conditions when calculating from *z*=0 to *L*. Test results also demonstrate that our multistep method can effectively improve the accuracy and stability in designing DRA. Compared with the direct integration method for the ordinary coupled differential equations, our proposed method can decrease the computing time over three orders of magnitude with excellent accuracy promises [39, 40]. Additionally, this PCM can improve the accuracy and stability in comparison with the four-step method based on Adams-Bashforth formula in the previous report [39]. The numerical results show that: ① the strong interactions between pump-to-pump, pump-to-signal and signal-to-signal are caused by SRS and Rayleigh scattering; ② SRS can not only increase the signal power but also create the noise; ③ the power of backscattering pumps and signals is lower by ~30 dB and 20 dB than their original power, respectively; ④ the power of forward and backward noises is less than that of input signals by ~30 dB. The results obtained have great significance for designing the ultrabroad-band Raman amplifier in the WDM systems and networks, and this technique can be applied to real system performance evaluations.

## Appendix

function CalFig2

dz=0.1; z=0:dz:4; % dz is equivalent to Δz

y0=ExaSol(z); % exact solution

y(1,:)=y0(1:4); y=repmat(y,[3,1]); % initial values

F(1,:)=dydzy(z(1:4),y(1,:)); F=repmat(F,[3,1]);

for i=4:(length(z)-1),

y(1,i+1)=y(1,i)×exp(F(1,i)×dz);

F(1,i+1)=dydzy(z(i+1),y(1,i+1)); % for one-step method

y(2,i+1)=y(2,i)×exp(dz/24×(55×F(2,i)-59×F(2,i-1)+37×F(2,i-2)-9×F(2,i-3)));

F(2,i+1)=dydzy(z(i+1),y(2,i+1)); % for 4-step method

y(3,i+1)=y(3,i)×exp(dz/24×(55×F(3,i)-59×F(3,i-1)+37×F(3,i-2)-9×F(3,i-3)));

F(3,i+1)=dydzy(z(i+1),y(3,i+1)); % for predictor

y(3,i+1)=y(3,i)×exp(dz/24×(9×F(3,i+1)+19×F(3,i)-5×F(3,i-1)+F(3,i-2)));

F(3,i+1)=dydzy(z(i+1),y(3,i+1)); % for corrector

end

semilogy(z,abs((y0-y(1,:))./y0),z,abs((y0-y(2,:))./y0),‘.:b’,z,abs((y0-y(3,:))./y0),‘-.y˄’);

legend(‘one-step’,‘4-step method’,‘PCM’);

function fy=dydzy(z,y) % y’/y=f(z,y)/y

fy=(-y+z+1)./y;

function y1=ExaSol(z) % exact solution

y1=exp(-z)+z;

## Acknowledgement

This work was supported by the Ministry of Education of Korea through the Brain Korea 21 (Information Technology) Program.

## Footnotes

^{1)} | The initial values for the power of backscattering pumps and forward noises at z=0 and of backscattering signals and backward noises at z=L are -70 dBm in this calculations (see Fig. 5). |

^{2)} | In the calculation, the starting four values for F(z
_{0, 1, 2, 3}) are from the exact expression. |

^{3)} | The exact solution is that y(z)=exp(-z)+z. |

^{4)} | See figure 1 and its descriptions in [39]. |

^{5)} | The number of coupled equations is 300 for noises, 100 for signals and 8 for pumps, and the same number for their backscattering waves. See Figs. 4 and 5. |

## References and links

**1. **M. N. Islam, “Raman amplifiers for telecommunications,” IEEE J. Sel. Top. Quantum Electron. **8**, 548–559 (2002). [CrossRef]

**2. **E. M. Dianov, “Advances in Raman fibers,” J. Lightwave Technol. **20**, 1457–1462 (2002). [CrossRef]

**3. **M. Karásek and M. Menif, “Channel addition/removal response in Raman fiber amplifiers: modeling and experimentation,” J. Lightwave Technol. **20**, 1680–1687 (2002). [CrossRef]

**4. **A. Carena, V. Curri, and P. Poggiolini, “On the optimization of hybrid Raman/erbium-doped fiber amplifiers,” IEEE Photon. Technol. Lett. **13**, 1170–1172 (2001). [CrossRef]

**5. **N. Kikuchi, K. K. Wong, K. Uesaka, K. Shimizu, S. Yam, E. S. Hu, M. Marhic, and L. G. Kazovsky, “Novel in-service wavelength-band upgrade scheme for fiber Raman amplifier,” IEEE Photon. Technol. Lett. **15**, 27–29 (2003). [CrossRef]

**6. **K. Suto, T. Saito, T. Kimura, J. I. Nishizawa, and T. Tanabe, “Semiconductor Raman amplifier for terahertz bandwidth optical communication,” J. Lightwave Technol. **20**, 705–711 (2002). [CrossRef]

**7. **H. S. Seo, K. Oh, and U. C. Paek, “Gain optimization of germanosilicate fiber Raman amplifier and its applications in the compensation of Raman-induced crosstalk among wavelength division multiplexing channels,” IEEE J. Quantum Electron. **37**, 1110–1116 (2001). [CrossRef]

**8. **B. Cuenot, “Comparison of engineering scenarios for N×160 Gb/s WDM transmission systems,” IEEE Photon. Technol. Lett. **15**, 864–866 (2003). [CrossRef]

**9. **A. Pizzinat, M. Santagiustina, and C. Schivo, “Impact of hybrid EDFA-distributed Raman amplification on a 4×40-Gb/s WDM optical communication system,” IEEE Photon. Technol. Lett. **15**, 341–343 (2003). [CrossRef]

**10. **V. E. Perlin and H. G. Winful, “Optimizing the noise performance of broad-band WDM systems with distributed Raman amplification,” IEEE Photon. Technol. Lett. **14**, 1199–1201 (2002). [CrossRef]

**11. **E. Poutrina and G. P. Agrawal, “Timing jitter in dispersion-managed soliton systems with distributed, lumped, and hybrid amplification,” J. Lightwave Technol. **20**, 762–769 (2002). [CrossRef]

**12. **D. Dahan and G. Eisenstein, “Numerical comparison between distributed and discrete amplification in a point-to-point 40-gb/s 40-WDM-based transmission system with three different modulation formats,” J. Lightwave Technol. **20**, 379–388 (2002). [CrossRef]

**13. **S. Radic, S. Chandrasekhar, P. Bernasconi, J. Centanni, C. Abraham, N. Copner, and K. Tan, “Feasibility of hybrid Raman/EDFA amplification in bidirectional optical transmission,” IEEE Photon. Technol. Lett. **14**, 221–223 (2002). [CrossRef]

**14. **A. G. Okhrimchuk, G. Onishchukov, and E. Lederer, “Long-haul soliton transmission at 1.3 µm using distributed Raman amplification,” J. Lightwave Technol. **19**, 837–841 (2001). [CrossRef]

**15. **N. Takachio and H. Suzuki, “Application of Raman-distributed amplification to WDM transmission systems using 1.55-µm dispersion-shifted fiber,” J. Lightwave Technol. **19**, 60–69 (2001). [CrossRef]

**16. **X. zhou, M. Birk, and S. Woodward, “Pump-noise induced FWM effect and its reduction in a distributed Raman fiber amplifier,” IEEE Photon. Technol. Lett. **14**,1686–1688 (2002). [CrossRef]

**17. **B. Min, P. Kim, and N. Park, “Flat amplitude equal spacing 798-channel Rayleigh-assisted Brillouin/Raman multiwavelength comb generation in dispersion compensating fiber,” IEEE Photon. Technol. Lett. **13**, 1352–1354 (2001). [CrossRef]

**18. **T. Okuno, T. Tsuzaki, and M. Nishimura, “Novel optical hybrid line configuration for quasi-lossless transmission by distributed Raman amplification,” IEEE Photon. Technol. Lett. **13**, 806–808 (2001). [CrossRef]

**19. **L. D. Garrett, M. Eiselt, R. W. Tkach, V. Dominic, R. Waarts, D. Giltner, and D. Mehuys, “Field demonstration of distributed Raman amplification with 3.8-dB Q-improvement for 5×120-km transmission,” IEEE Photon. Technol. Lett. **13**, 157–159 (2001). [CrossRef]

**20. **W. S. Wong, C. J. Chen, M. C. Ho, and H. K. Lee, “Phase-matched four-wave mixing between pumps and signals in a copumped Raman amplifier,” IEEE Photon. Technol. Lett. **15**, 209–211 (2003). [CrossRef]

**21. **K. Song and S. D. Dods, “Cross modulation of pump-signals in distributed Raman amplifiers, theory and experiment,” IEEE Photon. Technol. Lett. **13**, 1173–1175 (2001). [CrossRef]

**22. **H. Kim and R. J. Essiambre, “Transmission of 8×20 Gb/s DQPSK signals over 310-km SMF with 0.8-b/s/Hz spectral efficiency,” IEEE Photon. Technol. Lett. **15**, 769–771 (2003). [CrossRef]

**23. **K. Toge, K. Hogari, and T. Horiguchi, “Measurement of Raman gain distribution in optical fibers,” IEEE Photon. Technol. Lett. **14**, 974–976 (2002). [CrossRef]

**24. **P. B. Hansen, L. Eskildsen, S. G. Grubb, A. J. Stentz, T. A. Strasser, J. Judkins, J. J. DeMarco, R. Pedrazzani, and D. J. DiGiovanni, “Capacity upgrades of transmission systems by Raman amplification,” IEEE Photon. Technol. Lett. **9**, 262–264 (1997). [CrossRef]

**25. **H. Kidorf, K. Rottwitt, M. Nissov, M. Ma, and E. Rabarijaona, “Pump interactions in a 100-nm bandwidth Raman amplifier,” IEEE Photon. Technol. Lett. **11**, 530–532 (1999). [CrossRef]

**26. **M. Achtenhagen, G. G. Change, B. Nyman, and A. Hardy, “Analysis of a multiple-pump Raman amplifier,” Appl. Phys. Lett , **78**, 1322–1324 (2001). [CrossRef]

**27. **X. zhou, C. Lu, P. Shum, and T. H. Cheng, “A simplified model and optimal design of a multiwavelength backward-pumped Raman amplifier,” IEEE Photon. Technol. Lett. **13**, 945–947 (2001). [CrossRef]

**28. **P. C Xiao, Q. J zeng, J. Huang, and J. M. Liu, “A new optimal algorithm for multipump sources of distributed fiber Raman amplifier,” IEEE Photon. Technol. Lett. **15**, 206–208 (2003). [CrossRef]

**29. **L. Helczynski and A. Berntson, “Comparison of EDFA and bidirectionally pumped Raman amplifier in a 40-Gb/s Rz transmission system,” IEEE Photon. Technol. Lett. **13**, 669–671 (2001). [CrossRef]

**30. **P. M. Krummrich, R. E. Neuhauser, and C. Glingener, “Bandwidth limitations of broadband distributed Raman fiber amplifiers for WDM systems,” in Optical Fiber Communications Conference 2001, MI3-1, 2001.

**31. **Z. M. Liao and G. P. Agrawal, “Role of distributed amplification in designing high-capacity soliton systems.” Opt. Express **9**, 66–71 (2001), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-9-2-66 [CrossRef] [PubMed]

**32. **C. Finot, G. Millot, C. Billet, and J. M. Dudley, “Experimental generation of parabolic pulses via Raman amplification in optical fiber,” Opt. Express **11**, 1547–1552 (2003), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-11-13-1547 [CrossRef] [PubMed]

**33. **T. E. Murphy, “10-GHz 1.3-ps pulse generation using chirped soliton compression in a Raman gain medium,” IEEE Photon. Technol. Lett. **14**, 1424–1426 (2002). [CrossRef]

**34. **M. Nakazawa, “Rayleigh backscattering theory for single-mode optical fibers,” J. Opt. Soc. Amer. **73**, 1175–1181 (1983). [CrossRef]

**35. **P. B. Hansen, L. Eskildsen, A. J. Stentz, T. A. Strasser, J. Judkins, J. J. DeMarco, R. Pedrazzani, and D. J. DiGiovanni, “Rayleigh scattering limitations in distributed Raman pre-amplifiers,” IEEE Photon. Technol. Lett. **10**, 159–161 (1998). [CrossRef]

**36. **D. N. Christodoulides and R. B. Jander, “Evolution of stimulated Raman crosstalk in wavelength division multiplexed systems,” IEEE Photon. Technol. Lett. **8**, 1722–1724 (1996). [CrossRef]

**37. **C. M. McIntosh, A. G. Grandpierre, D. N. Christodoulides, J. Toulouse, and J. M. P. Delavaux, “Eliminating SRS channel depletion in massive WDM systems via optical filtering techniques,” IEEE Photon. Technol. Lett. **13**, 302–304 (2001). [CrossRef]

**38. **V. E. Perlin and H. G. Winful, “Optimal design of flat-gain wide-band fiber Raman amplifiers,” J. Lightwave Technol. **20**, 250–254 (2002). [CrossRef]

**39. **X. M Liu, H. Y zhang, and Y. L Guo, “A novel method for Raman amplifier propagation equations,” IEEE Photon. Technol. Lett. **15**, 392–394 (2003). [CrossRef]

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

**41. **X. M Liu and B Lee, “Effective shooting algorithm and its application to fiber amplifiers,” Opt. Express **11**, 1452–1461 (2003), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-11-12-1452 [CrossRef] [PubMed]

**42. **P. Kim, J. Park, H. Yoon, J. Park, and N. Park, “In situ design method for multichannel gain of a distributed Raman amplifier with multiwave OTDR,” IEEE Photon. Technol. Lett. **14**, 1683–1685 (2002). [CrossRef]

**43. **M. E. Marhic and D. E. Nikonov, “Low third-order glass-host nonlinearities in erbium-doped waveguide amplifiers,” Proceedings of SPIE , vol. **4645**, pp. 193 (2002). [CrossRef]

**44. **A. Quarteroni, R. Sacco, and F. Saleri, *Numerical Mathematics* (Springer-Verlag, New York, 2000).

**45. **B. F. Plybon, *An Introduction to Applied Numerical Analysis* (PWS-KENT Publishing Company, Boston, 1992), pp. 428–441.

**46. **J. H. Mathews, *Numerical Methods for Mathematics, Science, and Engineering* (Second edition, Prentice Hall, New Jersey, 1992), pp. 464–475.

**47. **J. D. Faires and R. L. Burden, *Numerical method*. Boston (PWS-KENT Publishing Company, Boston, 1992), pp. 168–179.

**48. **http://www.nr.com.