Numerical modeling considerations for an applied nonlinear Schrödinger equation Todd A. Pitts,* Mark R. Laine, Jens Schwarz, Patrick K. Rambo, Brenna M. Hautzenroeder, and David B. Karelitz Sandia National Laboratories, P.O. Box 5800, MS-0980, Albuquerque, New Mexico 87185, USA *Corresponding author: [email protected] Received 3 October 2014; revised 6 December 2014; accepted 31 December 2014; posted 6 January 2015 (Doc. ID 221360); published 17 February 2015
A model for nonlinear optical propagation is cast into a split-step numerical framework via a variable stencil-size Crank–Nicolson finite-difference method for the linear step and a choice of two different nonlinear integration schemes for the nonlinear step. The model includes Kerr, Raman scattering, and ionization effects (as well as linear and nonlinear shock, diffraction, and dispersion). We demonstrate the practical importance of numerical effects when interpreting computational studies of high-intensity optical pulse propagation in physical materials. Examples demonstrate the significant error that can arise in discrete, limited precision implementations as one attempts to improve practical operator accuracy through increased operator support size and sampling frequency. We also demonstrate the effect of the method used to obtain the finite-difference operator coefficients defining the equations ultimately used in the discrete model. Smooth, plausible, but incorrect solutions may result from these numerical effects. This implies the necessity of a complete, precise description of all numerical methods when reporting results of computational physics investigations in order to ensure proper interpretation and reproducibility. © 2015 Optical Society of America OCIS codes: (190.5530) Pulse propagation and temporal solitons; (190.7110) Ultrafast nonlinear optics; (350.5500) Propagation. http://dx.doi.org/10.1364/AO.54.001426
1. Introduction
Mathematical modeling of intense optical pulse propagation physics via the nonlinear Schrödinger equation has applications in an increasing number of fields. Optical parametric amplifiers and oscillators [1,2], second-harmonic generation [3] and third harmonic generation [4], Raman scattering [5], optical bistability [6], solitons [7–9], and many other phenomena have found numerous practical applications. In support of these applications, as well as investigations into the nature of the physical phenomena themselves, the scientific and engineering communities expend significant effort in numerical modeling [10–13]. It is important to realize that a numerical 1559-128X/15/061426-10$15.00/0 © 2015 Optical Society of America 1426
APPLIED OPTICS / Vol. 54, No. 6 / 20 February 2015
model tends to introduce its own contributions to the intended physics model, altering the description of pulse evolution. Numerical physics investigations must therefore attempt to understand and separate the numerical from the physical. The numerical examples in this paper demonstrate the importance of completely describing not only the physical but also the numerical model when reporting on such investigations. The literature on nonlinear propagation is very broad [13]. A variety of models and unit systems are used to describe various subsets of known physical effects (see Brabec and Krausz [14], Couairon and Mysyrowicz [13], or Zozulya et al. [15]). We summarize the complete propagation and material interaction model used in our work in Section 2. This model is derived in detail in [16], bringing all physical phenomena of interest under the same framework and
clarifying the necessary physical assumptions and approximations. The primary differences between the work described in this paper and that in well-known canonical references such as Agrawal [17] and Dudley and Taylor [18] are as follows. The physical problems we are interested in (free-space propagation) require a diffraction model represented by the transverse Laplacian in Eq. (1) of the model summary in Section 2. This operator necessitates additional numerical methods discussed in Section 3.A. Also, in addition to the Kerr effect and Raman scattering we have models for ionization [see Eqs. (8–14) in Section 2]. Many frameworks without ionization use a Taylor series expansion for the Raman response in Eqs. (6) and (7) (see Agrawal [17]), which simplifies the nonlinear propagation step. This approach is related to the simple integration described in Eq. (42). When combining both Raman and ionization effects, this method is problematic. For this reason we use the more general approach in Eq. (43), which allows inclusion of very general material interaction models. In Section 3 we develop a numerical model for the physics. We describe the overarching numerical method (split-step), which breaks the problem into linear and nonlinear operators applied successively over short intervals along the propagation coordinate axis. Split-step has a rigorous underpinning despite an often heuristic description. It can be implemented with a reasonable amount of effort. Perhaps most importantly for physicists and engineers, it is straightforward to modify existing or add new physical submodels that would otherwise be unwieldy. For the linear operator, the dispersion and linear shock are evaluated via Fourier transform in the frequency domain. This method is fast and provides a very accurate implementation of the continuous temporal operators. However, it does require regular sample spacing along the temporal axis and that we keep the pulse away from the initial and final grid times to prevent a physical scattering from the periodic boundary conditions. The transverse Laplacian (diffraction) is implemented via finite differences using a Crank–Nicolson (CN) method. This requires the solution of a series of large linear systems (one per frequency in the temporal FFT mentioned above). CN is an implicit, unconditionally stable method that typically admits larger steps than explicit methods (which must satisfy the Courant–Friedrichs–Levy criterion). It is nondissipative and therefore does not introduce artificial loss into the system. Its principle drawback for our application is that it is dispersive and can sometimes suffer from oscillations. We derive all these operators in Section 3.A.1 and set up the computational framework in Section 3.A.2 to easily vary the size of the operator support. We also discuss the use of a preconditioner and its effect on the accuracy and speed of the solution for the large linear solves in Section 3.A.3.
Section 3.B describes methods for solving the nonlinear operator portion of the split-step framework. For our physical model there are actually two equations that must be integrated in this step: a description of the evolution of electron number density and the nonlinear portion of the split-step method. Both may be solved via Runge–Kutta methods. For problems in which the nonlinear portion of the splitstep framework takes a particularly simple form, we show that a simple nonlinear interest model suffices. In each section, we give numerical examples showing the effect of different choices for various numerical model parameters, such as operator support size, sample spacing, z-direction step size, and use of a preconditioner. The results are summarized in Section 4. 2. Summary of the Physical Model
Here we summarize the physical model that we discretize and solve in Section 3. It includes paraxial diffraction, group velocity and third-order dispersion, linear and nonlinear shock, Kerr and Raman effects, and ionization. There are many models for nonlinear optical propagation in the literature [13,14,19–25]. Ours is similar to those in the aforementioned references. However, our method of including the impedance effects in the plasma defocusing term leads to a difference in the shock terms (see, for example, [13]). Reference [16] derives the model described here in detail with a clear explanation of all assumptions and approximations. The propagation equation for the electric field envelope E is ∂ζ E
i −1 2 T ∇⊥ E iDE TNE; 2k
(1)
with wavenumber k 2π∕λ and wavelength λ. The dispersion and shock operators are given by k00 k000 D ≡ − ∂2τ − i ∂3τ 2 6
and
i T ≡ 1 ∂τ ; ωo
(2)
respectively. The nonlinear operator may be decomposed as N N kr T −1 N pls T −1 N MPA :
(3)
Kerr and Raman effects are modeled with N kr iko n2 1 − αI αQR :
(4)
We have used the approximately time-harmonic relationship I≈
cno εo 2 jEj 2
(5)
to specify the intensity envelope I. Delayed Raman response is described by 20 February 2015 / Vol. 54, No. 6 / APPLIED OPTICS
1427
Z QR r; t; z
t −∞
RτIr; t − τ; zdτ;
(6)
(13)
If we do not wish to include ionization saturation effects, we may instead use
with Rt
γ 2 ω2R exp−γt sinωR t; ωR
(7)
where ωR and γ are, respectively, characteristic frequency and attenuation constants for the molecular species interacting with the pulse. The following quantities help us express the remaining ionization physics: τc Electron collision or relaxation time ωo Pulse center frequency no Background refractive index U i Energy required to free single electron ρnt Neutral atom number density αr Electron-positive ion recombination loss σ k Ionization cross-section: For plasma defocusing we use (8)
or, if we desire to include collision effects, (9)
with ρc ≡ εo me ω2o ∕e2 ;
and σ ≡
ko ωo τc : no ρc 1 ω2o τ2c
(10)
The energy loss in the pulse due to multiphoton absorption is included via N MPA −
βK K−1 I ; 2
(11)
with βK σ k ρnt Kℏωo :
We seek to solve an equation of the form ∂z A LA NAA:
(15)
u0 uo ;
(16)
where L may be a partial or ordinary differential operator. Use of the Baker–Campbell–Hausdorff formula leads to the development of higher order stepping patterns, but we consider only the firstand second-order frameworks here. For these, we solve the two equations
K Photons required to free single electron
σ N pls ≡ − 1 iωo τc ρ; 2
(14)
3. Numerical Model
ut Lu;
ko Pulse center wavenumber
iko T −1 ρE; N pls E − 2no ρc
dρ σ β 2 ρI K I K − αr ρ2 : Kℏωo dt no U i
The split-step approach (see Strang [26]) can be derived rigorously by considering an initial value problem of the form
me Mass of an electron
(12)
The evolution of the electron density requires an auxiliary equation: 1428
dρ σ 2 ρI I K σ k ρnt − ρ − αr ρ2 : dt no U i
APPLIED OPTICS / Vol. 54, No. 6 / 20 February 2015
∂z A LA;
(17)
∂z A NAA
(18)
in an alternating fashion. Solving each of these equations will require additional numerical methods. The numerical evolution of the pulse may take place along either the temporal or z coordinate (where z is the principal propagation direction). We describe the evolution of the pulse through a plane over time at successively increasing z locations. The numerical grid has transverse spatial coordinates r or x; y and a temporal coordinate. At any given step in the propagation algorithm, this grid exists at a particular location in the primary propagation direction and describes the spatial and temporal variations of the pulse over the transverse plane. Each iteration of the propagation algorithm transforms the transverse space–time description at the current location to the distribution at the next location in the propagation direction. It is a z-propagated description of the pulse evolution. We begin with a description of the initial pulse and solve Eq. (17) over a short distance (z step) along the direction of propagation. In other words, we linearly propagate the pulse a distance dz. Section 3.A describes the linear numerical model used to accomplish this task. We use the resulting pulse description as an initial condition for Eq. (18) as described in Section 3.B. This completes a single first-order z-direction step. The process of alternating solution of the linear and nonlinear operator
equations continues until the desired number of zdirection steps have been taken. A second-order solution consists of taking half a linear step, an entire nonlinear step, and then another linear half-step. Diffraction and self-focusing can balance one another. However, in radially symmetric models where only a single filament is present there is a limit on the amount of power in the pulse that diffraction can balance. The same is true of free-electron defocusing in problems that involve ionization. Eventually, either fully (3D 1) models must be used so that multiple filaments are allowed or the pulse will collapse. A.
Linear Operator
Our basic framework for solving the linear equation could take many forms. The complex exponential is an eigenfunction of linear shift-invariant (LSI) operators. This fact leads to the prevalence of spectral methods in the solution of LSI partial differential equations in circumstances where the mathematically implied and physical boundary conditions are compatible. The implied periodicity of the discrete Fourier transform leads to scattering off the edges of the domain. Waves hitting the boundary have corresponding sources from an adjacent domain radiating energy into our numerical grid. In many situations (converging beams, etc.) this would be acceptable. However, even relatively low amplitude disturbances near the edge of the pulse may cause problematic numerical disturbances if they propagate to the transverse edges of the domain and reflect back to the center. For these reasons we use a Fourier-transformbased method in time only. We transform the linear operator into the temporal frequency domain leaving only the transverse derivatives. We then solve the resulting linear algebraic problem for each temporal frequency in the pulse. This allows extension of the method to include absorbing boundary conditions and adaptive mesh refinement (see Mlejnek et al. [10] for the first reported simulations of nonlinear propagation via AMR). 1. Basic Differential Operators We begin with an approximation to the nth derivative at location x taken over a support of q samples. We may approximate the derivative [27] q X dn u ≈ γ u γ u γ u γ i ui : 1 1 2 2 q q dxn x i1
(19)
To obtain the necessary coefficients in Eq. (19) for an odd-length stencil of size q 2p 1, with a derivative of order n, evaluated at the center stencil location we may solve PΓ l; where
(20)
q−1 h Phk k − ; 2
2
Γ1 Γ2 Γ3 .. .
3
7 6 7 6 7 6 7 6 7; 6 Γ6 7 7 6 6Γ 7 4 q−1 5 Γq
h; k 0; …; q − 1;
(21)
2
3 0 . 6 .. 7 6 7 607 6 7 6 7 and l 6 n! 7; 6 7 607 6 . 7 4 .. 5 0
(22)
and the only non-zero entry n! on the right-hand side is in the m th position (the initial position is designated zero). We then form γ Δ−n Γ;
(23)
where Δ is the physical spacing between samples. By using Eqs. (20)–(23) we obtain the finite-difference coefficients for a given derivative order and sufficiently broad spatial support. Figure 1 shows plots of the error obtained by using finite-difference operators with coefficients obtained in this fashion via three different linear algebraic solution methods. We use a simple unit period sinusoid as our test function, computing the derivative at x 1∕4 via finite differences. We then compute ϵ, the absolute value of the difference between the finitedifference estimate and the analytical derivative 2π cos2πt at the same location. Each plot shows log10 ε as a function of the operator support and log10 dt, where dt is the finite-difference sample spacing. Operator coefficients obtained via Eqs. (21) and (22) are consistent approximations, meaning that as the finite-difference spacing goes to zero the approximations must improve [27]. Initially, we see this trend in all three plots. As we move from the smallest possible odd-length support and the largest shown sample spacing toward longer filter supports and smaller sample spacing, the accuracy of the operator improves. However, in both directions we eventually see that numerical inaccuracies in the coefficients lead to a reduction in the accuracy of the operator. This means that while the finite-difference operator is technically consistent, practical considerations such as round-off error caused by finite precision arithmetic do not always allow us to obtain improved accuracy by decreasing the sample spacing or increasing the support of the filter. As we move between the plots we see that the accuracy of the numerical method used to obtain the coefficients for the finite-difference operator is important. Using a numerically stable method gives clearly better coefficients. This is particularly important in AMR frameworks where variable grid spacing requires us to obtain coefficients for many different sampling patterns. 20 February 2015 / Vol. 54, No. 6 / APPLIED OPTICS
1429
−2 −4 −6
−2
−10
−2 0 −6−4
−8
−3
−12
−8
−2 −4 −6
−8
−10
−4 −8 log (dt)
−8
−2 −4 −6
LU Decomposition with Partial Pivoting −8−10 −10
−6
−8
−12
−8−10 −12
10
−5 −10
−6
−4
−7
−4 −5
−10 −10
−10
−2 −3
−12
−10
−12 −12
−1
−6 −10
−6
log10(dt)
Gaussian Elimination
Moore−Penrose Pseudo−Inverse −1
−7
−8 −8
−8
−6
−9
−2
−4
−8
−2 −4
−8 −9
0 −6
−10
3
5
7 9 11 13 Operator Support (samples)
15
3
5
7 9 11 13 Operator Support (samples)
(a)
3
15
5
7 9 11 13 Operator Support (samples)
(b)
15
−10
(c)
Fig. 1. Finite-difference coefficient and operator accuracies. Each contour plot shows the base 10 logarithm of the absolute value of the difference between the finite-difference operator estimate and the analytical first derivative of the unit period sinusoid sin2πx at x 1∕4. (a) Finite-difference coefficients obtained by solving Eqs. (20)–(22) via a Moore–Penrose pseudo-inverse. (b) Coefficients obtained via Gaussian elimination. (c) Coefficients obtained via lower–upper (LU) decomposition with partial pivoting.
Significantly reducing the radial sample spacing and z-step size (magenta curve) alters the solution only marginally. If we push the operator support and number of radial samples too high the solution begins to degrade (green curve). Any further increases in either operator support or the number of radial samples causes the numerical solution to break down entirely (likely due to round-off error).
Figure 2 shows the effects of varying operator support, radial sample spacing, and the number of z steps used to compute a solution. Results for numerical and analytical linear propagation of a focused Gaussian beam are shown in Fig. 2(a). We observe that with three-point operator support and 64 radial samples that the quality of the numerical solution is low. We may increase either operator support, the number of radial samples, or both in order to obtain a better solution. In Fig. 2(b) we see similar effects for nonlinear propagation. The red curve has minimal operator support, radial sampling, and a small number of z steps. Increasing only the operator support changes the solution dramatically without requiring significantly more computer resources (blue curve).
2. Discretization of Propagation Equation Here we describe the discretization of the linear spatial operator including the application of relevant boundary conditions. We do so in a fashion that allows the programmatic specification of differential Stencil size: 3, Radial samples: 90, z−steps: 128 Stencil size: 13, Radial samples: 90, z−steps: 128 Stencil size: 13, Radial samples: 256, z−steps: 256 Stencil size: 17, Radial samples: 700, z−steps: 256
Intensity (normalized)
Amplitude (Normalized)
3
2
Stencil size: 3, Radial samples: 64, z−steps: 200 Stencil size: 3, Radial samples: 256, z−steps: 200 Stencil size: 13, Radial samples: 64, z−steps: 200 Stencil size: 13, Radial samples: 256, z−steps: 200 Analytical solution
1
0
0
1
2
3
3
2
1
4
0
−150
−100
−50
0
Distance (m)
Time (fs)
(a)
(b)
50
100
150
Fig. 2. Variable stencil size. (a) Linear propagation of a focused Gaussian beam with focal distance f 3 m, e−1 beam radius re 1.5 mm, and wavelength λ 810 nm. The dotted line shows the analytical solution. Each numerical solution has identical z-step spacing. It is clear from the plots that a significant increase in accuracy may be achieved either by decreasing the radial sample spacing or increasing the operator support. Increasing the operator support requires much fewer computer resources. (b) Nonlinear propagation with linear and nonlinear shock using a linear algebraic preconditioner for solution of finite-difference equations. For this model the linear refractive index of the material is no 1.45, nonlinear refractive index n2 2.5 × 10−20 m2 ∕W, the group velocity dispersion is 3.6 × 10−26 s2 ∕m, the thirdorder dispersion coefficient is zero, the coefficient controlling the percentage of instantaneous Kerr to Raman effect is γ 0.15, and the characteristic Raman response time and frequency are, respectively, τR 50 × 10−15 s and ωR 84 × 1012 Hz. The beam is propagated 3 cm. The curves show the effect of moving from minimal operator support, number of radial samples, and number of z steps through increases in each. The blue and magenta curves are computed with the most accurate operators. Further increases in radial samples or operator support cause the solution to break down very quickly. 1430
APPLIED OPTICS / Vol. 54, No. 6 / 20 February 2015
operator order and, hence, stencil size. This provides a means of developing an empirical understanding of how operator order affects solution accuracy when combined with sample spacing in the presence of round-off error. In the frequency domain the linear portion of Eq. (1) takes the form ∂z u iq1 u iq2 ∇2⊥ u;
(24)
with q1 and q2 real. Define the operator dm n to be an nth order derivative in space (in our case the rˆ direction) with support of size m. For example, d31
1 −1; 0; 1: 2Δ
(26)
In cylindrical coordinates the Laplacian is 1 ∇2⊥ ≡ ∂2r ∂r : r
(27)
The CN approximations are then 1 u → un1;j un;j ; 2
∂2r u →
(28)
1 u − un;j ; δz n1;j
(29)
1 m d · un1;jm un;jm ; 2 2
(30)
∂z u →
1 1 ∂ u → dm · un1;jm un;jm ; r r 2r 1
(31)
where the notation 2
un;jm
3 un;j−m−1∕2 6 un;j−m−1∕21 7 6 7 .. 6 7 7 6 . 7 6 7 u 6 n;j 7 6 7 6 . .. 7 6 7 6 4u 5 n;jm−1∕2−1
un;j−m−1∕2
m pm j · un1;jm pj · un;jm ;
(33)
where pm j ≡
2 m 1 m m d d : d0 − i q1 dm q 2 0 2 δz r 1
(34)
Note that at j 0 we must use L’Hospital’s rule to arrive at
(25)
The operator in Eq. (25) includes the grid spacing Δ. In general, the operator is inversely proportional to the n th power of the grid spacing for an order n derivative. Section 3.A.1 describes how we arrive at the coefficients. Without loss of generality we concern ourselves only with odd support lengths. By definition, the dm 0 is a support centered Kronecker delta function: d50 0; 0; 1; 0; 0:
describes an m-length vector evaluated at step n in the z-coordinate direction, centered at grid node index j in the r-coordinate direction. Substituting into Eq. (24) gives (after some algebra)
pm 0 ≡
2 m m d − iq1 dm 0 2q2 d2 : δz 0
Equation (33) is an algebraic equation in m unknowns. We know the un; jm (they are the field values at the particular frequency specified by the values of q1 and q2 as a function of radial coordinate at z-step index n) and the entries of the row vectors pm j . The right-hand side of Eq. (33) is a scalar. The m unknowns are un1; jm (the field values at the particular frequency specified by the values of q1 and q2 as a function of radial coordinate at z-step index n). If we add another equation centered at an adjacent value of j (the radial coordinate index), then we have an additional equation and an additional unknown (two equations in six unknowns). Each time we add an adjacent equation we get an additional variable. If we add an adjacent equation until we reach a boundary, the new variable added by each boundary equation is known. At each boundary we add an equation with a known additional variable m − 1∕2 times. After reaching both boundaries we have exactly the number of equations we need to solve the system. Solving the system lets us update all the un;j to un1;j for all j at once. This method (CN) is implicit (uses both the current and future time to obtain the future time) and, hence, requires a solve (in this case a linear solve) at each step. It is unconditionally stable and will tolerate a larger step size than explicit methods (which are stable only if the step size is chosen according to the Courant stability criterion). The price for this unconditional stability is a global linear algebraic solve at each step. Practical implementation for a grid with J radial samples requires that we use Eq. (34) and the known current values of un to create a matrix equation: Lun1 L un ;
32
(35)
(36)
where L is a matrix of dimensions J × J and un is a known vector of length J. Row j of L is formed from Eq. (34). The top and bottom m − 1∕2 rows are special as they touch the boundaries. If we apply symmetric boundary conditions at the j 0 (r 0 20 February 2015 / Vol. 54, No. 6 / APPLIED OPTICS
1431
boundary) and a Dirichlet boundary condition with a value of zero at the outer boundary, noting that −l −l1 …p0j …pjl−1 plj ; pm j pj pj
(37)
where l m − 1∕2, then we may form the augmented matrix 2 6 6 6
p−l 0
p0−l1 p−l 1
0 … p0 p1−l1 …
… p01
p0l−1 … .. .
pl0 p1l−1
pl1 ; (38a)
.. p−l J−2
−l1 pJ−2 p−l J−1
. …
−l1 pJ−1
p0J−2 …
p1J−2 p0J−1
7 7 7 7 5;
38b
where Eq. (38a) describes the upper-left portion of the matrix and Eq. (38b) shows the lower right. The vertical line in Eq. (38a) represents the central propagation axis of our numerical grid. We must apply a symmetry boundary condition here. Mathematically, this implies that we fold the entries to the left of this line over onto those on the right-hand side and add. Equation (38b) shows a Dirichlet boundary condition with a value of zero on the right, meaning we do not need to extend beyond the matrix edge. We simply truncate when the center of the pm j reaches the edge (this will, of course, always be at row J − 1 when we start with index zero). With this understanding we can easily vary the stencil size via a parameter specified at run time. 3. Linear Solvers and Preconditioners Application of the finite-difference method generates large systems of linear equations for each frequency. Solving these systems implements the forward linear propagation at each frequency. Often the structure of the matrix is exploited to obtain a stable solution. For very large matrices, parallel computation may provide a means of propagating pulses requiring very large numerical grids. In these cases it is difficult to exploit the matrix structure. Instead, Krylov methods are used [28,29]. These methods find the solution using matrix-vector products only. These can be formed efficiently on distributed systems. For these examples, we use the biconjugate gradient stabilized (BICGStab) method [30]. In doing so we have found that preconditioning is necessary to obtain an accurate solution. The condition number of a square matrix is the magnitude of the ratio of the largest to the smallest eigenvalue. It is a measure of the sensitivity of the solution to a linear system of equations to errors in the data (which may result from the limited precision of the arithmetic). The purpose of preconditioning is to obtain a system of equations whose condition number is near unity (ideal case). This means that all of the eigenvalues 1432
are clustered near one another, which improves the accuracy and convergence of Krylov space methods. Instead of solving
APPLIED OPTICS / Vol. 54, No. 6 / 20 February 2015
Ax b;
(39)
P−1 Ax P−1 b;
(40)
we solve
where P−1 A has a better condition number than A. In practice this means that the matrix P is in some sense close to a scalar multiple of the system matrix A. In our numerical simulations we use the system matrix as our preconditioner. It would seem that multiplication by the inverse of the preconditioner would solve the system, but due to the numerical inaccuracies there remains some error, which the iterative solver then reduces below a specified tolerance. Figure 3 shows the effect of using a preconditioner for the linear solves in the case of a beam with critical power propagating in a material with a pure Kerr nonlinearity. The mathematical instability of this scenario (no ionization) creates an opportunity for evaluating the accuracy of our numerical methods. Numerical errors in the simulation will have the same effect as very small refractive index perturbations in a physical experiment. The less accurate the simulation, the sooner the pulse will lose confinement. Each example propagates the pulse the same physical distance. Figure 3(a) shows that 1000 z steps provide sufficient numerical accuracy with or without a preconditioner. Both curves show the expected confinement. Figure 3(b) compares the results obtained when only 500 z steps are used. With a preconditioner (blue curve) the linear solver is sufficiently accurate to allow diffraction and Kerr self-focusing to balance each other effectively. Without a preconditioner (green curve), inaccuracies in the linear propagator accrue, causing the beam to immediately diverge after initially self-focusing. Figure 4(a) shows the effect of the linear preconditioner with both Kerr and Raman effects in the presence of linear and nonlinear shock. The plot shows the normalized axial envelope intensity after propagating 3 cm in fused silica (compare with Zozulya et al. [15]). It is clear that preconditioning the linear solves has a significant effect on the solution, leading to more accurate results with fewer z steps and a lower number of radial samples. The BICGStab method is also significantly faster when using a preconditioner. B. Nonlinear Operator
In this section we describe the implementation of the nonlinear material response operator. This includes instantaneous Kerr effects and Raman scattering. These two effects have different physical origins but are often combined in numerical models due to the fact that they are both proportional to the cube of the electric field value. The nonlinear material model also includes ionization described by Eq. (13). If we
1
Nz=1000, Nr=512, with PC Nz=1000, Nr=512, No PC Normalized pulse width
Normalized pulse width
1
0.9
0.8
0.7
0.6
0
50
100
150
200
250
300
Nz=500, Nr=1024, with PC Nz=500, Nr=1024, No PC
0.9
0.8
0.7
0.6
0
50
100
150
200
Propagation distance (m)
Propagation distance (m)
(a)
(b)
250
300
Fig. 3. This simulation shows the effect of using a linear preconditioner when the propagation takes place in a simple Kerr material (no Raman scattering or ionization). Both linear and nonlinear shock are turned off. The pulse-center wavelength is λ 1024 nm. The linear and nonlinear refractive indices are n 1 and n2 4 × 10−23 m2 ∕W, respectively. Group velocity dispersion is 2.6 × 10−29 s2 ∕m, and thirdorder dispersion is given as 1 × 10−85 s3 ∕m. The initial beam has e−1 beam radius re 6 mm, an initially flat wavefront (focus at infinity), power at 95% of the critical power, and a propagation distance of 300 m. All simulations used finite-difference operators with a support of 3 and a single temporal sample (beam). (a) The blue curve shows the result of a simulation with preconditioner, 512 radial samples over 5re, and 1000 z steps. The green is the result of an identical simulation with no preconditioner. (b) The blue curve shows the result of a simulation with preconditioner, 1024 radial samples over 5re, and 500 z steps [half as many steps per unit distance as in (a)]. The green curve is the result of an identical simulation with no preconditioner.
do not know the neutral atom density we may use instead Eq. (14). We also include a loss term in the nonlinear portion of the propagator (discussed in [16]) that is designed to account for energy loss due to multiphoton ionization by removing the requisite amount of energy from the pulse every step. The nonlinear portion of our split-step operator method is an equation of the form ∂ A NAA: ∂z
(41)
When the nonlinear operator in Eq. (41) can be written in the form where NA contains no differential operators applying to the final envelope expression A [in this case NA reduces to a number], then we have an easier problem. This does not mean there are no differential operators in NA. Sometimes there are shock terms that would normally operate on the final envelope A except for a judicious manipulation of the operator expression. When this simpler circumstance prevails, Eq. (41) is essentially a nonlinear first-order interest equation. For initial condition Az propagated over a distance h, the solution to this problem may be written Z Az h; r exp
zh z
NAdz0 Az
≈ exphNAAz:
(42)
When we wish to include the shock operators in a straightforward fashion (or cannot manipulate the expression to exclude application of differential operators to the final envelope), we have a differential equation of the form
d Az; t f Az; t; z; dz
z1 ≤ z ≤ z2 z1 h; (43)
with initial condition Az1 ; t given. We may integrate the equation with an explicit one-step Runge–Kutta(2, 3) method from Bogacki and Shampine [31], which is useful for moderately stiff ordinary differential equations of the form in Eq. (43). We must integrate the nonlinear ionization model in Eq. (13) [or Eq. (14) if the neutral atom density is not known] over the temporal duration of the pulse each time we perform the nonlinear step. We may use the same RK-(2, 3) method from and Shampine as above. The computations shown in this paper use this method. However, the solution of these equations can go through sharp transitions and regions of much slower evolution depending on the structure of the optical pulse. For this reason it is sometimes necessary to use ODE solvers designed for stiff differential equations (for example, those described in Shampine and Reichelt [32]). Notice that the equations for the nonlinear portion of the split-step method and the electron number density are coupled and should technically be solved simultaneously. We make the approximation that they may be solved sequentially, solving first for the evolution of number density using the output of the previous linear step and then using the result in the integration of the nonlinear splitstep equation. Critical power. There are many expressions for critical power depending on the model for transverse beam or pulse profile, the type of nonlinear material, and the criteria for entrapment (see, for example, Diels and Rudolph [33], Ashcom [34], Zozulya et al. [15], Marburger [35], or Wright et al. [36]). For 20 February 2015 / Vol. 54, No. 6 / APPLIED OPTICS
1433
Axial Envelope at 3 cm (λ = 800 nm)
3.5
Nz: 101, Nr: 512, No PC Nz: 111, Nr: 256, No PC Nz: 101, Nr: 512, PC Nz: 101, Nr: 256, PC Nz: 501, Nr: 512, PC
2.5
Linear, Nt: 64 Nonlinear, Nt: 32 Nonlinear, Nt: 64 Nonlinear, Nt: 512
8 Amplitude (Normalized)
Normalized axial intensity
3
Normalized Maximum Axial Field Amplitude
10
2 1.5 1
6
4
2 0.5 0 −200
−150
−100
−50
0 Time (fs)
50
100
200
150
(a)
0
0
0.01
0.02 0.03 Distance (m)
0.04
0.05
(b)
Fig. 4. (a) Simulation pulse-center wavelength λ 800 nm. Material parameters no 1.45, n2 2.5 × 10−20 m2 ∕W, group velocity dispersion 3.6 × 10−26 s2 ∕m, and no third-order dispersion. The coefficient controlling the percentage of instantaneous Kerr to Raman effect is γ 0.15, the characteristic Raman response time and frequency are, respectively, τR 50 × 10−15 s and ωR 84 × 1012 Hz. Both linear and nonlinear shock are turned on. The initial beam has e−1 beam radius re 59.45 μm, an initially flat wavefront (focus at infinity), a temporal FWHM of 90 fs, and power Pin 4.72 mW (approximately 1.83Pcr;Z ). The propagation distance is 3 cm. The spatiotemporal extent of the numerical grid is 4 times the temporal pulse duration and 4 times the radial extent. All simulations have 128 temporal samples. As indicated in the plot legend, the number of z steps, the number of radial samples, and the use of a preconditioner vary among the simulations. The two extreme examples (red and dashed red curves) do not use a preconditioner. All three remaining curves use a preconditioner during the linear, frequency-domain finite-difference update for each frequency and produce results very close to one another despite differences in step number and sample spacing. (b) This simulation shows the effects of varying the temporal sample rate. The pulse propagates 5 cm with no linear or nonlinear shock using 300 z steps. The linear solver uses the linear system matrix as its own preconditioner. The central wavelength is λ 586 nm, and the linear and nonlinear refractive indices are, respectively, no 1.00027 and n2 4.9 × 10−23 m2 ∕W. The group velocity dispersion is 2.6 × 10−26 s2 ∕m and there is no third-order dispersion. The initial beam has e−1 beam radius re 200 μm and temporal extent 200 fs. The grid has temporal and spatial extents of 1800 fs and 1800 μm, respectively. Each simulation uses 256 samples along the radial coordinate. The geometric focus is 2.5 cm and the input power is Pcr .
the Gaussian beam in the simulations reported in Figs. 3 and 4(b) we use Pcr
λ2 : 2πno n2
(44)
For the simulations in Fig. 4(a) we compute per Zozulya et al. [15]: Pcr;Z
π0.61λ2 : 8no n2
(45)
4. Summary and Conclusions
In this paper we developed mathematical and numerical physics models for intense optical pulse propagation in homogeneous, isotropic, nonlinear media. We included linear dispersion, diffraction, nonlinear Kerr and Raman effects, and ionization, as well as linear and nonlinear shock terms. The mathematical physics propagator is a one-way paraxial scalar model. The numerical split-step algorithm was explicitly described. We showed how the differential operators of varying support are derived for the linear step. Implementation of the temporal derivatives (including shock terms) takes place in the frequency domain via Fourier transform. The transverse Laplacian is modeled via a CN finitedifference method using a BICGStab solver with a linear preconditioner. Ionization models describing 1434
APPLIED OPTICS / Vol. 54, No. 6 / 20 February 2015
the evolution of electron number density are integrated via a well-known explicit, one-step Runge– Kutta solver. The nonlinear split-step integration uses either a simple, nonlinear interest model or the same Runge–Kutta method used for the ionization. We showed that understanding the effects of the numerical framework is very important when using computation to gain physical insight. As seen in Fig. 1, even in mathematically consistent finite-difference approximations, practical operator accuracy does not necessarily increase with operator support size and sampling frequency due to significant truncation and round-off error effects. Also, when programmatically solving for the finite-difference coefficients themselves, the solution method can have a marked effect on their accuracy and hence the subsequent trade-off between truncation and round-off error. In some cases, significant improvements in accuracy or reductions in run time may be obtained with larger stencils instead of decreasing sample spacing [see Figs. 2(a) and 2(b)]. Clearly demonstrated in these figures is the ease with which we might obtain a smooth, plausible, but incorrect solution. This also implies that some care must be taken when reporting results to precisely and completely describe the numerical methods to ensure meaningful interpretation and reproducibility. Figure 3 shows the effect of using a preconditioner compared with decreasing the step size in the direction of propagation. We may obtain accurate results with fewer steps. Also advantageous is the fact that the preconditioner generally
appears to significantly accelerate the linear solution in most cases. This manuscript has been authored by Sandia Corporation under Contract No. DE-AC04-94AL85000 with the U.S. Department of Energy.
17. 18. 19.
References 1. J. Hansryd, P. A. Andrekson, M. Westlund, J. Li, and P. Hedekvist, “Fiber-based optical parametric amplifiers and their applications,” IEEE J. Sel. Top. Quantum Electron. 8, 506–520 (2002). 2. G. Cerullo and S. De Silvestri, “Ultrafast optical parametric amplifiers,” Rev. Sci. Instrum. 74, 1–18 (2003). 3. P. A. Franken, A. E. Hill, C. W. Peters, and G. Weinreich, “Generation of optical harmonics,” Phys. Rev. Lett. 7, 118–119 (1961). 4. R. W. Terhune, P. D. Maker, and C. M. Savage, “Optical harmonic generation in calcite,” Phys. Rev. Lett. 8, 404–406 (1962). First observation of third-order harmonic generation. 5. E. J. Woodbury and W. K. Ng, “Ruby laser operation in the near IR,” Proc. IRE 50, 2365–2383 (1962). First observation of optical Raman scattering. 6. H. M. Gibbs, S. L. McCall, and T. N. C. Venkatesan, “Differential gain and bistability using a sodium-filled Fabry–Perot interferometer,” Phys. Rev. Lett. 36, 1135–1138 (1976). First demonstration and explanation of optical bistability. 7. A. Hasegawa and F. Tappert, “Transmission of stationary nonlinear optical pulses in dispersive dielectric fibers. I. Anomalous dispersion,” Appl. Phys. Lett. 23, 142–144 (1973). Theoretical prediction of soliton. 8. A. Hasegawa and F. Tappert, “Transmission of stationary nonlinear optical pulses in dispersive dielectric fibers. II. Normal dispersion,” Appl. Phys. Lett. 23, 171–172 (1973). Theoretical prediction of soliton. 9. L. F. Mollenauer, “The non-linear Schrödinger equation and ordinary solitons,” Unpublished presentation (Lucent Technologies, 2002). 10. M. Mlejnek, M. Kolesik, J. V. Moloney, and E. M. Wright, “Optically turbulent femtosecond light guide in air,” Phys. Rev. Lett. 83, 2938–2941 (1999). 11. A. Couairon, E. Brambilla, T. Corti, D. Majus, O. de J. Rámirez-Góngora, and M. Kolesik, “Practitioners guide to laser pulse propagation models and simulation,” Eur. Phys. J. Spec. Top. 199, 5–76 (2011). 12. L. Berge and A. Couairon, “Modeling the filamentation of ultra-short pulses in ionizing media,” Phys. Plasmas 7, 193–209 (2000). 13. A. Couairon and A. Mysyrowicz, “Femtosecond filamentation in transparent media,” Phys. Rep. 441, 47–189 (2007). 14. T. Brabec and F. Krausz, “Nonlinear optical pulse propagation in the single-cycle regime,” Phys. Rev. Lett. 78, 3283–3285 (1997). 15. A. A. Zozulya, S. A. Diddams, and T. S. Clement, “Investigations of nonlinear femtosecond pulse propagation with the inclusion of Raman, shock, and third-order phase effects,” Phys. Rev. A 58, 3303–3310 (1998). 16. T. A. Pitts, M. R. Laine, J. Schwarz, P. K. Rambo, B. M. Hautzenroeder, and D. B. Karelitz, “Derivation of an applied
20. 21.
22. 23.
24.
25.
26. 27. 28.
29. 30.
31. 32. 33. 34. 35. 36.
nonlinear schrödinger equation,” Tech. Rep. (Sandia National Laboratories, 2014). G. P. Agrawal, Nonlinear Fiber Optics, 4th ed., Optics and Photonics (Academic, 2007). J. M. Dudley and J. R. Taylor, eds., Supercontinuum Generation in Optical Fibers (Cambridge University, 2010). T. Brabec and F. Krausz, “Intense few-cycle laser fields: Frontiers of nonlinear optics,” Rev. Mod. Phys. 72, 545–591 (2000). A. L. Gaeta, “Nonlinear propagation and continuum generation in microstructured optical fibers,” Opt. Lett. 27, 924–926 (2002). N. Aközbek, A. Iwasaki, A. Becker, M. Scalora, S. L. Chin, and C. M. Bowden, “Third-harmonic generation and selfchanneling in air using high-power femtosecond laser pulses,” Phys. Rev. Lett. 89, 143901 (2002). G. Fibich, “Self-focusing in the damped nonlinear Schrödinger equation,” SIAM J. Appl. Math. 61, 1680–1705 (2001). M. Mlejnek, M. Kolesik, E. M. Wright, and J. V. Moloney, “Recurrent femtosecond pulse collapse in air due to plasma generation: numerical results,” Math. Comput. Simul. 56, 563–570 (2001). L. Bergé, S. Skupin, R. Nuter, J. Kasperian, and J.-P. Wolf, “Ultrashort filaments of light in weakly ionized, optically transparent media,” Rep. Prog. Phys. 70, 1633–1713 (2007). A. A. Zozulya, S. A. Diddams, A. G. Van Engen, and T. S. Clement, “Propagation dynamics of intense femtosecond pulses: multiple splittings coalescence and continuum generation,” Phys. Rev. Lett. 82, 1430–1433 (1999). G. Strang, “On the construction and comparison of difference schemes,” SIAM J. Numer. Anal. 5, 506–517 (1968). M. A. Celia and W. G. Gray, Numerical Methods for Differential Equations (Prentice-Hall, 1992). A. N. Krylov, “On the numerical solution of the equation by which in technical questions frequencies of small oscillations of material systems are determined,” Otdel. Mat. i Estest. Nauk 12, 491–539 (1931) [in Russian]. H. A. Van der Vorst, P. G. Ciarlet, A. Iserles, R. V. Kohn, and M. H. Wright, Iterative Krylov Methods for Large Linear Systems (Cambridge University, 2003). H. A. van der Vorst, “BI-CGSTAB: a fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems,” SIAM J. Sci. Stat. Comput. 13, 631–644 (1992). P. Bogacki and L. F. Shampine, “A 3(2) pair of Runge-Kutta formulas,” Appl. Math. Lett. 2, 321–325 (1989). L. F. Shampine and M. W. Reichelt, “The MATLAB ODE suite,” SIAM J. Sci. Comput. 18, 1–22 (1997). J.-C. Diels and W. Rudolph, Ultrashort Laser Pulse Phenomena: Fundamentals, Techniques, and Applications on a Femtosecond Time Scale, Optics and Photonics (Academic, 1996). J. B. Ashcom, “The role of focusing in the interaction of femtosecond laser pulses with transparent materials,” Ph.D. dissertation (Harvard, 2003). J. H. Marburger, “Self-focusing: theory,” Prog. Quantum Electron. 4, 35–110 (1975). E. M. Wright, W. J. Firth, and I. Galbraith, “Beam propagation in a medium with a diffusive Kerr-type nonlinearity,” J. Opt. Soc. Am. B 2, 383–386 (1985).
20 February 2015 / Vol. 54, No. 6 / APPLIED OPTICS
1435
A model for nonlinear optical propagation is cast into a split-step numerical framework via a variable stencil-size Crank–Nicolson finite-difference method for the linear step and a choice of two different nonlinear integration schemes for the nonlinear step. The model includes Kerr, Raman scattering, and ionization effects (as well as linear and nonlinear shock, diffraction, and dispersion). We demonstrate the practical importance of numerical effects when interpreting computational studies of high-intensity optical pulse propagation in physical materials. Examples demonstrate the significant error that can arise in discrete, limited precision implementations as one attempts to improve practical operator accuracy through increased operator support size and sampling frequency. We also demonstrate the effect of the method used to obtain the finite-difference operator coefficients defining the equations ultimately used in the discrete model. Smooth, plausible, but incorrect solutions may result from these numerical effects. This implies the necessity of a complete, precise description of all numerical methods when reporting results of computational physics investigations in order to ensure proper interpretation and reproducibility. © 2015 Optical Society of America OCIS codes: (190.5530) Pulse propagation and temporal solitons; (190.7110) Ultrafast nonlinear optics; (350.5500) Propagation. http://dx.doi.org/10.1364/AO.54.001426
1. Introduction
Mathematical modeling of intense optical pulse propagation physics via the nonlinear Schrödinger equation has applications in an increasing number of fields. Optical parametric amplifiers and oscillators [1,2], second-harmonic generation [3] and third harmonic generation [4], Raman scattering [5], optical bistability [6], solitons [7–9], and many other phenomena have found numerous practical applications. In support of these applications, as well as investigations into the nature of the physical phenomena themselves, the scientific and engineering communities expend significant effort in numerical modeling [10–13]. It is important to realize that a numerical 1559-128X/15/061426-10$15.00/0 © 2015 Optical Society of America 1426
APPLIED OPTICS / Vol. 54, No. 6 / 20 February 2015
model tends to introduce its own contributions to the intended physics model, altering the description of pulse evolution. Numerical physics investigations must therefore attempt to understand and separate the numerical from the physical. The numerical examples in this paper demonstrate the importance of completely describing not only the physical but also the numerical model when reporting on such investigations. The literature on nonlinear propagation is very broad [13]. A variety of models and unit systems are used to describe various subsets of known physical effects (see Brabec and Krausz [14], Couairon and Mysyrowicz [13], or Zozulya et al. [15]). We summarize the complete propagation and material interaction model used in our work in Section 2. This model is derived in detail in [16], bringing all physical phenomena of interest under the same framework and
clarifying the necessary physical assumptions and approximations. The primary differences between the work described in this paper and that in well-known canonical references such as Agrawal [17] and Dudley and Taylor [18] are as follows. The physical problems we are interested in (free-space propagation) require a diffraction model represented by the transverse Laplacian in Eq. (1) of the model summary in Section 2. This operator necessitates additional numerical methods discussed in Section 3.A. Also, in addition to the Kerr effect and Raman scattering we have models for ionization [see Eqs. (8–14) in Section 2]. Many frameworks without ionization use a Taylor series expansion for the Raman response in Eqs. (6) and (7) (see Agrawal [17]), which simplifies the nonlinear propagation step. This approach is related to the simple integration described in Eq. (42). When combining both Raman and ionization effects, this method is problematic. For this reason we use the more general approach in Eq. (43), which allows inclusion of very general material interaction models. In Section 3 we develop a numerical model for the physics. We describe the overarching numerical method (split-step), which breaks the problem into linear and nonlinear operators applied successively over short intervals along the propagation coordinate axis. Split-step has a rigorous underpinning despite an often heuristic description. It can be implemented with a reasonable amount of effort. Perhaps most importantly for physicists and engineers, it is straightforward to modify existing or add new physical submodels that would otherwise be unwieldy. For the linear operator, the dispersion and linear shock are evaluated via Fourier transform in the frequency domain. This method is fast and provides a very accurate implementation of the continuous temporal operators. However, it does require regular sample spacing along the temporal axis and that we keep the pulse away from the initial and final grid times to prevent a physical scattering from the periodic boundary conditions. The transverse Laplacian (diffraction) is implemented via finite differences using a Crank–Nicolson (CN) method. This requires the solution of a series of large linear systems (one per frequency in the temporal FFT mentioned above). CN is an implicit, unconditionally stable method that typically admits larger steps than explicit methods (which must satisfy the Courant–Friedrichs–Levy criterion). It is nondissipative and therefore does not introduce artificial loss into the system. Its principle drawback for our application is that it is dispersive and can sometimes suffer from oscillations. We derive all these operators in Section 3.A.1 and set up the computational framework in Section 3.A.2 to easily vary the size of the operator support. We also discuss the use of a preconditioner and its effect on the accuracy and speed of the solution for the large linear solves in Section 3.A.3.
Section 3.B describes methods for solving the nonlinear operator portion of the split-step framework. For our physical model there are actually two equations that must be integrated in this step: a description of the evolution of electron number density and the nonlinear portion of the split-step method. Both may be solved via Runge–Kutta methods. For problems in which the nonlinear portion of the splitstep framework takes a particularly simple form, we show that a simple nonlinear interest model suffices. In each section, we give numerical examples showing the effect of different choices for various numerical model parameters, such as operator support size, sample spacing, z-direction step size, and use of a preconditioner. The results are summarized in Section 4. 2. Summary of the Physical Model
Here we summarize the physical model that we discretize and solve in Section 3. It includes paraxial diffraction, group velocity and third-order dispersion, linear and nonlinear shock, Kerr and Raman effects, and ionization. There are many models for nonlinear optical propagation in the literature [13,14,19–25]. Ours is similar to those in the aforementioned references. However, our method of including the impedance effects in the plasma defocusing term leads to a difference in the shock terms (see, for example, [13]). Reference [16] derives the model described here in detail with a clear explanation of all assumptions and approximations. The propagation equation for the electric field envelope E is ∂ζ E
i −1 2 T ∇⊥ E iDE TNE; 2k
(1)
with wavenumber k 2π∕λ and wavelength λ. The dispersion and shock operators are given by k00 k000 D ≡ − ∂2τ − i ∂3τ 2 6
and
i T ≡ 1 ∂τ ; ωo
(2)
respectively. The nonlinear operator may be decomposed as N N kr T −1 N pls T −1 N MPA :
(3)
Kerr and Raman effects are modeled with N kr iko n2 1 − αI αQR :
(4)
We have used the approximately time-harmonic relationship I≈
cno εo 2 jEj 2
(5)
to specify the intensity envelope I. Delayed Raman response is described by 20 February 2015 / Vol. 54, No. 6 / APPLIED OPTICS
1427
Z QR r; t; z
t −∞
RτIr; t − τ; zdτ;
(6)
(13)
If we do not wish to include ionization saturation effects, we may instead use
with Rt
γ 2 ω2R exp−γt sinωR t; ωR
(7)
where ωR and γ are, respectively, characteristic frequency and attenuation constants for the molecular species interacting with the pulse. The following quantities help us express the remaining ionization physics: τc Electron collision or relaxation time ωo Pulse center frequency no Background refractive index U i Energy required to free single electron ρnt Neutral atom number density αr Electron-positive ion recombination loss σ k Ionization cross-section: For plasma defocusing we use (8)
or, if we desire to include collision effects, (9)
with ρc ≡ εo me ω2o ∕e2 ;
and σ ≡
ko ωo τc : no ρc 1 ω2o τ2c
(10)
The energy loss in the pulse due to multiphoton absorption is included via N MPA −
βK K−1 I ; 2
(11)
with βK σ k ρnt Kℏωo :
We seek to solve an equation of the form ∂z A LA NAA:
(15)
u0 uo ;
(16)
where L may be a partial or ordinary differential operator. Use of the Baker–Campbell–Hausdorff formula leads to the development of higher order stepping patterns, but we consider only the firstand second-order frameworks here. For these, we solve the two equations
K Photons required to free single electron
σ N pls ≡ − 1 iωo τc ρ; 2
(14)
3. Numerical Model
ut Lu;
ko Pulse center wavenumber
iko T −1 ρE; N pls E − 2no ρc
dρ σ β 2 ρI K I K − αr ρ2 : Kℏωo dt no U i
The split-step approach (see Strang [26]) can be derived rigorously by considering an initial value problem of the form
me Mass of an electron
(12)
The evolution of the electron density requires an auxiliary equation: 1428
dρ σ 2 ρI I K σ k ρnt − ρ − αr ρ2 : dt no U i
APPLIED OPTICS / Vol. 54, No. 6 / 20 February 2015
∂z A LA;
(17)
∂z A NAA
(18)
in an alternating fashion. Solving each of these equations will require additional numerical methods. The numerical evolution of the pulse may take place along either the temporal or z coordinate (where z is the principal propagation direction). We describe the evolution of the pulse through a plane over time at successively increasing z locations. The numerical grid has transverse spatial coordinates r or x; y and a temporal coordinate. At any given step in the propagation algorithm, this grid exists at a particular location in the primary propagation direction and describes the spatial and temporal variations of the pulse over the transverse plane. Each iteration of the propagation algorithm transforms the transverse space–time description at the current location to the distribution at the next location in the propagation direction. It is a z-propagated description of the pulse evolution. We begin with a description of the initial pulse and solve Eq. (17) over a short distance (z step) along the direction of propagation. In other words, we linearly propagate the pulse a distance dz. Section 3.A describes the linear numerical model used to accomplish this task. We use the resulting pulse description as an initial condition for Eq. (18) as described in Section 3.B. This completes a single first-order z-direction step. The process of alternating solution of the linear and nonlinear operator
equations continues until the desired number of zdirection steps have been taken. A second-order solution consists of taking half a linear step, an entire nonlinear step, and then another linear half-step. Diffraction and self-focusing can balance one another. However, in radially symmetric models where only a single filament is present there is a limit on the amount of power in the pulse that diffraction can balance. The same is true of free-electron defocusing in problems that involve ionization. Eventually, either fully (3D 1) models must be used so that multiple filaments are allowed or the pulse will collapse. A.
Linear Operator
Our basic framework for solving the linear equation could take many forms. The complex exponential is an eigenfunction of linear shift-invariant (LSI) operators. This fact leads to the prevalence of spectral methods in the solution of LSI partial differential equations in circumstances where the mathematically implied and physical boundary conditions are compatible. The implied periodicity of the discrete Fourier transform leads to scattering off the edges of the domain. Waves hitting the boundary have corresponding sources from an adjacent domain radiating energy into our numerical grid. In many situations (converging beams, etc.) this would be acceptable. However, even relatively low amplitude disturbances near the edge of the pulse may cause problematic numerical disturbances if they propagate to the transverse edges of the domain and reflect back to the center. For these reasons we use a Fourier-transformbased method in time only. We transform the linear operator into the temporal frequency domain leaving only the transverse derivatives. We then solve the resulting linear algebraic problem for each temporal frequency in the pulse. This allows extension of the method to include absorbing boundary conditions and adaptive mesh refinement (see Mlejnek et al. [10] for the first reported simulations of nonlinear propagation via AMR). 1. Basic Differential Operators We begin with an approximation to the nth derivative at location x taken over a support of q samples. We may approximate the derivative [27] q X dn u ≈ γ u γ u γ u γ i ui : 1 1 2 2 q q dxn x i1
(19)
To obtain the necessary coefficients in Eq. (19) for an odd-length stencil of size q 2p 1, with a derivative of order n, evaluated at the center stencil location we may solve PΓ l; where
(20)
q−1 h Phk k − ; 2
2
Γ1 Γ2 Γ3 .. .
3
7 6 7 6 7 6 7 6 7; 6 Γ6 7 7 6 6Γ 7 4 q−1 5 Γq
h; k 0; …; q − 1;
(21)
2
3 0 . 6 .. 7 6 7 607 6 7 6 7 and l 6 n! 7; 6 7 607 6 . 7 4 .. 5 0
(22)
and the only non-zero entry n! on the right-hand side is in the m th position (the initial position is designated zero). We then form γ Δ−n Γ;
(23)
where Δ is the physical spacing between samples. By using Eqs. (20)–(23) we obtain the finite-difference coefficients for a given derivative order and sufficiently broad spatial support. Figure 1 shows plots of the error obtained by using finite-difference operators with coefficients obtained in this fashion via three different linear algebraic solution methods. We use a simple unit period sinusoid as our test function, computing the derivative at x 1∕4 via finite differences. We then compute ϵ, the absolute value of the difference between the finitedifference estimate and the analytical derivative 2π cos2πt at the same location. Each plot shows log10 ε as a function of the operator support and log10 dt, where dt is the finite-difference sample spacing. Operator coefficients obtained via Eqs. (21) and (22) are consistent approximations, meaning that as the finite-difference spacing goes to zero the approximations must improve [27]. Initially, we see this trend in all three plots. As we move from the smallest possible odd-length support and the largest shown sample spacing toward longer filter supports and smaller sample spacing, the accuracy of the operator improves. However, in both directions we eventually see that numerical inaccuracies in the coefficients lead to a reduction in the accuracy of the operator. This means that while the finite-difference operator is technically consistent, practical considerations such as round-off error caused by finite precision arithmetic do not always allow us to obtain improved accuracy by decreasing the sample spacing or increasing the support of the filter. As we move between the plots we see that the accuracy of the numerical method used to obtain the coefficients for the finite-difference operator is important. Using a numerically stable method gives clearly better coefficients. This is particularly important in AMR frameworks where variable grid spacing requires us to obtain coefficients for many different sampling patterns. 20 February 2015 / Vol. 54, No. 6 / APPLIED OPTICS
1429
−2 −4 −6
−2
−10
−2 0 −6−4
−8
−3
−12
−8
−2 −4 −6
−8
−10
−4 −8 log (dt)
−8
−2 −4 −6
LU Decomposition with Partial Pivoting −8−10 −10
−6
−8
−12
−8−10 −12
10
−5 −10
−6
−4
−7
−4 −5
−10 −10
−10
−2 −3
−12
−10
−12 −12
−1
−6 −10
−6
log10(dt)
Gaussian Elimination
Moore−Penrose Pseudo−Inverse −1
−7
−8 −8
−8
−6
−9
−2
−4
−8
−2 −4
−8 −9
0 −6
−10
3
5
7 9 11 13 Operator Support (samples)
15
3
5
7 9 11 13 Operator Support (samples)
(a)
3
15
5
7 9 11 13 Operator Support (samples)
(b)
15
−10
(c)
Fig. 1. Finite-difference coefficient and operator accuracies. Each contour plot shows the base 10 logarithm of the absolute value of the difference between the finite-difference operator estimate and the analytical first derivative of the unit period sinusoid sin2πx at x 1∕4. (a) Finite-difference coefficients obtained by solving Eqs. (20)–(22) via a Moore–Penrose pseudo-inverse. (b) Coefficients obtained via Gaussian elimination. (c) Coefficients obtained via lower–upper (LU) decomposition with partial pivoting.
Significantly reducing the radial sample spacing and z-step size (magenta curve) alters the solution only marginally. If we push the operator support and number of radial samples too high the solution begins to degrade (green curve). Any further increases in either operator support or the number of radial samples causes the numerical solution to break down entirely (likely due to round-off error).
Figure 2 shows the effects of varying operator support, radial sample spacing, and the number of z steps used to compute a solution. Results for numerical and analytical linear propagation of a focused Gaussian beam are shown in Fig. 2(a). We observe that with three-point operator support and 64 radial samples that the quality of the numerical solution is low. We may increase either operator support, the number of radial samples, or both in order to obtain a better solution. In Fig. 2(b) we see similar effects for nonlinear propagation. The red curve has minimal operator support, radial sampling, and a small number of z steps. Increasing only the operator support changes the solution dramatically without requiring significantly more computer resources (blue curve).
2. Discretization of Propagation Equation Here we describe the discretization of the linear spatial operator including the application of relevant boundary conditions. We do so in a fashion that allows the programmatic specification of differential Stencil size: 3, Radial samples: 90, z−steps: 128 Stencil size: 13, Radial samples: 90, z−steps: 128 Stencil size: 13, Radial samples: 256, z−steps: 256 Stencil size: 17, Radial samples: 700, z−steps: 256
Intensity (normalized)
Amplitude (Normalized)
3
2
Stencil size: 3, Radial samples: 64, z−steps: 200 Stencil size: 3, Radial samples: 256, z−steps: 200 Stencil size: 13, Radial samples: 64, z−steps: 200 Stencil size: 13, Radial samples: 256, z−steps: 200 Analytical solution
1
0
0
1
2
3
3
2
1
4
0
−150
−100
−50
0
Distance (m)
Time (fs)
(a)
(b)
50
100
150
Fig. 2. Variable stencil size. (a) Linear propagation of a focused Gaussian beam with focal distance f 3 m, e−1 beam radius re 1.5 mm, and wavelength λ 810 nm. The dotted line shows the analytical solution. Each numerical solution has identical z-step spacing. It is clear from the plots that a significant increase in accuracy may be achieved either by decreasing the radial sample spacing or increasing the operator support. Increasing the operator support requires much fewer computer resources. (b) Nonlinear propagation with linear and nonlinear shock using a linear algebraic preconditioner for solution of finite-difference equations. For this model the linear refractive index of the material is no 1.45, nonlinear refractive index n2 2.5 × 10−20 m2 ∕W, the group velocity dispersion is 3.6 × 10−26 s2 ∕m, the thirdorder dispersion coefficient is zero, the coefficient controlling the percentage of instantaneous Kerr to Raman effect is γ 0.15, and the characteristic Raman response time and frequency are, respectively, τR 50 × 10−15 s and ωR 84 × 1012 Hz. The beam is propagated 3 cm. The curves show the effect of moving from minimal operator support, number of radial samples, and number of z steps through increases in each. The blue and magenta curves are computed with the most accurate operators. Further increases in radial samples or operator support cause the solution to break down very quickly. 1430
APPLIED OPTICS / Vol. 54, No. 6 / 20 February 2015
operator order and, hence, stencil size. This provides a means of developing an empirical understanding of how operator order affects solution accuracy when combined with sample spacing in the presence of round-off error. In the frequency domain the linear portion of Eq. (1) takes the form ∂z u iq1 u iq2 ∇2⊥ u;
(24)
with q1 and q2 real. Define the operator dm n to be an nth order derivative in space (in our case the rˆ direction) with support of size m. For example, d31
1 −1; 0; 1: 2Δ
(26)
In cylindrical coordinates the Laplacian is 1 ∇2⊥ ≡ ∂2r ∂r : r
(27)
The CN approximations are then 1 u → un1;j un;j ; 2
∂2r u →
(28)
1 u − un;j ; δz n1;j
(29)
1 m d · un1;jm un;jm ; 2 2
(30)
∂z u →
1 1 ∂ u → dm · un1;jm un;jm ; r r 2r 1
(31)
where the notation 2
un;jm
3 un;j−m−1∕2 6 un;j−m−1∕21 7 6 7 .. 6 7 7 6 . 7 6 7 u 6 n;j 7 6 7 6 . .. 7 6 7 6 4u 5 n;jm−1∕2−1
un;j−m−1∕2
m pm j · un1;jm pj · un;jm ;
(33)
where pm j ≡
2 m 1 m m d d : d0 − i q1 dm q 2 0 2 δz r 1
(34)
Note that at j 0 we must use L’Hospital’s rule to arrive at
(25)
The operator in Eq. (25) includes the grid spacing Δ. In general, the operator is inversely proportional to the n th power of the grid spacing for an order n derivative. Section 3.A.1 describes how we arrive at the coefficients. Without loss of generality we concern ourselves only with odd support lengths. By definition, the dm 0 is a support centered Kronecker delta function: d50 0; 0; 1; 0; 0:
describes an m-length vector evaluated at step n in the z-coordinate direction, centered at grid node index j in the r-coordinate direction. Substituting into Eq. (24) gives (after some algebra)
pm 0 ≡
2 m m d − iq1 dm 0 2q2 d2 : δz 0
Equation (33) is an algebraic equation in m unknowns. We know the un; jm (they are the field values at the particular frequency specified by the values of q1 and q2 as a function of radial coordinate at z-step index n) and the entries of the row vectors pm j . The right-hand side of Eq. (33) is a scalar. The m unknowns are un1; jm (the field values at the particular frequency specified by the values of q1 and q2 as a function of radial coordinate at z-step index n). If we add another equation centered at an adjacent value of j (the radial coordinate index), then we have an additional equation and an additional unknown (two equations in six unknowns). Each time we add an adjacent equation we get an additional variable. If we add an adjacent equation until we reach a boundary, the new variable added by each boundary equation is known. At each boundary we add an equation with a known additional variable m − 1∕2 times. After reaching both boundaries we have exactly the number of equations we need to solve the system. Solving the system lets us update all the un;j to un1;j for all j at once. This method (CN) is implicit (uses both the current and future time to obtain the future time) and, hence, requires a solve (in this case a linear solve) at each step. It is unconditionally stable and will tolerate a larger step size than explicit methods (which are stable only if the step size is chosen according to the Courant stability criterion). The price for this unconditional stability is a global linear algebraic solve at each step. Practical implementation for a grid with J radial samples requires that we use Eq. (34) and the known current values of un to create a matrix equation: Lun1 L un ;
32
(35)
(36)
where L is a matrix of dimensions J × J and un is a known vector of length J. Row j of L is formed from Eq. (34). The top and bottom m − 1∕2 rows are special as they touch the boundaries. If we apply symmetric boundary conditions at the j 0 (r 0 20 February 2015 / Vol. 54, No. 6 / APPLIED OPTICS
1431
boundary) and a Dirichlet boundary condition with a value of zero at the outer boundary, noting that −l −l1 …p0j …pjl−1 plj ; pm j pj pj
(37)
where l m − 1∕2, then we may form the augmented matrix 2 6 6 6
p−l 0
p0−l1 p−l 1
0 … p0 p1−l1 …
… p01
p0l−1 … .. .
pl0 p1l−1
pl1 ; (38a)
.. p−l J−2
−l1 pJ−2 p−l J−1
. …
−l1 pJ−1
p0J−2 …
p1J−2 p0J−1
7 7 7 7 5;
38b
where Eq. (38a) describes the upper-left portion of the matrix and Eq. (38b) shows the lower right. The vertical line in Eq. (38a) represents the central propagation axis of our numerical grid. We must apply a symmetry boundary condition here. Mathematically, this implies that we fold the entries to the left of this line over onto those on the right-hand side and add. Equation (38b) shows a Dirichlet boundary condition with a value of zero on the right, meaning we do not need to extend beyond the matrix edge. We simply truncate when the center of the pm j reaches the edge (this will, of course, always be at row J − 1 when we start with index zero). With this understanding we can easily vary the stencil size via a parameter specified at run time. 3. Linear Solvers and Preconditioners Application of the finite-difference method generates large systems of linear equations for each frequency. Solving these systems implements the forward linear propagation at each frequency. Often the structure of the matrix is exploited to obtain a stable solution. For very large matrices, parallel computation may provide a means of propagating pulses requiring very large numerical grids. In these cases it is difficult to exploit the matrix structure. Instead, Krylov methods are used [28,29]. These methods find the solution using matrix-vector products only. These can be formed efficiently on distributed systems. For these examples, we use the biconjugate gradient stabilized (BICGStab) method [30]. In doing so we have found that preconditioning is necessary to obtain an accurate solution. The condition number of a square matrix is the magnitude of the ratio of the largest to the smallest eigenvalue. It is a measure of the sensitivity of the solution to a linear system of equations to errors in the data (which may result from the limited precision of the arithmetic). The purpose of preconditioning is to obtain a system of equations whose condition number is near unity (ideal case). This means that all of the eigenvalues 1432
are clustered near one another, which improves the accuracy and convergence of Krylov space methods. Instead of solving
APPLIED OPTICS / Vol. 54, No. 6 / 20 February 2015
Ax b;
(39)
P−1 Ax P−1 b;
(40)
we solve
where P−1 A has a better condition number than A. In practice this means that the matrix P is in some sense close to a scalar multiple of the system matrix A. In our numerical simulations we use the system matrix as our preconditioner. It would seem that multiplication by the inverse of the preconditioner would solve the system, but due to the numerical inaccuracies there remains some error, which the iterative solver then reduces below a specified tolerance. Figure 3 shows the effect of using a preconditioner for the linear solves in the case of a beam with critical power propagating in a material with a pure Kerr nonlinearity. The mathematical instability of this scenario (no ionization) creates an opportunity for evaluating the accuracy of our numerical methods. Numerical errors in the simulation will have the same effect as very small refractive index perturbations in a physical experiment. The less accurate the simulation, the sooner the pulse will lose confinement. Each example propagates the pulse the same physical distance. Figure 3(a) shows that 1000 z steps provide sufficient numerical accuracy with or without a preconditioner. Both curves show the expected confinement. Figure 3(b) compares the results obtained when only 500 z steps are used. With a preconditioner (blue curve) the linear solver is sufficiently accurate to allow diffraction and Kerr self-focusing to balance each other effectively. Without a preconditioner (green curve), inaccuracies in the linear propagator accrue, causing the beam to immediately diverge after initially self-focusing. Figure 4(a) shows the effect of the linear preconditioner with both Kerr and Raman effects in the presence of linear and nonlinear shock. The plot shows the normalized axial envelope intensity after propagating 3 cm in fused silica (compare with Zozulya et al. [15]). It is clear that preconditioning the linear solves has a significant effect on the solution, leading to more accurate results with fewer z steps and a lower number of radial samples. The BICGStab method is also significantly faster when using a preconditioner. B. Nonlinear Operator
In this section we describe the implementation of the nonlinear material response operator. This includes instantaneous Kerr effects and Raman scattering. These two effects have different physical origins but are often combined in numerical models due to the fact that they are both proportional to the cube of the electric field value. The nonlinear material model also includes ionization described by Eq. (13). If we
1
Nz=1000, Nr=512, with PC Nz=1000, Nr=512, No PC Normalized pulse width
Normalized pulse width
1
0.9
0.8
0.7
0.6
0
50
100
150
200
250
300
Nz=500, Nr=1024, with PC Nz=500, Nr=1024, No PC
0.9
0.8
0.7
0.6
0
50
100
150
200
Propagation distance (m)
Propagation distance (m)
(a)
(b)
250
300
Fig. 3. This simulation shows the effect of using a linear preconditioner when the propagation takes place in a simple Kerr material (no Raman scattering or ionization). Both linear and nonlinear shock are turned off. The pulse-center wavelength is λ 1024 nm. The linear and nonlinear refractive indices are n 1 and n2 4 × 10−23 m2 ∕W, respectively. Group velocity dispersion is 2.6 × 10−29 s2 ∕m, and thirdorder dispersion is given as 1 × 10−85 s3 ∕m. The initial beam has e−1 beam radius re 6 mm, an initially flat wavefront (focus at infinity), power at 95% of the critical power, and a propagation distance of 300 m. All simulations used finite-difference operators with a support of 3 and a single temporal sample (beam). (a) The blue curve shows the result of a simulation with preconditioner, 512 radial samples over 5re, and 1000 z steps. The green is the result of an identical simulation with no preconditioner. (b) The blue curve shows the result of a simulation with preconditioner, 1024 radial samples over 5re, and 500 z steps [half as many steps per unit distance as in (a)]. The green curve is the result of an identical simulation with no preconditioner.
do not know the neutral atom density we may use instead Eq. (14). We also include a loss term in the nonlinear portion of the propagator (discussed in [16]) that is designed to account for energy loss due to multiphoton ionization by removing the requisite amount of energy from the pulse every step. The nonlinear portion of our split-step operator method is an equation of the form ∂ A NAA: ∂z
(41)
When the nonlinear operator in Eq. (41) can be written in the form where NA contains no differential operators applying to the final envelope expression A [in this case NA reduces to a number], then we have an easier problem. This does not mean there are no differential operators in NA. Sometimes there are shock terms that would normally operate on the final envelope A except for a judicious manipulation of the operator expression. When this simpler circumstance prevails, Eq. (41) is essentially a nonlinear first-order interest equation. For initial condition Az propagated over a distance h, the solution to this problem may be written Z Az h; r exp
zh z
NAdz0 Az
≈ exphNAAz:
(42)
When we wish to include the shock operators in a straightforward fashion (or cannot manipulate the expression to exclude application of differential operators to the final envelope), we have a differential equation of the form
d Az; t f Az; t; z; dz
z1 ≤ z ≤ z2 z1 h; (43)
with initial condition Az1 ; t given. We may integrate the equation with an explicit one-step Runge–Kutta(2, 3) method from Bogacki and Shampine [31], which is useful for moderately stiff ordinary differential equations of the form in Eq. (43). We must integrate the nonlinear ionization model in Eq. (13) [or Eq. (14) if the neutral atom density is not known] over the temporal duration of the pulse each time we perform the nonlinear step. We may use the same RK-(2, 3) method from and Shampine as above. The computations shown in this paper use this method. However, the solution of these equations can go through sharp transitions and regions of much slower evolution depending on the structure of the optical pulse. For this reason it is sometimes necessary to use ODE solvers designed for stiff differential equations (for example, those described in Shampine and Reichelt [32]). Notice that the equations for the nonlinear portion of the split-step method and the electron number density are coupled and should technically be solved simultaneously. We make the approximation that they may be solved sequentially, solving first for the evolution of number density using the output of the previous linear step and then using the result in the integration of the nonlinear splitstep equation. Critical power. There are many expressions for critical power depending on the model for transverse beam or pulse profile, the type of nonlinear material, and the criteria for entrapment (see, for example, Diels and Rudolph [33], Ashcom [34], Zozulya et al. [15], Marburger [35], or Wright et al. [36]). For 20 February 2015 / Vol. 54, No. 6 / APPLIED OPTICS
1433
Axial Envelope at 3 cm (λ = 800 nm)
3.5
Nz: 101, Nr: 512, No PC Nz: 111, Nr: 256, No PC Nz: 101, Nr: 512, PC Nz: 101, Nr: 256, PC Nz: 501, Nr: 512, PC
2.5
Linear, Nt: 64 Nonlinear, Nt: 32 Nonlinear, Nt: 64 Nonlinear, Nt: 512
8 Amplitude (Normalized)
Normalized axial intensity
3
Normalized Maximum Axial Field Amplitude
10
2 1.5 1
6
4
2 0.5 0 −200
−150
−100
−50
0 Time (fs)
50
100
200
150
(a)
0
0
0.01
0.02 0.03 Distance (m)
0.04
0.05
(b)
Fig. 4. (a) Simulation pulse-center wavelength λ 800 nm. Material parameters no 1.45, n2 2.5 × 10−20 m2 ∕W, group velocity dispersion 3.6 × 10−26 s2 ∕m, and no third-order dispersion. The coefficient controlling the percentage of instantaneous Kerr to Raman effect is γ 0.15, the characteristic Raman response time and frequency are, respectively, τR 50 × 10−15 s and ωR 84 × 1012 Hz. Both linear and nonlinear shock are turned on. The initial beam has e−1 beam radius re 59.45 μm, an initially flat wavefront (focus at infinity), a temporal FWHM of 90 fs, and power Pin 4.72 mW (approximately 1.83Pcr;Z ). The propagation distance is 3 cm. The spatiotemporal extent of the numerical grid is 4 times the temporal pulse duration and 4 times the radial extent. All simulations have 128 temporal samples. As indicated in the plot legend, the number of z steps, the number of radial samples, and the use of a preconditioner vary among the simulations. The two extreme examples (red and dashed red curves) do not use a preconditioner. All three remaining curves use a preconditioner during the linear, frequency-domain finite-difference update for each frequency and produce results very close to one another despite differences in step number and sample spacing. (b) This simulation shows the effects of varying the temporal sample rate. The pulse propagates 5 cm with no linear or nonlinear shock using 300 z steps. The linear solver uses the linear system matrix as its own preconditioner. The central wavelength is λ 586 nm, and the linear and nonlinear refractive indices are, respectively, no 1.00027 and n2 4.9 × 10−23 m2 ∕W. The group velocity dispersion is 2.6 × 10−26 s2 ∕m and there is no third-order dispersion. The initial beam has e−1 beam radius re 200 μm and temporal extent 200 fs. The grid has temporal and spatial extents of 1800 fs and 1800 μm, respectively. Each simulation uses 256 samples along the radial coordinate. The geometric focus is 2.5 cm and the input power is Pcr .
the Gaussian beam in the simulations reported in Figs. 3 and 4(b) we use Pcr
λ2 : 2πno n2
(44)
For the simulations in Fig. 4(a) we compute per Zozulya et al. [15]: Pcr;Z
π0.61λ2 : 8no n2
(45)
4. Summary and Conclusions
In this paper we developed mathematical and numerical physics models for intense optical pulse propagation in homogeneous, isotropic, nonlinear media. We included linear dispersion, diffraction, nonlinear Kerr and Raman effects, and ionization, as well as linear and nonlinear shock terms. The mathematical physics propagator is a one-way paraxial scalar model. The numerical split-step algorithm was explicitly described. We showed how the differential operators of varying support are derived for the linear step. Implementation of the temporal derivatives (including shock terms) takes place in the frequency domain via Fourier transform. The transverse Laplacian is modeled via a CN finitedifference method using a BICGStab solver with a linear preconditioner. Ionization models describing 1434
APPLIED OPTICS / Vol. 54, No. 6 / 20 February 2015
the evolution of electron number density are integrated via a well-known explicit, one-step Runge– Kutta solver. The nonlinear split-step integration uses either a simple, nonlinear interest model or the same Runge–Kutta method used for the ionization. We showed that understanding the effects of the numerical framework is very important when using computation to gain physical insight. As seen in Fig. 1, even in mathematically consistent finite-difference approximations, practical operator accuracy does not necessarily increase with operator support size and sampling frequency due to significant truncation and round-off error effects. Also, when programmatically solving for the finite-difference coefficients themselves, the solution method can have a marked effect on their accuracy and hence the subsequent trade-off between truncation and round-off error. In some cases, significant improvements in accuracy or reductions in run time may be obtained with larger stencils instead of decreasing sample spacing [see Figs. 2(a) and 2(b)]. Clearly demonstrated in these figures is the ease with which we might obtain a smooth, plausible, but incorrect solution. This also implies that some care must be taken when reporting results to precisely and completely describe the numerical methods to ensure meaningful interpretation and reproducibility. Figure 3 shows the effect of using a preconditioner compared with decreasing the step size in the direction of propagation. We may obtain accurate results with fewer steps. Also advantageous is the fact that the preconditioner generally
appears to significantly accelerate the linear solution in most cases. This manuscript has been authored by Sandia Corporation under Contract No. DE-AC04-94AL85000 with the U.S. Department of Energy.
17. 18. 19.
References 1. J. Hansryd, P. A. Andrekson, M. Westlund, J. Li, and P. Hedekvist, “Fiber-based optical parametric amplifiers and their applications,” IEEE J. Sel. Top. Quantum Electron. 8, 506–520 (2002). 2. G. Cerullo and S. De Silvestri, “Ultrafast optical parametric amplifiers,” Rev. Sci. Instrum. 74, 1–18 (2003). 3. P. A. Franken, A. E. Hill, C. W. Peters, and G. Weinreich, “Generation of optical harmonics,” Phys. Rev. Lett. 7, 118–119 (1961). 4. R. W. Terhune, P. D. Maker, and C. M. Savage, “Optical harmonic generation in calcite,” Phys. Rev. Lett. 8, 404–406 (1962). First observation of third-order harmonic generation. 5. E. J. Woodbury and W. K. Ng, “Ruby laser operation in the near IR,” Proc. IRE 50, 2365–2383 (1962). First observation of optical Raman scattering. 6. H. M. Gibbs, S. L. McCall, and T. N. C. Venkatesan, “Differential gain and bistability using a sodium-filled Fabry–Perot interferometer,” Phys. Rev. Lett. 36, 1135–1138 (1976). First demonstration and explanation of optical bistability. 7. A. Hasegawa and F. Tappert, “Transmission of stationary nonlinear optical pulses in dispersive dielectric fibers. I. Anomalous dispersion,” Appl. Phys. Lett. 23, 142–144 (1973). Theoretical prediction of soliton. 8. A. Hasegawa and F. Tappert, “Transmission of stationary nonlinear optical pulses in dispersive dielectric fibers. II. Normal dispersion,” Appl. Phys. Lett. 23, 171–172 (1973). Theoretical prediction of soliton. 9. L. F. Mollenauer, “The non-linear Schrödinger equation and ordinary solitons,” Unpublished presentation (Lucent Technologies, 2002). 10. M. Mlejnek, M. Kolesik, J. V. Moloney, and E. M. Wright, “Optically turbulent femtosecond light guide in air,” Phys. Rev. Lett. 83, 2938–2941 (1999). 11. A. Couairon, E. Brambilla, T. Corti, D. Majus, O. de J. Rámirez-Góngora, and M. Kolesik, “Practitioners guide to laser pulse propagation models and simulation,” Eur. Phys. J. Spec. Top. 199, 5–76 (2011). 12. L. Berge and A. Couairon, “Modeling the filamentation of ultra-short pulses in ionizing media,” Phys. Plasmas 7, 193–209 (2000). 13. A. Couairon and A. Mysyrowicz, “Femtosecond filamentation in transparent media,” Phys. Rep. 441, 47–189 (2007). 14. T. Brabec and F. Krausz, “Nonlinear optical pulse propagation in the single-cycle regime,” Phys. Rev. Lett. 78, 3283–3285 (1997). 15. A. A. Zozulya, S. A. Diddams, and T. S. Clement, “Investigations of nonlinear femtosecond pulse propagation with the inclusion of Raman, shock, and third-order phase effects,” Phys. Rev. A 58, 3303–3310 (1998). 16. T. A. Pitts, M. R. Laine, J. Schwarz, P. K. Rambo, B. M. Hautzenroeder, and D. B. Karelitz, “Derivation of an applied
20. 21.
22. 23.
24.
25.
26. 27. 28.
29. 30.
31. 32. 33. 34. 35. 36.
nonlinear schrödinger equation,” Tech. Rep. (Sandia National Laboratories, 2014). G. P. Agrawal, Nonlinear Fiber Optics, 4th ed., Optics and Photonics (Academic, 2007). J. M. Dudley and J. R. Taylor, eds., Supercontinuum Generation in Optical Fibers (Cambridge University, 2010). T. Brabec and F. Krausz, “Intense few-cycle laser fields: Frontiers of nonlinear optics,” Rev. Mod. Phys. 72, 545–591 (2000). A. L. Gaeta, “Nonlinear propagation and continuum generation in microstructured optical fibers,” Opt. Lett. 27, 924–926 (2002). N. Aközbek, A. Iwasaki, A. Becker, M. Scalora, S. L. Chin, and C. M. Bowden, “Third-harmonic generation and selfchanneling in air using high-power femtosecond laser pulses,” Phys. Rev. Lett. 89, 143901 (2002). G. Fibich, “Self-focusing in the damped nonlinear Schrödinger equation,” SIAM J. Appl. Math. 61, 1680–1705 (2001). M. Mlejnek, M. Kolesik, E. M. Wright, and J. V. Moloney, “Recurrent femtosecond pulse collapse in air due to plasma generation: numerical results,” Math. Comput. Simul. 56, 563–570 (2001). L. Bergé, S. Skupin, R. Nuter, J. Kasperian, and J.-P. Wolf, “Ultrashort filaments of light in weakly ionized, optically transparent media,” Rep. Prog. Phys. 70, 1633–1713 (2007). A. A. Zozulya, S. A. Diddams, A. G. Van Engen, and T. S. Clement, “Propagation dynamics of intense femtosecond pulses: multiple splittings coalescence and continuum generation,” Phys. Rev. Lett. 82, 1430–1433 (1999). G. Strang, “On the construction and comparison of difference schemes,” SIAM J. Numer. Anal. 5, 506–517 (1968). M. A. Celia and W. G. Gray, Numerical Methods for Differential Equations (Prentice-Hall, 1992). A. N. Krylov, “On the numerical solution of the equation by which in technical questions frequencies of small oscillations of material systems are determined,” Otdel. Mat. i Estest. Nauk 12, 491–539 (1931) [in Russian]. H. A. Van der Vorst, P. G. Ciarlet, A. Iserles, R. V. Kohn, and M. H. Wright, Iterative Krylov Methods for Large Linear Systems (Cambridge University, 2003). H. A. van der Vorst, “BI-CGSTAB: a fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems,” SIAM J. Sci. Stat. Comput. 13, 631–644 (1992). P. Bogacki and L. F. Shampine, “A 3(2) pair of Runge-Kutta formulas,” Appl. Math. Lett. 2, 321–325 (1989). L. F. Shampine and M. W. Reichelt, “The MATLAB ODE suite,” SIAM J. Sci. Comput. 18, 1–22 (1997). J.-C. Diels and W. Rudolph, Ultrashort Laser Pulse Phenomena: Fundamentals, Techniques, and Applications on a Femtosecond Time Scale, Optics and Photonics (Academic, 1996). J. B. Ashcom, “The role of focusing in the interaction of femtosecond laser pulses with transparent materials,” Ph.D. dissertation (Harvard, 2003). J. H. Marburger, “Self-focusing: theory,” Prog. Quantum Electron. 4, 35–110 (1975). E. M. Wright, W. J. Firth, and I. Galbraith, “Beam propagation in a medium with a diffusive Kerr-type nonlinearity,” J. Opt. Soc. Am. B 2, 383–386 (1985).
20 February 2015 / Vol. 54, No. 6 / APPLIED OPTICS
1435
