General vector auxiliary differential equation finite-difference time-domain method for nonlinear optics

The auxiliary differential equation finite-difference tim edomain method for modeling electromagnetic wave propagati on in dispersive nonlinear materials is applied to problems wher e the electric field is not constrained to a single vector component. A full-vect or Maxwell’s equations solution incorporating multiple-pole linear Lo rentz, nonlinear Kerr, and nonlinear Raman polarizations is presented. The a pplic tion is illustrated by modeling a spatial soliton having two orthog onal electric field components. To the best of our knowledge, the numerical technique presented here is the first to model electromagnetic wave pro pagation with two or three orthogonal vector components in dispersiv nonlinear materials. This technique offers the possibility of modeli ng sub-wavelength interactions of vector spatial solitons. © 2006 Optical Society of America OCIS codes:(190.5530) Pulse propagation and solitons; (190.3270) Ker r effect; (190.4420) Nonlinear optics, transverse effects in; (190.5650) Raman e ffect References and links 1. T. Kashiwa and I. Fukai, “A treatment by FDTD method of disper ive characteristics associated with electronic polarization,” Microwave Optics Tech. Lett. 3, 203-205 (1990). 2. R. M. Joseph, S. C. Hagness, and A. Taflove, “Direct time inte gration of Maxwell’s equations in linear dispersive media with absorption for scattering and propagation of femto sec nd electromagnetic pulses,” Opt. Lett. 16, 1412-1414 (1991). 3. S. Nakamura, Y. Koyamada, N. Yoshida, N. Karasawa, H. Sone, M . Ohtani, Y. Mizuta, R. Morita, H. Shigekawa, and M. Yamashita, “Finite-difference time-domain calculatio n with all parameters of Sellmeier’s fitting equation for 12 fs laser pulse propagation in a silica fiber,” IEEE Phot onics Technol. Lett. 14, 480-482 (2002). 4. S. Nakamura, N. Takasawa, and Y. Koyamada, “Comparison betwe en finite-difference time-domain calculation with all parameters of Sellmeier’s fitting equation and experi m ntal results for slightly chirped 12 fs laser pulse propagation in a silica fiber,” IEEE J. of Lightwave Technol. 23, 855-863 (2005). 5. M. Fujii, M. Tahara, I. Sakagami, W. Freude, and P. Russer, “ High-order FDTD and auxiliary differential equation formulation of optical pulse propagation in 2D Kerr and Raman nonlinear dispersive media,” IEEE J. Quantum Electron. 40, 175-182 (2004). 6. G. P. Agrawal,Nonlinear Fiber Optics, 3rd ed. Academic Press, San Diego, CA, 2001. 7. R. W. Boyd,Nonlinear Optics, 2nd ed. Academic Press, San Diego, CA, 2003. 8. A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Ti me-Domain Method, 3rd ed.Artech House, Norwood, MA, 2005. 9. A T. Ryan and G. P. Agrawal, “Spatiotemporal coupling in dis persive nonlinear planar waveguides,” J. Opt. Soc. Am. B12, 2382-2389 (1995). #72299 $15.00 USD Received 26 June 2006; revised 22 August 2006; accepted 24 August 2006 (C) 2006 OSA 4 September 2006 / Vol. 14, No. 18 / OPTICS EXPRESS 8305


Introduction
Emerging applications in nanophotonics involve electromagnetic wave interactions with materials having frequency-dependent and intensity-dependent polarizations.Auxiliary differential equation (ADE) methods extend the finite-difference time-domain (FDTD) method to incorporate polarization by time-stepping auxiliary differential equations with Maxwell's curl equations [1,2].Nakamura et al. validated experimentally an ADE FDTD model of nonlinear femtosecond ultrabroadband-pulse propagation in silica fiber [3,4].The Fujii ADE FDTD method [5] is an efficient reformulation which eliminates the need to solve a system of N equations at each time-step, where N is the number of poles of the chromatic dispersion.
In this paper, we provide details of how to apply Fujii's method to nonlinear optics problems where the electric field has two or three orthogonal vector components.We designate this technique the general vector auxiliary differential equation (GVADE) FDTD method.This method reproduces published results for temporal soliton propagation in a dispersive nonlinear material [5].The GVADE FDTD method is then used to model spatial soliton propagation with two orthogonal electric field components in fused silica.We use the same chromatic dispersion, nonlinear Kerr polarization, and nonlinear Raman polarization used by Nakamura et al., which models realistic physical effects in silica.

