The nonlinear Schrödinger equation (NLS) is the canonical model in nonlinear optics for propagation of intense laser beams in isotropic Kerr media. In the case of propagation through a bulk medium, Kelley [1] used the 2D NLS to predict the possibility of a catastrophic collapse of beams whose input power is above the critical power for collapse. In the case of propagation through planar waveguides, the 1D NLS was used to predict the existence of spatial solitons [2]. Both beam collapse in bulk medium and spatial solitons in planar waveguides were observed in experiments [3, 4]. More recently, configurations of two counter-propagating beams were modeled by two coupled NLS equations [5].

In nonlinear optics, the NLS is derived from the nonlinear Maxwell equations via a series of approximations. First, if the electric field is monochromatic and third harmonic generation is neglected, Maxwell’s equations reduce to the vector nonlinear Helmholtz equation (NLH). If the field is also linearly polarized, the vector NLH reduces to the scalar NLH [6]. Finally, the NLS is derived from the scalar NLH using the paraxial approximation, which is valid when the beam radius is sufficiently large compared with the wavelength. As, however, the 2D NLS predicts that the beam radius shrinks to zero at collapse, the paraxial approximation breaks down at this point. In the case of spatial 1D solitons, the paraxial approximation sets a lower limit on the soliton radius.

The singular behavior of the 2D NLS solutions for collapsing beams is non-physical. Therefore, an important question is whether the singularity formation is already arrested by taking one step back in the aforementioned series of approximations and employing the scalar NLH model, or only in a more comprehensive model. Both the mathematical analysis and simulations of the scalar NLH have proved to be considerably more difficult than for the NLS, since for the NLH one solves a nonlinear boundary-value problem, whereas the NLS requires solving an initial value problem. An additional computational obstacle is that unlike the NLS, which governs the slowly varying envelope, the NLH has to be approximated with sub-wavelength resolution. For these reasons, the question of collapse in the scalar NLH model was not fully answered for over 40 years.

Previously, numerical simulations and asymptotic analysis [7, 8, 9] suggested that nonparaxiality arrests the collapse in bulk medium. These studies, however, applied various simplifying approximations to the scalar NLH. In particular, they considered only forward traveling waves and completely neglected the backscattered field. Even though backscattering is generally believed to be “small”, it may still significant affect the overall propagation, because collapse dynamics in the 2D cubic NLS is extremely sensitive to small perturbations [10].

To study the arrest of collapse in the scalar NLH with no simplifying assumptions (and in particular, with the backscattering included), Fibich and Tsynkov developed a fixed-point iterative numerical method for solving the NLH as genuine boundary value problem [11, 12], which is based on freezing the nonlinearity at each iteration. This method converged for input powers below the critical power for collapse P_cr, but diverged for input powers higher than P_cr. It was unclear, however, whether the divergence above P_cr was due to limitations of the numerical method, or because collapse is not arrested in the scalar NLH model. Subsequently, the method of [11, 12] was used to show numerically the arrest of collapse in the linearly-damped scalar NLH [13]. In that work, however, the magnitude of damping was much larger than in actual physical settings, and could not be reduced to zero. More recently, Sever proved existence of solutions (and hence arrest of collapse) in the scalar NLH with self-adjoint boundary conditions [14]. The proof in [14], however, relies heavily on self-adjointness, whereas propagating fields satisfy radiation boundary conditions (BCs), which are non self-adjoint. Therefore, until now, there has been no conclusive evidence that the collapse is arrested in the scalar NLH model.

 

Fig. 1. The physical setup: A: single beam, B: counter-propagating beams. C: A schematic of the upstream BC at z=Z_max+δ, which freely admits all forward propagating waves (red). D: A schematic of the downstream BC at z=-δ, which freely admits all the backward propagating waves (blue) to pass, and also specifies the (forward moving) incoming beam.

Download Full Size | PDF

In [15], we studied numerically the (0+1)D NLH, which models the propagation of plane waves in a Kerr medium. In this case, the solution always exists, but becomes non-unique (bistable) above a certain input power threshold [16]. Numerically, we observed that the fixedpoint frozen nonlinearity method of [11, 12] converges for low input powers, but diverges for higher powers which are still below the threshold for non-uniqueness. This indicates that the divergence of the fixed-point frozen nonlinearity method is due to the numerical methodology itself, rather than to non-uniqueness or non-existence of the solutions. Therefore, an alternative iterative solver, based on Newton’s method, was constructed and shown to have much better convergence properties. In this Letter, we extend the Newton-based method of [15] to the multi-dimensional case. The resulting technique enables us to solve the (2+1)D NLH for input powers above P_cr. Hence, we obtain the first ever computational evidence that the collapse of the beam is indeed arrested in the scalar NLH model. We also calculate the field backscattered from the domain. Moreover, we solve the (1+1)D NLH for a “nonparaxial” soliton with radius equal to half a wavelength, and observe that it propagates virtually unchanged over 40 diffraction lengths. This indicates that such beams are still in the paraxial regime. Finally, we solve the (1+1)D NLH for two counter-propagating beams and compare the results to those obtained using the coupled NLS model.

