We improve the accuracy of numerical simulations for short fiber optical parametric amplifiers (OPAs). Instead of using the usual coarse-step method, we adopt a model for birefringence and dispersion which uses fine-step variations of the parameters. We also improve the split-step Fourier method by exactly treating the nonlinear ellipse rotation terms. We find that results obtained this way for two-pump OPAs can be significantly different from those obtained by using the usual coarse-step fiber model, and/or neglecting ellipse rotation terms.
©2008 Optical Society of America
Numerical simulation of propagation along optical fibers has been developed over the past twenty years. Since the main application has been for optical communication over long fiber links, spanning hundreds or even thousands of kilometers, approximations appropriate to that regime have been used extensively  In particular, in such simulations “coarse-step” variations of the linear parameters of the fiber are used, making approximations in two areas:
(1) The linear properties of the fiber are modeled by dividing it into segments of length dz, assumed to be uniform. This segment is much sorter than the fiber length, but larger than or of the order of the correlation length of the linear parameters. The birefringence and dispersion parameters in these segments are obtained by assuming that they are random variables with suitable probability density functions (pdf’s). Because dz is larger than the correlation length of the linear parameters of the fiber, this approach is called the “coarse-step” method .
It is common that typical fiber lengths used for fiber optical parametric amplifiers (OPAs) are much smaller, ranging from a few kilometers, to as little as a few meters . In this situation, the assumptions made for simulation of propagation in long fibers using the coarse-step method may very well become inappropriate, and it is therefore important to examine their range of validity, and to develop alternate models appropriate for short fiber OPAs.
In this paper, we develop an alternative model for fibers with longitudinal variations of birefringence and dispersion, applying realistic fine-step variations of the random parameters, as well as including the ER terms in the basic equations. In particular, we apply the so-called second Wai-Menyuk model (WMM)  to emulate the fine-step variations of the linear parameters.
In , Section V, it is shown that when the fiber length is much larger than the correlations length of the linear parameters, then simulations using the fine-step variations of the parameters given by the WMM give the same results, in a statistical sense, than the simulations using the coarse-step method, provided that a correct interpretation is made of how to relate the parameters used in the different models. But since OPAs use fiber lengths not necessarily larger than the correlation length, the equivalence between both models fails and realistic fine-step variations of the liner parameters must be used. Moreover, neglecting ER terms leads to inaccurate results, which in short-fiber OPAs can be orders of magnitude different than when properly considering all the terms of the equations. The main motivation of this article is to demonstrate these last two points by developing a model suitable for short-fiber OPAs simulations.
This article is organized as follows: in Section 2 we introduce the model for accurately describing the linear properties of short fibers. In Section 3 we present a method for solving the nonlinear Schrodinger equation in a birefringent fiber by means of the split-step Fourier method (SSFM) when the ER terms are kept. In Section 4 we present simulation results for the gain spectra of short fiber OPAs, illustrating the importance of the new model for the linear properties, and of keeping the ER terms. We conclude in Section 6.
2. Modeling the linear properties of fibers with longitudinal variations.
It has been shown that both zero-dispersion-wavelength (ZDW) fluctuations and polarization-mode dispersion (PMD) have a strong impact on OPA gain behavior [6,7]. Hence it is important to have a suitable model for describing random longitudinal variations of birefringence and dispersion in fibers. When dealing with fibers with such variations, one can define two characteristic lengths, respectively associated with birefringence and dispersion:
- Birefringence. One defines the birefringence correlation length L b as the maximum distance along z between two points after which it is no longer possible to predict, within any level of certainty, the birefringence strength and orientation at the second point when knowing those at the first one. L b is determined by the longitudinal variations of the linear birefringence, Δn(z), and the orientation of the axes of linear birefringence, characterized by the angle θ(z) with respect to a fixed x axis. The random fluctuations of these parameters give rise to the well-known phenomenon of PMD, extensively studied because of its impact on optical-fiber communication systems.
- Dispersion. Similarly, the dispersion is characterized by the variations of the ZDW, λ 0(z). The autocorrelation of λ 0(z) has a width L d, the dispersion correlation length, which is the scale length of the dispersion fluctuations.
For short OPAs, the fiber length L can be of the same order as, or even shorter than L d and L b. Then it is clear that the usual coarse-step model for describing birefringence and dispersion in numerical calculations is no longer applicable. The reason is that the standard coarse-step model considers uniform segments of fiber of length dz≈L d, L b≪L, in which the dispersion and birefringence are chosen at random, according to some pdf’s (for instance a Gaussian pdf is used in ). Also, a random rotation angle is used between segments. The resulting fiber structure is not very physical, because a real fiber should exhibit continuously varying parameters, not completely discontinuous ones. Nevertheless this model yields accurate results for long fibers, because of the efficient averaging occurring over many correlation lengths, provided that one correctly interprets the relation between the input parameters of the simulation and the real fiber parameters (see , section V).
To obtain smooth variations of the parameters of the different segments in short fibers, we use dz≪Ld, Lb. For each segment of length dz the values of Δn(z), θ(z) and λ 0(z) are obtained by applying the second Wai-Menyuk model  with a given L b. This model is described in the Appendix.
As an example, in Fig. 1 we show θ(z) and Δn(z) in two different situations: when they change abruptly at each Lb, or when the Wai-Menyuk model is applied. In both cases we use Lb = 50 m and an average birefringence Δn = 10-7.
In the simulations of Section 4 we will use fine-scale variations of both Δn(z) and θ(z), as shown in Fig. 1, demonstrating that this more accurate modeling significantly changes the results given by the coarse-step method.
A similar procedure can be used to obtain λ 0(z) Although not shown in this paper, we have performed preliminary simulations showing that the same happens if a fine-step Langevin process is applied for the modeling of ZDW fluctuations. We will not study ZDW fluctuations further in this paper; however the preceding reasoning already indicates that a fine-step model should also be applied to ZDW when modeling its fluctuations.
3. Keeping the ellipse rotation terms in the nonlinear propagation equations
A powerful approach for studying light propagation in fibers is the coupled nonlinear Schrodinger equation (C-NLSE) and its generalizations, which are generally solved by means of the split-step Fourier method (SSFM). This method has the advantage that it can handle simultaneously a large number of frequency components, which might be present in communication systems with several optical carriers, modulated carriers, and/or modulated pumps.
In fibers that exhibit constant linear birefringence the x and y components of the electric fields can be written as
where Cx (z, t) and Cy (z, t) are the slowly-varying envelopes (SVEs) of the electric field components, and
(see the Appendix for definitions). The total SVE of the electric field vector is obtained by letting C⃗ = Cxx̂+Cyŷ, where x̂ and ŷ are the unit vectors along the x and y axes.
Then the time-domain C-NLSE takes the form
where Px = CxCx* and Py = CyCy*. γ is the fiber nonlinearity coefficient. β (n) x,y represents the nth derivative of βx,y at the center of the optical spectrum, ω c; and Δβxy = βx(ωc)-βy(ωc).
The C-NLSE is generally numerically integrated by the SSFM, schematically represented in Fig. 2
The fiber is divided into a large number Nz of equal segments, with length dz = L/Nz. In each segment, the elementary contributions due to dispersion (including birefringence if present) and nonlinearity are calculated sequentially. For this reason, this is known as a split-step procedure. If birefringence and dispersion have random variations along the fiber, these can be introduced by varying the parameters in consecutive segments, as described in Section 2 We assume that C⃗(z, t) is known at the end of the preceding section. Since linear birefringence and dispersion are best handled in the frequency domain, the Fourier transform of C⃗(z, t), D⃗(z, Ω), is calculated by means of a fast Fourier transform (FFT). (ω is the angular frequency measured from the center frequency, i.e. Ω=ω-ωc.) We then have
where Φ0 = Δnωc dz/2c = πΔndz/λc. β and βc can be expanded as usual in terms of the desired orders of dispersion, generally up to the fourth order for fiber OPAs. (See the Appendix for more details.)
In the section with Kerr nonlinearity alone, we have
In long fiber transmission the terms containing (Cy)2 Cx* and (Cx)2 Cy* can be neglected, because they lead to phase factors of the form exp[±2i(βx-βy)z]. These terms oscillate rapidly due to the difference in phase velocity between the two axes, and essentially average out to zero over several beat lengths. These terms are also associated with the phenomenon of ellipse rotation, and so neglecting them amounts to neglecting ellipse rotation. (These terms are also known as coherent-coupling terms). In such situations, we then keep only those terms on the right-hand sides that are synchronous. We then obtain the well-known equations
Thus with this approach the nonlinear effect on each linear field component is simply the multiplication by a phase factor. These exponential transformations are valid for any length dz, and are easy to implement in numerical simulations. For these reasons they are widely used in code for simulating propagation in optical fibers.
However, neglecting (Cy)2 Cx* and (Cx)2 Cy* in Eq. (8) may not be justified, or may not be a very good approximation, as in the following cases. For example, if one wants the numerical method to work in the limit of an isotropic fiber, then these terms must be kept, because otherwise an incorrect limit is obtained when using circularly-polarized (CP) states of polarization (SOPs). Also, in a fiber with linear birefringence, if dz is small, and/or the birefringence is weak, it may not be possible to argue that these terms will average out to zero. For OPAs in particular, if a high pump power and a fiber with very high γ are used, only a short fiber is needed, and neglecting these terms is probably not justified. It is then necessary to keep them, and to use them in the numerical procedure.
One way to solve Eq. (8) exactly is to switch to a basis of circularly-polarized (CP) states; this leads to an equation which is easy to integrate. Specifically, the unit Jones vectors for right- and left-circularly polarized light are respectively r̂ = (x̂+iŷ)/√2 and l̂ = (x̂-iŷ)/√2. In that basis C⃗ has the components Cr = (Cx+iCy)/√2 and Cl = (Cx-iCy)/√2. Eq. (4) then leads to
where Pr = Cr Cr* and Pl = Cl Cl*. The terms Pr and Pl remain constant. Therefore Eq. (11) can be integrated, with the results
This shows that the circular field components simply experience phase shifts θr and θl when they pass through the nonlinear section. When the phase difference between them, Δθ = θr-θl = 2(Pl-Pr)γdz/3, is not equal to zero, this leads to the well-known effect of ellipse rotation, whereby the major axis of the ellipse of the input SOP is rotated by δθ = Δθ/2.
In practice, the implementation of Eqs. (12) and (13) in an SSFM algorithm is simple. Cx and Cy are used in the step where linear birefringence and dispersion are dealt with. Then one switches to Cr and Cl. Finally one returns to Cx and Cy to work on the next section by means of Cx = (Cr+Cl)/√2 and Cy = (Cr-Cl)/(i√2).
As an alternative, it is also possible to do all the calculations of the Kerr effect in terms of Cx and Cy, without any approximation. Specifically, one can show that
where P 0 = Pz+Py = Pr+Pl is the total power. And
Equation (15) says that the field after the section of length dz is obtained by rotating the input fields by δθ, and by multiplying the result by a phase factor which is calculated as a self phase modulation (SPM) phase shift, which involves the total power regardless of how it is distributed among the SOPs.
In practice, one then has a choice between two approaches for implementing the ellipse rotation: (i) by going from LP to CP states, and back, one only needs to calculate powers and to introduce two phase shifts in the CP basis; (ii) by staying in the LP basis, one also needs to calculate powers, calculate the SPM phase shift and the ellipse rotation angle, and rotate the input fields accordingly. Either way, these operations are all relatively fast compared to those required for the FFT when typical values for the number N of points are used, and so either approach can be implemented without introducing significant computational overhead. For all the simulations we will present in this article we used 214 points, and the simulation time to obtain each gain spectrum was ~35 seconds, independently of the method used.
Finally, after each fiber segment a rotation by a random angle is introduced, to describe the random orientation of the axes of refraction (not shown in Fig. 2). We have described in Section 2 how such rotations should be modeled, for different fiber lengths.
4. Numerical simulations
Here we present the results of numerical simulations showing the separate effects of (i) the fine-step modeling of birefringence, and (ii) keeping the ellipse rotation terms in the SSFM. All the Matlab codes used to generate the following figures, as well as Fig. 1, can be downloaded at .
4.1. Effect of fine-step modeling of birefringence
As an example of the importance of the fine-step modeling proposed here, we simulated an OPA pumped by two orthogonal CP pumps (2P-OPA) . We simulated two cases, with θ(z) and Δn(z) changing either smoothly or abruptly over each segment of length Lb, corresponding to the dashed red and solid blue lines in Fig.1, respectively. We used 1 W of input power for each CW pump (at 1502.6 and 1600.6 nm), γ = 10 W-1km-1 and λ 0x(z)=λ 0y(z)=constant=1550 nm. We also used β (3)=0.1 ps3/km, β (4)=10-4 ps4/km and β (n)=0 for n>4. From Δn, shown for example in Figure 1(b), it is also possible to calculate b (0) and b (1) from (A.7). From these relations, since the mean value of Δn is 10-7, the mean value of b (0) is ~0.2 m-1 and the mean value of b (1) is ~0.001 ps/m. From the relation Dp = 2(2Lb)0.5 b (0) (Dp is the PMD coefficient), and since we used Lb = 50 m, we obtain Dp ~0.6 ps/km0.5. This value is higher than in state of the art highly-nonlinear fiber (HNLF), with Dp ~0.2 ps/km0.5, because we purposely choose a high value of Lc to show that our method is valid under any circumstance. The gain spectra are shown in Fig. 3 (they were obtained by applying the SSFM described in Section 3). Clearly, a more realistic modeling of θ(z) and Δn(z) changes the final result.
In Fig. 1, we used a longitudinal step size dz = 0.5 m for the generation of the birefringence parameters according to the WMM and for the integration of the coupled equations. When applying the coarse-step method, we also use dz = 0.5 m for the integration, but the birefringence parameters vary as shown by the red curves in Fig. 1. This is because we use a low-order SSFM, and the numerical errors associated with it do not enable one to use dz = 50 m. We could use dz = 50 m if we had implemented a higher-order SSFM, or instead of the SSFM used a model considering only four waves in the time domain (as done in  and ), and would obtain the same red curve as shown in Fig. 3 (with no dips near the pumps , which do not appear with the four-wave model).
The important point to be noted is that when numerical errors can be neglected, the different ways of modeling the random birefringence parameters give completely different results. In fact, in all the simulations presented in this paper we checked that the numerical error was always below 0.2 dB at any point of the spectra. In Fig. 3, for example, we needed to use 214 point in a spectral window of 25 THz (time increment of 0.04 ps) around the central frequency (in other words, increasing the number of sample points, or increasing the spectral window or shortening dz do not change the results by an amount larger than 0.2 dB).
We note that the fine-step method yields a gain spectrum which is lower, but more uniform than that predicted by the coarse-gain method. To show that this is not a rare probabilistic event, we performed 100 other realizations of the random process that generated the parameters shown in Fig. 1, and used these new parameters to simulate new gain curves (as we did with Fig. 3). Three random realizations are shown in Fig. 4(a), (b) and (c), while the average gain curve over 50 realizations is shown in Fig. 4(d).
It can be seen that the coarse-step modeling of the linear parameters gives an average gain curve with larger gain variation and polarization dependent gain (PDG). This fact can be understood by noting that: (i) coarse-step methods artificially increase the amount of birefringence, as noted by Marcuse, Menyuk and Wai in , and, consequently also artificially increases the PMD coefficient Dp; (ii) the PMD coefficient is inversely proportional to the depolarization length which characterizes how fast the pumps lose their original orthogonality [1,12]. The depolarization length Ldep is related to Dp by Ldep = 3/(DpΔω)2, where Δω is the angular frequency separation between the pumps. So, in the coarse-step modeling of the linear parameters the pumps lose their original orthogonally faster than when performing proper fine-step modeling, and the gain variation and PDG of the simulated orthogonally-pumped OPA are larger.
4.2. Effect of keeping the ellipse rotation terms
We then performed numerical simulations of a number of 2P-OPAs with orthogonal CP pumps, to illustrate the need for keeping the ER terms. We implemented the SSFM algorithm including the ER terms, and compared the results to those obtained by neglecting these terms. We first simulated propagation in an isotropic fiber, to clearly show the effect of the ER terms, and then included random birefringence, to show how it weakens the influence of the ER terms, particularly in long fibers.
4.2.1. Isotropic fiber
In Fig. 5 we considered the same input and fiber parameters as in the simulations shown in Fig. 3, but used an isotropic fiber [b (0)(z)=b (1)(z)=b (2)(z)=0], where b is the birefringence parameter. As expected, we found that (away from the pumps) the SSFM with ER terms yielded exactly the same gain spectrum as the analytical solution using the four-wave model , while the SSFM without ER terms yielded a significantly inaccurate result.
The gain dip observed around the shortest pump wavelength appears because the signal has the same SOP as this pump; it can be explained through a six-wave model including the ER terms .
4.2.2. Fiber with random birefringence
We then concentrated on the case of OPAs with random birefringence. In sufficiently long fibers, the birefringence tends to average out the ER terms to zero. We considered the same type of OPA as in the preceding section, but added random birefringence, and considered different lengths and pump powers (to keep the same gain). In all cases we used the model introduced in Section 2 for describing the longitudinal variations of birefringence in short fiber OPAs. We used a fiber with θ(z) given by Fig. 1(a), and Δn(z) given by Fig. 1(b).
In Fig. 6 we show the gain spectra of three different two-pump OPAs, with or without the ER terms (solid blue and dashed red curves, respectively), with three different fiber lengths L and total pump powers P 0, namely: L = 200 m and P 0 = 2 W; L = 20 m and P 0 = 20 W; L = 2 m and P 0 = 200 W. Fig. 6 clearly shows again the importance of the ER terms for short fibers: in the case of the 2-m long fiber [Fig. 6(c)] the calculation without the ER terms underestimates the gain in decibels by about 7 dB, which is a very significant difference. Observe from Fig. 1 that in this case [Fig. 6(c)] the birefringence of the fiber is almost constant, since the changes of the birefringence parameters are small in the first 2 m of fiber. This result, together with the one shown in Fig. 5, shows the importance of ER terms in modeling the light propagation not just in isotropic fibers, but also in fibers with constant birefringence. As expected the difference becomes smaller as the length is increased, because of the averaging of ER due to random birefringence. As a result, when the length reaches 200 m [Fig. 6(a)] there is little difference between the results obtained with and without the ER terms.
We have attempted to develop an improved simulation technique for short fiber OPAs based on the WMM model, using a fine longitudinal step model (the second Wai-Menyuk model ) for representing the longitudinal variations of linear birefringence and dispersion, and keeping the nonlinear ellipse rotation terms in the basic equations. We have shown that the ellipse rotation terms can be handled fairly simply when calculating the nonlinear contribution in an SSFM algorithm for light propagation in optical fibers. Keeping these terms ensures that all the effects of the nonlinearity are accounted for. This can be important in situations where the effects of these terms do not necessarily average out to zero, such as in relatively short fibers, including fiber OPAs. We have applied this new SSFM algorithm to a two-pump fiber OPA with orthogonal CP pumps, and shown that it yields results significantly different from those obtained by the conventional SSFM.
In a fiber with linear birefringence, we let x and y denote the directions of the principal axes. βx and βy are the propagation constants for field components polarized along x and y, respectively; they are generally independent functions of the angular frequency ω.
We then define the average propagation constant
β is used to represent dispersion, just as in an isotropic fiber. It is generally approximated by a power-series expansion about ω c. Letting ω=ω-ωc,
where β (m) is the mth derivative of β with respect to ω, evaluated at ω c.
When the fiber has random dispersion variations, these are modeled by providing the pdf for the ZDW, which is directly related to that of β (2). The parameters β (3) and β (4) are usually considered to be constant.
We rewrite βx and βy as
We see that b characterizes the birefringence, because when b = 0 the propagation constants along the two axes are identical, and therefore there is no birefringence. More precisely, if we introduce the effective refractive indices for the two polarizations, nx,y = cβx,y/ω, we see that
where Δn = nx-ny is by definition the fiber birefringence.
We can then expand b in power series about the center frequency ω c. This yields
where b (m) is the mth derivative of b with respect to ω, evaluated at ω c.
If Δn is independent of ω, Eq. (A.5) shows that
Hence under these conditions the birefringence coefficient b and its derivatives do not contribute to chromatic dispersion along the x or y axes, because from Eq. (A.2) the derivatives of βx and βy beyond the first order are not affected by b at all. In particular this means that the fiber then has the same ZDW along both axes. To avoid ZDW fluctuations along z one can also assume that the changes in nx(ω) and ny(ω) at each Lb are given by the mere multiplication of nx(ω) and ny(ω) by some wavelength-independent factor, which changes the values of b (0) and b (1), but does not change β (2), and thus does not produce ZDW fluctuations. This modeling enabled the study, in references  and , of the separate effects of PMD and ZDW fluctuations on OPA gain. We use the same approach in this paper.
[However, if Δn is wavelength-dependent, or the spectral shapes nx(ω) and ny(ω) change along z, both effects (PMD and ZDW fluctuations) are coupled, and proper simulations should consider the influence of PMD on ZDW fluctuations.]
The work of M. Marhic and C. Braimiotis was supported in part by the European Regional Development Fund, through the Welsh Assembly Government, and by a grant from the UK’s EPSRC.
References and links
1. D. Marcuse, C. R. Menyuk, and P. K. A. Wai, “Application of the Manakov-PMD Equation to Studies of Signal Propagation in Optical Fibers with Randomly Varying Birefringence,” J. Lightwave Technol. 15, 1735–1746 (1997). [CrossRef]
2. C. R. Menyuk, “Nonlinear pulse propagation in birefringent optical fibers,” IEEE J. Quantum. Electron . QE-23, 174–176 (1987). [CrossRef]
3. A. Carena, V. Curri, R. Gaudino, P. Poggiolini, and S. Benedetto, “A time-domain optical transmission system simulation package accounting for nonlinear and polarization-related effects in fiber,” IEEE J. Sel. Areas Commun. 15, 751–765 (1997). [CrossRef]
4. M. E. Marhic, G. M. Williams, L. Goldberg, and J. M. P. Delavaux, “Tunable fiber optical parametric wavelength converter with 900 mW of CW output power at 1665 nm,” Proc. SPIE , 6103 0W-1-0W-12 (2006). [CrossRef]
5. P. K. A. Wai and C. R. Menyuk, “Polarization mode dispersion, decorrelation, and diffusion in optical fibers with randomly varying birefringence,” J. Lightwave Technol. 14, 148–157 (1996). [CrossRef]
6. F. Yaman, Q. Lin, S. Radic, and G. P. Agrawal, “Impact of dispersion fluctuations on dual-pump fiber-optic parametric amplifiers,” IEEE Photon. Technol. Lett. 16, 1292–1294 (2004). [CrossRef]
7. F. Yaman, Q. Lin, and G. P. Agrawal, “Effects of polarization-mode dispersion in dual-pump fiber-optic parametric amplifiers,” IEEE Photon. Technol. Lett. 16, 431–433 (2004). [CrossRef]
10. M. E. Marhic, K. K. Y. Wong, and L. G. Kazovsky, “Fiber optical parametric amplifiers with linearly or circularly polarized waves,” J. Opt. Soc. Amer. B 20, 2425–2433 (2003). [CrossRef]
11. M. E. Marhic, A. A. Rieznik, and H. L. Fragnito, “Investigation of the gain spectrum near the pumps of two-pump fiber-optic parametric amplifiers,” J. Opt. Soc Amer. B 25, 22–30 (2008). [CrossRef]
12. M. Karlsson and J. Brentel, “Autocorrelation function of the polarization-mode dispersion vector,” Opt. Lett. 24, 939–941 (1999). [CrossRef]