Electromagnetic wave propagation in dispersive nonlinear materials
To derive the GVADE method, we formulate Maxwell's equations using polarization current: where E and H are the electric and magnetic field vectors and J is the polarization current, ∂ P ∂t .In materials with multiple-pole linear Lorentz polarization, instantaneous Kerr nonlinear polarization, and Raman nonlinear polarization, the polarization current is J = J Lorentz + J Kerr + J Raman .The linear Lorentz polarization models the chromatic dispersion and contains a contribution from multiple resonances: where each J Lorentz p corresponds to the polarization current due to a single pole of the Sellmeier expansion where the phasor polarization is given by where β p and ω p are the strength and frequency, respectively, of the pth resonance [6].
In general, the third-order nonlinear polarization is where χ (3) is the third-order susceptibility tensor [7].For a simple model of the electron re- sponse accounting for nonresonant incoherent (intensity-dependent) nonlinear effects, the thirdorder nonlinear polarization can be described by the Born-Oppenheimer approximation, where g(t) is the causal response function and the induced polarization is assumed to lie in the same direction as the electric field [6].To model the Kerr and the Raman polarizations, where α represents the relative strengths of the Kerr and Raman polarizations, δ (t) is a Dirac delta function that models the instantaneous Kerr nonresonant virtual transitions, and where g Raman (t) is an approximation of the Raman response function with parameters τ 1 and τ 2 chosen to fit the Raman-gain spectrum [6], and U(t) is the Heaviside step function.From Eq. ( 6) and Eq. ( 7), the polarization and polarization current from the Kerr nonlinearity are Writing Eq. ( 6) as a convolution, the polarization from the Raman effect is where χ (3)

General vector auxiliary differential equation FDTD method
To implement the Fujii ADE method for the general electric field vector case, we use the formulation of Maxwell's equations which emphasizes the polarization current, J, instead of the polarization, P, because this eliminates the need to store and update the electric displacement field.The GVADE method can be used for any electric field having two or three orthogonal vector components.However, for a concise derivation, we consider the case where derivatives with respect to the z-coordinate are zero and the electromagnetic field is composed of E x , E y , and H z .Eq. ( 2) will be solved using FDTD combined with a semi-implicit method in which E n+1 will be updated in terms of H n+1/2 and J n+1/2 , where the superscript indicates the time-step.

Linear Lorentz polarization
Expression of Eq. ( 4) in the form of an equivalent polarization current phasor gives: Multiplying both sides of Eq. ( 13) by (ω 2 p − ω 2 ), and transforming to the time domain, yields Finite-differencing Eq. ( 14) centered at time-step n and solving for J n+1 Lorentz p yields where Using Eq. ( 15), and averaging across step n and n + 1, gives

Nonlinear Kerr polarization
The nonlinear Kerr polarization at step n + 1 2 is obtained from Eq. ( 10).The finite-difference expression for the time derivative of the Kerr polarization centered at step n + 1 2 is