The propagation of linearly polarized, continuous wave beams in isotropic Kerr media is governed by the scalar nonlinear Helmholtz equation:

where E is the electric field, k ₀ is the linear wavenumber, n ₀ is the linear index of refraction and n ₂ is the Kerr coefficient. In the bulk medium (2+1)D case x _⊥=(x,y) and Δ_⊥=∂² _x+∂² _y; in the planar waveguide (1+1)D case x _⊥=x and Δ_⊥=∂² _x. We consider an incoming beam traveling in the positive z direction (henceforth “forward” or “right”) impinging on a finite-length Kerr material slab at the z=0 interface and exiting the Kerr medium at the z=Z _max interface, see Fig. 1(A). A portion of the field may be reflected by the interfaces at z=0 or z=Z _max, or backscattered inside the Kerr medium, because of the variations of the index of refraction induced by the forward-propagating beam. To derive the NLS, the standard approach is to represent the field as $E={\mathrm{Ae}}^{{\mathrm{ik}}_{0}z}$ , where the envelope A is assumed slowly varying. Using the standard rescaling x̃ _⊥=x _⊥/r ₀, z̃=z/2L_DF and $\tilde{A}(\tilde{z},{\tilde{\mathbf{x}}}_{\perp })=\sqrt{\frac{{2n}_2}{n_0}}{r}_0{k}_0·A(z,{\mathbf{x}}_{\perp })$ , where r ₀ is the input beam radius and L_DF=k ₀ r ² ₀ is the diffraction length, the NLH can be written in the dimensionless form

where $f^2={\left(r_0{k}_0\right)}^{-2}={\left(\frac{{\lambda }_0}{2\pi r_0}\right)}^2$ is the nonparaxiality parameter. Typically λ ₀≪r ₀ so that f ²≪1 and f ² Ã _z̃z̃≪Ã_z̃. Therefore, the paraxial approximation, which consists of neglecting f ² Ã_z̃z̃, leads to the NLS

In our simulations, the (2+1)D NLH with cylindrical symmetry, i.e., E=E(z, r) where $r=\mid {\mathbf{x}}_{\perp }\mid =\sqrt{x^2+y^2}$ , is approximated with a fourth order finite-difference scheme. The solution is computated for -δ≤z≤Z_max+δ in order to implement the BCs in the linear regions. At the material interfaces z=0 and z=Z _max where the index of refraction is discontinuous, Maxwell equations for a normal-incident field imply that E and E_z are continuous across the interfaces. At z=Z _max+δ we imposed the radiation BC that the field does not have any left-going component for z>Z _max, see Fig. 1(C). Similarly, at z=-δ we implement the two-way radiation BC that for z<0 the field does not have right-going components except for the prescribed incoming beam which impinges on the interface z=0 with a transverse profile E _inc(r), see Fig. 1(D). Because z=-δ and z=Z _max+δ are outside the Kerr slab, the field propagation there is linear, which simplifies the implementation of the radiation BCs, see [11, 12] for more details. The discretized system of nonlinear algebraic equations is solved using Newton’s method [15].

In order to focus on the effects of the Kerr nonlinearity, the values of n ₀ in the Kerr medium (0≤z≤Z _max) and in the surrounding linear medium (z<0 and at z>Z _max) are chosen to be equal, so that to eliminate the reflections due to discontinuity of n ₀ at the interfaces. However, discontinuities in the nonlinear coefficient are not eliminated, and are a source of reflections at z=0,Z _max. Our numerical method can be applied to the case of different n ₀ with no change [15]. Note that, since we solve the NLH in non-dimensional form, simulations in this Letter are valid for any physical value of k₀, n₀, n₂ that corresponds to the same dimensionless quantities f ² and P/P_cr.

