Numerical simulations of self-focusing of ultrafast laser pulses Gadi Fibich* School of Mathematical Sciences, Tel Aviv University, Tel Aviv 69978, Israel Weiqing Ren Courant Institute of Mathematical Science, New York University, New York, New York 10012 Xiao-Ping Wang Department of Mathematics, The Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong Received 14 November 2002; published 7 May 2003 Simulation of nonlinear propagation of intense ultrafast laser pulses is a hard problem, because of the steep spatial gradients and the temporal shocks that form during the propagation. In this study we adapt the iterative grid distribution method of Ren and Wang J. Comput. Phys. 159, 246 2000 to solve the two-dimensional nonlinear Schro ¨dinger equation with normal time dispersion, space-time focusing, and self-steepening. Our simulations show that, after the asymmetric temporal pulse splitting, the rear peak self-focuses faster than the front one. As a result, the collapse of the rear peak is arrested before that of the front peak. Unlike what has sometimes been conjectured, however, collapse of the two peaks is not arrested through multiple splittings, but rather through temporal dispersion. DOI: 10.1103/PhysRevE.67.056603 PACS numbers: 42.25.Bs, 42.65.Sf, 42.65.Jx I. INTRODUCTION The nonlinear Schro ¨ dinger equation NLS i z z , x , y + +| | 2 =0 1 is the model equation for the propagation of cw continuous wavelaser beams in Kerr media. Here, is the electric field envelope, z is the axial distance in the direction of beam propagation, x and y are the coordinates in the transverse plane, and = xx + yy is the diffraction term. It is well known that when the power, or L 2 norm, of the input beam is sufficiently high, solutions of Eq. 1can self-focus and be- come singular in a finite distance z 1,2. Because of the infinitely large gradients that exist at the singularity, standard numerical methods break down after the solution undergoes moderate focusing 3. Therefore, as part of the research ef- fort during the 1980s to find the blowup rate of the NLS, McLaughlin et al. developed the numerical method of dy- namical rescaling 4, which can resolve the solution near the singularity at extremely high amplitudes. This method exploits the known self-similar structure of the collapsing part of the solution near the singularity, which relates the shrinking transverse width of the solution to the increase in its norm. Therefore, the solution is computed on a fixed com- putational grid, which in physical space corresponds to a grid that shrinks uniformly toward the singularity. The focusing rate of the grid points is determined dynamically from some norm of the solution ( | | | 2 dxdy , max x,y ||, etc.. Be- cause the focusing rate of the grid points can be chosen to be the same as the physical focusing rate, in the transformed variables the solution remains smooth and can thus be solved using ‘‘standard’’ methods. The method of dynamic rescaling works extremely well for solutions of the NLS with radially symmetric initial con- ditions, in which case focusing by 10 10 or more can easily be realized see, e.g., Fig. 3.5 in 1. Although the method of dynamic rescaling has been extended to NLS’s with noniso- tropic initial conditions 5and to perturbed NLS’s e.g., NLS’s with normal time dispersion 6, in such cases dy- namic rescaling is considerably less efficient, because the solution does not focus uniformly and/or it is not clear how to extract the physical focusing rate from the solution. The iterative grid redistribution IGRmethod, developed by Ren and Wang, overcomes these difficulties by allowing the grid points to move independently rather than uniformlyaccord- ing to a general variational minimization principle. This method has been showed to be highly effective for solving partial differential equations PDE’swith singular behavior such as the NLS 1and the Keller-Segal equations with multiple blowup points 7. As we shall see, however, apply- ing the IGR method to nonstationary NLS models that de- scribe the propagation of ultrashort pulses turns out to be considerably more demanding. A. Simulations of ultrashort pulses The NLS 1does not include temporal effects. These effects become important in the case of ultrashort laser pulses, whose propagation can be modeled by the dimension- less nonstationary NLS 8 i z z , x , y , t + +| | 2 + 1 zz +i 2  | | 2 t - t - 3 tt =0, 2 where now is also a function of time t. The dimensionless parameters are given by *Electronic address: fibich@math.tau.ac.il PHYSICAL REVIEW E 67, 056603 2003 1063-651X/2003/675/0566039/$20.00 ©2003 The American Physical Society 67 056603-1