1. 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.

2. 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, $\frac{\partial \mathbf{P}}{\partial 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 _Lorentzp 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 χ ⁽³⁾ is the third-order susceptibility tensor [7]. For a simple model of the electron response accounting for nonresonant incoherent (intensity-dependent) nonlinear effects, the third-order 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 τ ₁ and τ ₂ 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 ⁿ⁺¹ will be updated in terms of H ^n+1/2 and J ^n+1/2, where the superscript indicates the time-step.

3.1. Linear Lorentz polarization

Expression of Eq. (4) in the form of an equivalent polarization current phasor gives:

Multiplying both sides of Eq. (13) by (${\omega }_p^2$ -ω ²), and transforming to the time domain, yields

Finite-differencing Eq. (14) centered at time-step n and solving for J ⁿ⁺¹ _Lorentzp yields

where

Using Eq. (15), and averaging across step n and n+1, gives

3.2. Nonlinear Kerr polarization

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

3.3. Nonlinear Raman polarization

Eq. (11) is solved by introducing a scalar auxiliary variable for the convolution,

with Fourier transform

where

Inserting Eq. (21) into Eq. (20), multiplying by (${\omega }_{\text{Raman}}^2$ + 2jω δ _Raman -ω ²), and transforming to the time domain yields the auxiliary differential equation,

Because Eq. (23) contains a second derivative, it is finite-differenced, centered at step n:

Using Eq. (11) and Eq. (19), the finite-difference expression for the time derivative of the Raman polarization centered at step $n+\frac{1}{2}$ is

3.4. Solution for electric field

Now that all terms in Eq. (2) are known at step $n+\frac{1}{2},{\mathbf{E}}^{n+1}$, E ⁿ⁺¹ is determined from

Because E is a vector quantity, solving Eq. (26) for E ⁿ⁺¹ 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_x^{\mathrm{n}+1}$ and $E_y^{n+1}$ by substituting Eq. (17), Eq. (18), and Eq. (25). We then define the objective function vector [X Y]^T as

Next, we define $G_x^g$ and $G_y^g$ to be the gth guesses for $E_x^{n+1}$ and $E_y^{n+1}$. 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.

4. 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 +x-directed 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: β ₁ = 0.69617, β ₂ = 0.40794, β ₃ = 0.89748, ω ₁ = 2.7537 × 10¹⁶ rad/s, ω ₂ = 1.6205 × 10¹⁶ rad/s, and ω ₃ = 1.90342 × 10¹⁴ rad/s. The third-order electric susceptibility scalar is ${\chi }_0^{\left(3\right)}$ = 1.89 × 10^-22 m²/V². The relative strengths of the Kerr and Raman polarizations are given by the parameter α = 0.7. The Raman polarization parameters are τ ₁ = 12.2 fs and τ ₂ = 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^-18s. 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

 

Fig. 1. GVADE simulation results for temporal soliton propagation in a dispersive nonlinear material. These results reproduce those of [5].

Download Full Size | PDF

where H ₀ = 4.77×10⁷ A/m, ω_c = 4.35×10¹⁵ rad/s (λ ₀ = 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 .

Download Full Size | PDF

5. 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.