In many important practical situations, beam propagation through a nonlinear media can be accurately modeled by the following nonlinear Schrödinger equation (NLS)
i Q_{η} + Δ_{⊥}Q + Q^{2} Q = 0  (NLS) 
with a function Q_{0}(ρ) as an incident radially symmetric pulse Q(ρ, η=0)= Q_{0}(ρ) (here Q_{η} = ∂Q / ∂_{η} and the Laplacian Δ_{⊥} is given by Δ_{⊥}= ∂^{2}/∂ρ^{2} + (1/ρ)(∂/∂ρ), η=z / L_{df} is a normalized propagation distance, and ρ= (x^{2}+y^{2})^{1/2} / R_{0} is the normalized distance from the center of the pulse in the transverse twodimensional plane, L_{df} is the diffraction length and R_{0} is the beam’s FW1/eM radius).
However, this model is prone to exponential growth or even singularity formation of a numerical solution in the focusing regions when selffocusing is a dominating process. Another model for nonlinear propagation, which proved to better handle selffocusing regimes [6], is a nonparaxial approximation (nopaxNLS) which, compared to NLS, has an extra term of the pulse second derivative
i Q_{η} + ē Q_{ηη} + Δ_{⊥}Q + Q^{2} Q = 0  (nonpaxNLS) 
with the following coefficient of nonparaxiality ē = (λ_{0} / 2πR_{0})^{2} for the incident pulse Q(ρ, η=0)= Q_{0}(ρ) with the radius R_{0} and wavelength λ_{0}. Both equations are derived from the following nonlinear Helmholtz (NHL) equation [7,8].
E_{zz}(x,y,z,t) + Δ_{⊥}E + (k_{0})^{2}(1+(2n_{2}/n_{0})E^{2}) E = 0  (NHL) 
by assuming no backscattering happens while pulse propagate the nonlinear sample. The model for such pulse can be represented by the following SVEA form with the complex conjugate term Q*(x,y,z,t), representing the backscattered part of the pulse, left out
E(x,y,z,t) = R_{0} k_{0} (2n_{2}/n_{0})^{1/2} Q(x,y,z,t) exp[ i (ω_{0}t– k_{0}z) ] 
If λ_{0} << R_{0} , which is the case in many situations, the nonparaxiality coefficient ē is negligibly small, so that nonparaxial approximation turns to NLS – paraxial approximation. NLS can be solved numerically as an example of an evolution problem (eta acts as the time variable in this case).
Problems with a numerical solution of paraxial NLS arise when the initial power of the incident pulse N_{0} = (Q_{0}_{2})^{2}=∫ Q_{0}(ρ)^{2} ρ dρ exceeds a threshold power N_{c} = 1.86 [1,2], which happen in selffocusing regimes. What may happen in such cases is that the nonparaxial term ēQ_{ηη} increases substantially to become as important as other terms in nonparaxial NLS, so that disregarding this term is no longer safe. Nonparaxial term plays as a balancing term, without which the Laplacian (diffraction) and the cubic nonlinearity (Kerr factor) become misbalanced in paraxial NLS, which leads to possible numerical collapses [3,4,5].
One solution would be to relax the nonparaxiality coefficient, introducing a linear damping ε which “turns on” the nonparaxial term near selffocusing regions only [6,7,8]. An extension of splitstep method based on Padé approximation operators [9] allows solving nonparaxial NLS, where the nonparaxiality coefficient ε is adapted at each propagation step based on the behavior of the current solution.
See also: Paraxial approximation, Slowly varying envelope approximation (SVEA).
References  
[1]  G. Fibich, A. Gaeta, Critical power for selffocusing in bulk media and in hollow waveguides, Opt. Lett. 25 (2000) 335–337. 
[2]  M.I. Weinstein, Nonlinear Schrödinger equations and sharp interpolation estimate, Comm. Math. Phys. 87 (1983) 567–576. 
[3]  N. Akhmediev, A. Ankiewicz, J.M. SotoCrespo, Does the nonlinear Schrödinger equation correctly describe beam propagation?, Opt. Lett. 18 (1993) 411–413. 
[4]  M.D. Feit, J.A. Fleck, Beam nonparaxiality, filament formation, and beam breakup in the selffocusing of optical beams, J. Opt. Soc. Amer. B Opt. Phys. 5 (1988) 633–640. 
[5]  G. Fibich, Small beam nonparaxiality arrests selffocusing of optical beams, Phys. Rev. Lett. 76 (1996) 4356–4359. 
[6]  G. Fibich, B. Ilan, S. Tsynkov, Backscattering and nonparaxiality arrest collapse of damped nonlinear waves, SIAM J. Appl. Math. 63 (5), (2003), 17181736. 
[7]  G. Fibich A and S. Tsynkov B, Numerical solution of the nonlinear Helmholtz equation using nonorthogonal expansions, J. Comp. Phys. 210, (2005), 183224. 
[8]  G. Baruch, G. Fibich, and Semyon Tsynkov, Simulations of the nonlinear Helmholtz equation: arrest of beam collapse, nonparaxial solitons and counterpropagating beams, Opt. Express 16, 1332313329 (2008). 
[9]  K. Malakuti and E. Parilov, A splitstep difference method for nonparaxial nonlinear Schrödinger equation at critical dimension, Appl. Num. Math. 61(7), 891899 (2011). 
SimphoSOFT® propagation mathematical model is based on paraxaial approximation. 
SimphoSOFT® can be purchased as a single program and can be also configured with Energy Transfer addon , MultiBeam addon , Optimization addon , Zscan addon , and MPA Info+ addon for an additional charge. Please, contact our sales staff for more information


You can request a remote evaluation of our software by filling out online form  no questions asked to get it, no need to install it on your computer. 