Background
The vector version of the SSPROP solves the coupled nonlinear Schrödinger equations for propagation in a birefringent fiber. The code can model birefringence, differential group delay (PMD), polarizationdependent dispersion, and polarization dependent loss, all in the context of nonlinear propagation.
The user may choose from two different algorithms, depending on whether the birefringent beat length is shorter or longer than the nonlinear length.
Birefringence
In general, the birefringent axes of an optical fiber may not be oriented in the x and y directions, but in some other arbitrary direction ψ. Moreover, the two orthogonal eigenstates of the fiber may not even be linearly polarized — they could be circularly or even elliptically polarized. This would be the case, for example, in fiber that is twisted or spun during or after fabrication. To handle the most general case, SSPROP allows the user to separately specify not only the dispersion β(ω) and loss (α) for each of the two eigenstates, but also the exact polarization states to which these coefficients apply.
The most general elliptical polarization state can be described by two angular parameters, ψ and χ. As depicted in the figure below, ψdescribes the angle that the polarization ellipse makes with the xaxis and χ is an anglar quantity that describes the degree of ellipticity.
Positive values of χ correspond to righthanded polarization states while negative values of χ are lefthanded polarization states. χ = 0 corresponds to linear polarization while χ = π/4 is circularly polarized. On the Poincaré sphere, 2ψ and 2χ describe the longitude and lattitude of the principal eigenstate, respectively.
When specifying the eigenstates of the fiber, it is sufficient to give ψ and χ for one eigenstate because the second eigenstate is known to be orthogonal to the first.
Elliptical Basis Method
Any polarization state may be decomposed into a linear combination of the two orthogonal eigenstates of the fiber, which we label “a” and “b”. If u_{x} and u_{x} represent the two components of the electric field vector in the xy basis (i.e., the Jones vector), then the corresponding components u_{a} and u_{b} can be calculated using the unitary transformation:
where ψ and χ describe the principal eigenstate (“a”) of the fiber. In this new basis, the linear portion of the wave equations for u_{a} and u_{b} are decoupled. The linear portion of the propagation can therefore be performed separately on u_{a} and u_{b} in the spectral domain, using a technique analogous to that used in the scalar case.
where h_{a} and h_{b} are given by:
When performing the nonlinear part of the propagation, the appropriate coupled nonlinear equations (with linear terms omitted) are [Menyuk, JQE 1989]:
where χ quantifies the degree of ellipticity of the eigenstates as described above. The (…) terms in the above expression denote additional nonlinear terms that average to zero when the birefringent beat length is much shorter than the nonlinear length. These additional terms are also identically zero when the eigenstates are circularly polarized.
After propagating through the desired number of steps, the final solution can be rotated from the elliptical basis (u_{a}, u_{a}) back into the Jones basis (u_{x},u_{y}) by using the inverse transformation:
Circular Basis Method
When the fiber birefringence is small, i.e., when the beat length is comparable to or larger than the nonlinear length, the additional terms in the nonlinear equations cannot be neglected. In this case, it is necessary to decompose the field into left and righthand circular polarization components before computing the nonlinear propagation. The circular components u_{+} and u_{–} can be computed from u_{x} and u_{y} using the following unitary transformation:
With this transformation, the coupled nonlinear equations for u_{+} and u_{–} become (again omitting linear terms):
where in this case no additional nonlinear terms have been neglected. Because the eigenstates of the fiber are not in general circularly polarized, the linear portion of the propagation is not as simple in the circular basis. After some algebra, one finds that the linear propagation can be be computed in the spectral domain, using the following matrix multiplication:
where the matrix elements h_{nm} are given by:
and h_{a} and h_{b} are the same quantities given earlier in the context of the elliptical basis method.
After propagating through the desired number of steps, the final solution can be rotated from the circular basis (u_{+}, u_{–}) back into the Jones basis (u_{x},u_{y}) by using the inverse transformation:
The circular basis method is more accurate than the elliptical basis method because it does not neglect any nonlinear terms. The disadvantage of the circular method is that the stepsize dz must always be much smaller than the beat length in order to produce meaningful results. If the beat length is smaller than the nonlinear length, this requirement forces one to use a stepsize that is much smaller than the nonlinearity would otherwise dictate.
Syntax
A summary of the syntax and usage can be obtained from Matlab by typing “help sspropv” or “help sspropvc“.
The compiled mex file (sspropvc) can be invoked from Matlab using one of the following forms:
u1 = sspropvc(u0x,u0y,dt,dz,nz,alphaa,alphab,betapa,betapb,gamma);
u1 = sspropvc(u0x,u0y,dt,dz,nz,alphaa,alphab,betapa,betapb,gamma,psp,method);
u1 = sspropvc(u0x,u0y,dt,dz,nz,alphaa,alphab,betapa,betapb,gamma,psp,method,maxiter);
u1 = sspropvc(u0x,u0y,dt,dz,nz,alphaa,alphab,betapa,betapb,gamma,psp,method,maxiter,tol);
The last four arguments assume a default value if they are left unspecified. The corresponding Matlab mfile can be invoked using a similar syntax by replacing sspropvc with sspropv.
sspropvc may also be invoked with a single input argument, to specify options specific to the FFTW routines (discussed below):
sspropvc option
Input Arguments
u0x, u0y