Solution for electric field
Now that all terms in Eq. ( 2) are known at step n + 1 2 , E n+1 is determined from Because E is a vector quantity, solving Eq. ( 26) for E n+1 requires solving a nonlinear system of coupled equations.This is efficiently solved using a multi-dimensional Newton's method, where zeroes of objective functions are determined by iteration.For the field components {E x , E y , H z } considered here, Eq. ( 26) is solved for E n+1 x and E n+1 y by substituting Eq. ( 17), Eq. ( 18), and Eq. ( 25).We then define the objective function vector Next, we define G g x and G g y to be the gth guesses for E n+1 x and E n+1 y .Newton's method updates the guesses until the objective functions X and Y are sufficiently close to zero.Each subsequent guess is made from the current guess by: where J is the Jacobian matrix, ∂ (X,Y )/∂ (G x , G y ), with elements and | g indicates evaluation using the values from the gth guess.

GVADE FDTD simulation of temporal and spatial solitons
We applied the GVADE method to model temporal and spatial solitons in dispersive nonlinear materials.First, we modeled temporal soliton propagation in a material with single pole linear Lorentz dispersion and Kerr and Raman nonlinear polarizations.The same pulse, grid, and material from [5] were used; the results, shown in Fig. 1, reproduce those of [5]. Figure 2 illustrates the full-vector capability of the GVADE method.Here, we model the +xdirected propagation of a higher-order spatial soliton having the electromagnetic field vector components {E x , E y , H z }.We assume for the vector electric field the three-pole set of linear Sellmeier dispersions and the instantaneous Kerr and dispersive Raman nonlinearities published in [3,4].The strengths and resonant frequencies of the linear dispersions from Sellmeier's equation are given by: β 1 = 0.69617, β 2 = 0.40794, β 3 = 0.89748, ω 1 = 2.7537 × 10 16 rad/s, Fig. 1.GVADE simulation results for temporal soliton propagation in a dispersive nonlinear material.These results reproduce those of [5].
ω 2 = 1.6205 × 10 16 rad/s, and ω 3 = 1.90342 × 10 14 rad/s.The third-order electric susceptibility scalar is χ (3) 0 = 1.89 × 10 −22 m 2 /V 2 .The relative strengths of the Kerr and Raman polarizations are given by the parameter α = 0.7.The Raman polarization parameters are τ 1 = 12.2 fs and τ 2 = 32 fs.The soliton was modeled on a modified Yee grid with collocated electric field components [8].The grid was 6000 by 2615 cells with a spatial resolution of ∆x = ∆y = 8nm and a temporal resolution of ∆t = 3.34 × 10 −18 s.The spatial and temporal resolutions were chosen empirically based on convergence tests.The chosen time-step was one quarter of the time-step required by the Courant limit for standard FDTD [8].The magnetic field was excited by a hard source at the far-left side of the grid (x = 0) having an initial transverse profile where H 0 = 4.77 × 10 7 A/m, ω c = 4.35 × 10 15 rad/s (λ 0 = 433nm), and w = 667 nm.From Fig. 2, we see that the calculated solution exhibits a periodic expansion and contraction that is characteristic of higher-order solitons [9].Fig. 2. GVADE simulation results of a +x-directed higher-order spatial soliton with field components {E x , E y , H z } in a material with a three-pole Sellmeier linear dispersion, an instantaneous Kerr nonlinearity, and a dispersive Raman nonlinearity: magnitude of H z .

Conclusion
The importance of modeling electromagnetic wave propagation in dispersive nonlinear materials will increase with modern engineering of nanophotonic devices.In this paper, we showed how the Fujii ADE FDTD method can be applied to nonlinear optics problems where the electric field has two or three orthogonal vector components.The GVADE method is, to the best of our knowledge, the first numerical technique to model electromagnetic wave propagation with two or three orthogonal electric field vector components in dispersive nonlinear materials.We are currently exploring the relationships of predictions made by nonlinear Schrödinger (NLS) equation theory and the GVADE technique for a variety of scalar and vector spatial solitons.Our initial aim is to determine the transition between the the two approaches, i.e., where a computation based upon the full-vector Maxwell's equations is required to properly model the vector wave physics.Furthermore, we intend to explore a regime where the NLS equation is arguably inapplicable -the scattering of a spatial soliton by a sub-wavelength particle.