The (2+1)D NLH (1) was solved for an incoming collimated Gaussian beam $E_{\mathrm{inc}}={\left(\sqrt{\frac{{2n}_2}{n_0}}{r}_0{k}_0\right)}^{-1}{e}^{-{\left(\frac{r}{r_0}\right)}^2}$ of radius r ₀=1.27λ ₀, corresponding to nonparaxiality parameter of f ²=(k ₀ r ₀)^-2=1/64, and input power of P=1.29P_cr. The NLH solution initially self-focuses, until z≈0.8L_DF where the collapse is arrested, after which the solution defocuses, see Fig. 2(A). The corresponding NLS solution collapses at z_c=0.68L_DF, see Fig. 2(C). This comparison of the NLH and NLS provides a direct numerical evidence that collapse is arrested in the scalar NLH model.

The fast oscillations of |E|² in the z direction in Fig. 2(A) are not a numerical artifact, but rather account for the actual physics. Indeed, let us first note that a part of the forward-propagating wave is reflected backwards by the material interfaces at z=0 and z=Z _max. In addition, since the forward propagating beam induces changes in the refraction index, part of the beam may be backscattered inside the Kerr medium. The presence of both forward and backward traveling fields, i.e,

implies that ${\mid E\mid }^2\approx {\mid A\mid }^2+{\mid B\mid }^2+2\mathrm{Re}\left({\mathrm{AB}}^{*}{e}^{i2k_{0}z}\right)$ . Hence, |E|² should undergo oscillations with wavenumber ~2k ₀. Note that the analytical solutions of the (0+1)D NLH also exhibit these 2k ₀ intensity oscillations [16]. The prediction that the intensity undergoes 2k ₀ oscillations implies that the index of refraction also oscillates. In other words, the backward traveling field induces a 2k ₀ Bragg grating. This prediction may be tested by pump-probe experiments. In order to find a smoother representation of the solution, recall that for the NLS (3) the conserved beam power is P_NLS=∫|Ã|² d x̃ _⊥. For the NLH (1), however, the conserved beam power is P_NLH=∫S_z d x _⊥, where S=k ₀Im(E*∇E) is the energy flux, or Poynting vector, and $S_z=k_0\mathrm{Im}\left(E^{*}\frac{\partial E}{\partial z}\right)$ is its z-component. Specifically, for the field (4) the value of S_z reduces to the flux difference S_z≈k ² ₀ (|A|²-|B|²). It is therefore much smoother than |E|², and provides a “more natural” depiction of the NLH solution, as confirmed by comparing S_z of Fig. 2(B) with |E|² of Fig. 2(A). The energy flux S_z shows the arrest of collapse and the focusing-defocusing dynamics more clearly, see also Fig. 2(C).

 

Fig. 2. (color online) Arrest of collapse in the (2+1)D NLH. A: |E|². B: S_z. C: comparison of normalized on-axis |E|² (blue solid), S_z (red dashed), and NLS solution (black dotted)

Download Full Size | PDF

 

Fig. 3. (color online) NLH solutions with $\frac{r_0}{{\lambda }_0}=\frac{3}{\pi }$ (blue, dots), $\frac{4}{\pi }$ (red, dash) and $\frac{6}{\pi }$ (green, solid). Solid black line is the NLS solution. A: Normalized on-axis intensity |E/E(z=0)|². B: Normalized on-axis Poynting flux S_z/S_z(z=0). C: Transverse profile of the backward field at z=0-.

Download Full Size | PDF

In order to analyze the effect of the nonparaxiality parameter f ², in Fig. 3 we fix the wavelength and vary the input beam radius r ₀ (while keeping the power unchanged) so that f ^-2=36, 64, and 144. All the NLH solutions initially follow the collapsing NLS solution, but later the collapse is arrested and the solution defocuses. As expected, for a wider input beam (lower nonparaxiality), the deviations from the NLS solution and the arrest of collapse occur later, and the maximum self-focusing is higher. Again we see that |E|² has 2k ₀ oscillations (whose magnitude increases as the input beam becomes more nonparaxial), while the energy flux S_z is smooth.

Our numerical algorithm for solving the NLH also enables the computation of backscattering from the Kerr medium. In Fig. 3(C), we present the backward propagating field profiles for the previous three solutions, just before the material interface at z=0. To the best of our knowledge, this is the first ever calculation of the backscattered field of collapsing beams, which is due to backscattering from inside the Kerr medium and reflections from the nonlinear interface. As the input beam radius r ₀ decreases, the power of the backscattered field increases from 0.46% to 0.63%, to 2.1% of the incoming beam power. This, as well as a comparison of magnitudes of oscillations in Figure 3(A), shows that the backscattered field increases as the input beam becomes more nonparaxial.

 

Fig. 4. (1+1)D NLH soliton with $r_0=\frac{{\lambda }_0}{2}$ propagating over 40L_DF. A: S_z. B: On-axis |E|².