vector (N)  Input optical field, specified by two lengthN vector time sequences. u0x represents the xcomponent of the complex, slowlyvarying envelope of the optical field, and u0y represents the corresponding ycomponent. The fields should be normalized so that u0x^2 + u0y^2 is the optical power. 
dt  scalar  The time increment between adjacent points in the vector u0. 
dz  scalar  The stepsize to use for propagation 
nz  scalar (int)  The number of steps to take. The total distance propagated is therefore L = nz*dz 
alphaa, alphab  scalar or vector (N)  The linear power attenuation coefficients for the two eigenstates of the fiber. Here we use the labels “a” and “b” to denote the two eigenstates, which need not coincide with the xy axes. Polarization dependent loss is modeled by using different numbers for alphaa and alphab.The loss coefficient may optionally be specified as a vector of the same length as u0x, in which case it will be treated as vector that describes a wavelengthdependent loss coefficient α(ω) in the frequency domain. (The function wspace.m in the tools subdirectory can be used to construct a vector with the corresponding frequencies.) 
betapa, betapb  vector  Realvalued vectors that specify the dispersion for each eigenstate (a, b) of the fiber. The dispersion can be specified to any polynomial order by using a betap vector of the appropriate length.Birefringence is accomodated by making the first elements betapa(1) and betapb(1) unequal. Differential group delay, or polarization mode dispersion is likewise treated by making the second elements betapa(2) and betapb(2) different. (See note below for a more complete discussion.)The propagation constant can also be specified directly by replacing the polynomial argument betap with a vector of the same length as u0x. In this case, the argument betap is treated as a vector describing propagation constant β(ω) in the frequency domain. (The function wspace.m in the tools subdirectory can be used to construct a vector with the corresponding frequencies.) 
gamma  scalar  A real number that describes the nonlinear coefficient of the fiber, which is related to the mode effective area and the nonlinear refractive index n_{2}. 
psp  scalar or vector (2)  Principal eigenstate of the fiber, specified as a 2vector containing the angles ψ and χ (see discussion above), psp = [ψ ,χ].If psp is a scalar, it is interpreted to be ψ, and χ is then taken to be zero. This corresponds to a linearlybirefringent fiber whose axes are oriented at an angle χ with respect to the xy axes.If psp is left completely unspecified, it assumes a default value of [0,0], which means that the fiber eigenstates are linearly polarized along the x and y directions. 
method  string  String that specifies which method to use when performing the splitstep calculations. The following methods are recognized “elliptical” or “circular”.When method = “elliptical”, sspropv will solve the equations by decomposing the input field into the (in general) elliptical eigenstates of the fiber. This method is appropriate only in fibers where the birefringent beat length is much shorter than the nonlinear length.When method = “circular”, sspropv will instead solve the equations by decomposing the input field into a right and leftcircular basis. This method is more accurate, but requires that the step size be small compared to the beat length. 
maxiter  scalar (int)  The maximum number of iterations to make per step. If the solution does not converge to the desired tolerance within this number of iterations, a warning message will be generated. Usually this means that the chosen stepsize was too small. (default = 4) 
tol  scalar  Convergence tolerance: controls to what level the solution must converge when performing the symmetrized splitstep iterations in each step. (default = 10^{–5}.) 
Output Arguments
u1x, u1y

