Accurate numerical simulation of short fiber optical parametric amplifiers.

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.


Introduction
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 [1] 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 [1].
(2) The nonlinear terms corresponding to ellipse rotation (ER) are dropped as they essentially average out to zero in each segment [2,3]. (These terms are also known as coherent coupling terms.) 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 [4]. In this situation, the assumptions made for simulation of propagation in long fibers using the coarsestep 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) [5] to emulate the fine-step variations of the linear parameters.
In [5], 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 shortfiber 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.

Modeling the linear properties of fibers with longitudinal variations.
It has been shown that both zero-dispersion-wavelength (ZDW) fluctuations and polarizationmode 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:

1)
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.

2)
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 [8]). 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 [1], section V).
To obtain smooth variations of the parameters of the different segments in short fibers, we use dz << L d , L b . For each segment of length dz the values of 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.

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 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 where ˆ x and ˆ y are the unit vectors along the x and y axes. Then the time-domain C-NLSE takes the form where . γ is the fiber nonlinearity coefficient. The fiber is divided into a large number N z of equal segments, with length dz = L/N z . 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 ) , ( t z C is known at the end of the preceding section. Since linear birefringence and dispersion are best handled in the frequency domain, the Fourier transform , is calculated by means of a fast Fourier transform (FFT). (Ω is the angular frequency measured from the center frequency, i.e. . β 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 * 2 ) ( 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 Eq. (9) implies that P x and P y remain constant. Then Eq. (9) can easily be integrated, which 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 * 2 ) ( 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ˆ y )/ 2 and ˆ l = (ˆ x −iˆ y )/ 2 . In that basis C has the components Eq. (4) then leads to where * r r r C C P = and * l l l C C P = . The terms P r and P l 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 P l − P r ( )γ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)  As an alternative, it is also possible to do all the calculations of the Kerr effect in terms of x C and y C , without any approximation. Specifically, one can show that where P 0 = P z + P y = P r + P l 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 2 14 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.

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 [9].

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) [10]. We simulated two cases, with θ (z) and ∆ n(z) changing either smoothly or abruptly over each segment of length L b , 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 -1 km -1 and λ 0x (z) = λ 0y (z) = constant = 1550 nm. We also used β (3) = 0.1 ps 3 /km, β (4) = 10 -4 ps 4 /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 D p = 2(2L b ) 0.5 b (0) (D p is the PMD coefficient), and since we used L b = 50 m, we obtain D p ~ 0.6 ps/km 0.5 . This value is higher than in state of the art highly-nonlinear fiber (HNLF), with D p ~ 0.2 ps/km 0.5 , because we purposely choose a high value of L c 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 [6] and [7]), and would obtain the same red curve as shown in Fig. 3 (with no dips near the pumps [11], which do not appear with the four-wave model). Fig. 3. Gain spectra for an OPA with two orthogonal CP pumps, using the total angle rotations shown in Fig. 1(a). The colors match the corresponding θ (z) and 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 2 14 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 [1], and, consequently also artificially increases the PMD coefficient D p ; (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 L dep is related to D p by L dep = 3/(D p ∆ ω ) 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.

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. Fig. 5. Gain spectra for an OPA with two orthogonal CP pumps using an isotropic fiber with the same parameters as in Fig.3. We show the analytical result using the four-wave model (dotted gray), and the results obtained using the SSFM with (solid blue) and without (dashed red) the ER terms. The input signal SOP was the same as that of the short-wavelength pump.

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 [ (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 [10], 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 [11].

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.
(a) (b) (c) Fig. 6. Gain spectra for three different two-pump OPAs with orthogonal CP pumps, with or without the ER terms (solid blue and dashed red curves, respectively). The three OPAs use the same fiber as in Fig. 3, whose birefringence parameters are shown in Fig. 1(a) Fig. 3. Also, the horizontal scale is different for each graph, since the gain bandwidth increases considerably for high pump powers.

Conclusion
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 [5]) 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.

Appendix
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 2 y x β β β + = . (A.1) β is used to represent dispersion, just as in an isotropic fiber. It is generally approximated by a power-series expansion about ω c . Letting c ω ω Ω = − , 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 ( 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 n x (ω) and n y (ω) at each L b are given by the mere multiplication of n x (ω) and n y (ω) by some wavelength-independent factor, which changes the values of (0) b and (1) b , but does not change (2) β , and thus does not produce ZDW fluctuations. This modeling enabled the study, in references [6] and [7], 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 n x (ω) and n y (ω) change along z, both effects (PMD and ZDW fluctuations) are coupled, and proper simulations should consider the influence of PMD on ZDW fluctuations.]