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

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 Δn(z), θ(z) and λ ₀(z) are obtained by applying the second Wai-Menyuk model [5] 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 L_b, or when the Wai-Menyuk model is applied. In both cases we use L_b = 50 m and an average birefringence Δn = 10^-7.

 

Fig. 1. (a) Total angle θ(z) along a fiber with L = 200 m and L_b = 50 m, using a coarse-step model (dashed red) or the Wai-Menyuk model (solid blue); both curves are such that their values are the same after each L_b. (b) Variation of Δn(z) along the same fiber, obtained in a similar manner.

Download Full Size | PDF

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 λ ₀(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 C_x (z, t) and C_y (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⃗ = C_xx̂+C_yŷ, where x̂ and ŷ are the unit vectors along the x and y axes.

Then the time-domain C-NLSE takes the form

where P_x = C_xC_x* and P_y = C_yC_y*. γ is the fiber nonlinearity coefficient. β ⁽ⁿ⁾ _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

 

Fig. 2. Schematic representation of the structure of the SSFM algorithm in one fiber segment. D is the dispersion; LBR=linear birefringence; NL=nonlinearity; FT=Fourier transform; IFT=inverse FT. L=>C and C=>L indicate transformation from a linear to a circular basis and vice-versa, respectively.

Download Full Size | PDF

The fiber is divided into a large number N_z of equal segments, with length d_z = 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 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 Φ₀ = Δ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 (C_y)² C_x* and (C_x)² C_y* 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

Eq. (9) implies that P_x and P_y remain constant. Then Eq. (9) can easily be integrated, which yields

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 (C_y)² C_x* and (C_x)² C_y* 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 C_r = (C_x+iC_y)/√2 and C_l = (C_x-iC_y)/√2. Eq. (4) then leads to

where P_r = C_r C_r* and P_l = C_l C_l*. 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) in an SSFM algorithm is simple. C_x and C_y are used in the step where linear birefringence and dispersion are dealt with. Then one switches to C_r and C_l. Finally one returns to C_x and C_y to work on the next section by means of C_x = (C_r+C_l)/√2 and C_y = (C_r-C_l)/(i√2).

As an alternative, it is also possible to do all the calculations of the Kerr effect in terms of C_x and C_y, without any approximation. Specifically, one can show that

where P ₀ = 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¹⁴ 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 [9].

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) [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^-1km^-1 and λ _0x(z)=λ _0y(z)=constant=1550 nm. We also used β ⁽³⁾=0.1 ps³/km, β ⁽⁴⁾=10^-4 ps⁴/km and β ⁽ⁿ⁾=0 for n>4. From Δn, shown for example in Figure 1(b), it is also possible to calculate b ⁽⁰⁾ and b ⁽¹⁾ from (A.7). From these relations, since the mean value of Δn is 10^-7, the mean value of b ⁽⁰⁾ is ~0.2 m^-1 and the mean value of b ⁽¹⁾ is ~0.001 ps/m. From the relation D_p = 2(2L_b)^0.5 b ⁽⁰⁾ (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 d_z = 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 Δn(z) curves in Figs. 1(a) and 1(b).

Download Full Size | PDF

 

Fig. 4. (a), (b) and (c): Gain spectra for an OPA with two orthogonal CP pumps, as Fig. 3 but with three other realizations of the random process for the determination of the linear parameters of the fiber. (d): Average gains after 50 realizations of the random process. Solid blue: fine-step method for the linear parameters. dashed red: coarse-step method.

Download Full Size | PDF

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¹⁴ 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Δω)², 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.

 

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.

Download Full Size | PDF

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 ⁽⁰⁾(z)=b ⁽¹⁾(z)=b ⁽²⁾(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 [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].

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 ₀, namely: L = 200 m and P ₀ = 2 W; L = 20 m and P ₀ = 20 W; L = 2 m and P ₀ = 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.

 

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) and 1(b) (solid blue lines). The input signal has the same SOP as the short-wavelength pump. (a): L = 200 m, P ₀ = 2 W; (b): L = 20 m, P ₀ = 20 W; (c) L = 2 m, P ₀ = 200 W. Observe that the solid blue curve in (a) is the same as in Fig. 3. Also, the horizontal scale is different for each graph, since the gain bandwidth increases considerably for high pump powers.

Download Full Size | PDF

5. 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.