Download Full Size | PDF

In (1+1)D configurations, the NLS possesses stable soliton solutions. It is generally believed that the paraxial approximation breaks down when the beam width becomes comparable to λ₀ and that, therefore, no solitons of such narrow width exist. To see that this is not the case, the (1+1)D NLH (1) is solved for the incoming NLS-soliton profile $E_{\mathrm{inc}}={\left(\sqrt{\frac{{2n}_2}{n_0}}{r}_0{k}_0\right)}^{-1}$ sech (x/r ₀) with width r ₀=λ₀/2, impinging on a Kerr slab of finite length Z _max=40L_DF. As in the (2+1)D case, we impose continuity of E and E_z at the material interfaces z=0 and z=Z _max, and apply the radiation BCs in the linear regions at z=-δ and z=Z _max +δ. The solution inside the Kerr-slab resembles a “nonparaxial soliton” which propagates virtually unchanged, see Fig. 4(A). We note that even for such a narrow beam, the nonparaxiality parameter is still moderate, as f ²=1/π ²≈0.1, which may explain why there there still exist soliton-like solutions. Similarly to the (2+1)D case, because a part of the forward propagating beam is backscattered, |E|² exhibits the fast 2k ₀ oscillations (Fig. 4(B)), while S_z is smooth. In this case, the backscattered field leads to 10% oscillations in |E|².

Posada, McDonald and New [17, 18] studied solutions of the (1+1)D Helmholtz equation over a semi-infinite Kerr medium, of the form $A(z,x)={\left(\sqrt{\frac{{2n}_2}{n_0}}{r}_0{k}_0\right)}^{-1}\mathrm{sech}\left(\frac{x}{r_0}\right)·e^{ic·z}$ . In later works they found similar stationary states for different nonlinearities [19, 20]. These solutions do not have any backward propagating components. In contrast, for the finite-length Kerr medium simulation of Fig 4, some backward moving waves must exist, because of reflections from the material discontinuity at z=Z _max, and the full NLH as a boundary-value problem must be solved.

Another (1+1)D configuration of recent interest is that of counter-propagating beams, when a right traveling soliton impinges at the left interface and a left traveling beam impinges at the right interface (Fig. 1(B)). This configuration was analyzed numerically by Cohen et al.[5] using a coupled NLS system, which is derived from the NLH by employing the paraxial approximation and further assuming that asynchronous terms of the Kerr nonlinearity can be neglected. In doing so, the BCs should simultaneously account for the coupled incoming and outgoing fields at each interface. As noted in [5], these BCs can only be approximately accommodated in the coupled NLS model. In contrast, they can be fully implemented in the NLH model, without any approximation. Figure 5 presents our solution of the NLH for counter-propagating beams of radius r ₀=λ ₀ that enter a Kerr material slab at the opposite interfaces with a transverse displacement of d=4.4λ ₀, and propagate over 10L_DF. It shows that the beams are slightly attracted toward each other and also become wider as they propagate. The results are in close agreement with the coupled NLS model, see Fig. 5(B). Therefore, the more comprehensive NLH model confirms the validity of the coupled NLS model for counter-propagating beams even for the “extreme” parameters of r ₀=λ ₀ and d=4.4λ ₀.

In this work, we solve the scalar NLH (1), which is the simplest model for the propagation of light in Kerr media that incorporates nonparaxial effects, backscattering and reflections from material interfaces. Moreover, unlike the NLS, it accurately models the propagation of oblique beams and reflections from interfaces at arbitrary angles. This model neglects vectorial effects, i.e., linear and nonlinear coupling between the three components of the electric field. We note that the vectorial effects scale as f ², and hence are of the same order of magnitude as nonparaxiality. In bulk media, they have been shown to have the same effect as nonparaxiality, which is to arrest the collapse [21, 6]. In contrast, in planar waveguides, the solitons are stable. Hence, when r ₀=λ ₀/2, f ² is small, and so we expect that both nonparaxial and vectorial effects are likely to have a secondary effect on the propagation dynamics. Therefore, the sub-wavelength solitons predicted in the Letter are likely to remain stable also in the more comprehensive vector NLH model.

 

Fig. 5. (color online) Energy flux (S_z) of the (1+1)D NLH with counter-propagating beams. A: Positive (forward) flux is red, negative (backward) flux is blue. B: The right-going beam at its incoming (blue dashed) and outgoing interface (red solid). Green dotted line is the coupled-NLS solution at the outgoing interface.

Download Full Size | PDF