vector (N)  Output optical field, specified as two lengthN vectors. 
Options
Several internal options of the routine can be controlled by separately invoking sspropvc with a single argument:
sspropvc savewisdom
sspropvc forgetwisdom
sspropvc loadwisdom
The first command will save the accumulated FFTW wisdom to a file that can be later used. The second command causes sspropc to forget all of the accumulated wisdom. The last command forces FFTW to load the wisdom file from the current directory. The wisdom file (if it exists) is automatically loaded the first time sspropvc is executed. The name of the wisdom file is “fftwwisdom.dat” for the doubleprecision version of the program and “fftwfwisdom.dat” for the singleprecision version. This can be changed by recompiling the code. The wisdom files can and are shared between the vector and scalar versions of SSPROP. Note that the wisdom files are platform and machinespecific. You should not expect optimal performance if you use wisdom files that were generated on a different computer.
The following four commands can be used to designate the planner method used by the FFTW routines in subsequent calls to sspropc.
sspropvc estimate
sspropvc measure
sspropvc patient
sspropvc exhaustive
The default method is patient. These settings are reset when the function is cleared or when Matlab is restarted.
These options are only available in the compiled version of the routine.
Notes
Slowlyvarying Envelope: In the scalar version of SSPROP, is is customary to factor out the rapidly oscillating terms exp(i(β_{0}z – ωt)) from the field in order to obtain an equation for the slowlyvarying envelope. In SSPROP, this is achieved by setting the first argument of the dispersion polynomial betap(1) equal to 0. In a fiber that has birefringence, it is no longer clear how to factor out these rapid oscillations: should we use β_{0x} or β_{0y}? One approach is to factor out exp(iβ_{0x}z) from the xcomponent of the field and exp(iβ_{0y}z) from the ycomponent of the field. However, with this definition we can no longer regard u_{0x} and u_{0y} as a Jones vector that describes the polarization state. Therefore, we instead choose to factor out a common phase exp(iβ_{0avg}z) variation from both components of the field. Provided we choose β_{0avg }to be the average of β_{0x} and β_{0y}, the resulting fields u_{x} and u_{y} will still be slowlyvarying envelopes that describe the instantaneous Jones vector of the optical signal. In SSPROP, this is accomplished numerically by choosing betapa(1) and betapb(1) to be equal and opposite such that betapa(1) – betapb(1) = Δβ_{0}.
Moving Reference Frames: A similar consideration applies to the difference in group velocity. In a birefringent fiber, the group velocities can be different for the x and y polarizations. Therefore we solve the nonlinear Schrodinger equations in a reference frame that is moving at a velocity in between v_{x} and v_{y}. This amounts to making a change of varibles T = t – β_{1avg}z, where β_{1avg} is the average value of β_{1}. In SSPROP, this is accomplished numerically by choosing betapa(2) and betapb(2) to be equal and opposite such that betapa(2) – betapb(2) = Δβ_{1}.
Units and Dimensions: The dimensions of the input and output quantities are arbitrary, as long as they are self consistent. For examle, if u0^{2} has dimensions of Watts and dz has dimensions of meters, then the nonlinearity parameter gamma should be specified in W^{1}m^{1}. Similarly, if dt is given in picoseconds, and dz is given in meters, then the dispersion polynomial coefficients betap(n) should have dimensions of ps^{(n1)}/m. It is of course possible to solve the normalized dimensionless nonlinear Schrödinger equation by setting some of the input terms to 1 or –1 as appropriate.
Periodicity: SSPROP uses the FFT (DFT) to calculate the spectrum, which implies that the input and output signals are periodic in time. The periodicity is determined by the time increment and the length of the input vector, T = dt*length(u0x). Because of the periodic boundary conditions used by the DFT, care must be taken to ensure that if the optical field at the edges of the window is not negligible it must be continuous in both magnitude and phase.
Iterations and Tolerance: The last two optional parameters, maxiter and tol are related to the symmetrized splitstep iteration algorithm. The algorithm uses a trapazoid integration equation to approximate the effect of the nonlinearity over a distance dz, but this approximation requires knowledge of the field at the subsequent distancestep. This problem is solved by using an iterative approach. maxiter represents the maximum number of iterations performed per step, and tol is a positive dimensionless number that tells the algorithm what level of convergence is required before the iteration stops.